Genetic and Morphologic Variation in a Potential Mosquito Biocontrol Agent, Hydrochara Affinis (Coleoptera: Hydrophilidae)

Hydrochara affinis (Coleoptera: Hydrophilidae), a water scavenger beetle, was recently identified as a natural and effective agent for biological mosquito control; it was reported to exhibit high rates of mosquito larvae predation. However, maintaining the quality (i.e., natural ecological attributes, such as genetic variation) of laboratory-reared populations is essential for ensuring the long-term success of biological control programs. Accordingly, here, we aimed to use mitochondrial Cytochrome c oxidase subunit I (COI) sequences to document the genetic diversity, population structure, and phylogenetic position of natural and lab-reared H. affinis populations in South Korea and use geometric morphometric analysis to investigate the populations’ morphological divergence. The natural H. affinis populations possessed high genetic diversity and numerous COI haplotypes, suggesting that these populations were healthy and could be directly applied to mosquito habitats without alterations to their natural genetic attributes. The lab-reared populations also possessed high genetic diversity and, thus, the potential for high adaptive capacity to new environments. Although no distinct population genetic structures were observed, quantitative variation was observed in the body shape of both the natural and lab-reared populations. The high levels of genetic and morphologic variation observed in the H. affinis populations examined here indicate the species’ favorable conservation status, genetic diversity, adaptive capacity, and, thus, “suitability” for field application as an effective mosquito control agent.


Introduction
The Hydrophilidae (i.e., water scavenger beetles) is one of the largest families of aquatic insects [1], and the larvae of some hydrophilid species have been reported as effective predators of mosquito larvae [2]. Hydrochara affinis (Hydrophilidae: Hydrophilinae), for example, was recently described as an effective biological control agent for mosquitos in South Korea, and both the efficacy and preference of H. affinis for mosquito (Culex pipiens molestus and Ochlerotatus togoi) larvae were investigated under laboratory conditions [3]. Third instar H. affinis consumed more mosquito larvae than other Sustainability 2020, 12, 5481 2 of 12 developmental stages, and prey handling time decreased dramatically with the progression of H. affinis development (i.e., instar number). The maximum numbers of consumed mosquito larvae, which were estimated using predation curves, were 926 and 304 larvae per day for C. pipiens molestus and O. togoi, respectively, and H. affinis exhibited a preference for C. pipiens molestus when both C. pipiens molestus and O. togoi were offered simultaneously [3].
For successful application, the introduction of predators into natural fields for biological pest control must consider several details, including predator efficiency and longevity and the level of expected disturbances to naturally existing populations. With regard to mosquito control, however, biocontrol efforts have focused mainly on the larvicidal efficiency of predator species [2][3][4], and the genetic compositions of the target populations have rarely been considered. Rudimentary population genetic information has been incorporated into biological control programs, before release, in terms of safety and efficacy [5][6][7]. For example, genetic analysis of the root-crown mining weevil (Ceutorhynchus scrobicollis), which was being considered as a biological control agent for garlic mustard (Alliaria petiolata), provided an estimate of the number of individuals required from a focal area to obtain the high level of genetic diversity of control agent [7], and analysis of cultured (quarantined) and naturalized (after release) populations of the mirid bug Eccritotarsus catarinensis, a water hyacinth biological control agent, revealed that genetic bottleneck events did not affect the genetic diversity of introduced biocontrol populations [8]. Thus, population genetic data can be useful for evaluating the distribution and differentiation of candidate biological control agent populations [7,8]. Mass rearing systems must be established and maintained for the successful and repeated release of biocontrol agents into natural habitats, and the quality (i.e., natural ecological attributes, such as genetic variation) of laboratory-reared populations is essential for ensuring the long-term success of biological control programs. However, systems for the mass rearing of insects in the lab are likely to incur high rates of inbreeding and, thus, reductions in fitness (e.g., fecundity), growth rate, and viability [9][10][11]. Indeed, such issues have been noted since the advent of biocontrol research in the 1970s [9], and biocontrol efforts should also aim to prevent any negative effects of agent release on the genetic attributes of natural populations.
Morphological variation indicates the ecological characteristics of populations as well as their capacity to adapt to different local environments [12][13][14][15][16]; as such, geometric morphometric analysis has been used to investigate the roles of both environmental factors and genetic components in shaping morphological variation [17][18][19][20]. Indeed, differences in the morphologies of populations under different environmental conditions can indicate the capacities of these populations to adapt to new environments [21,22]. The introduction of H. affinis from mass laboratory production to various mosquito habitats will entail novel selection pressures, owing to differences in environmental conditions. Thus, the differences in the morphology of the natural and laboratory populations might indicate whether those variations are the result of genetic divergence or phenotypic plasticity [23]. Furthermore, in the case of long-term biological control, geometric morphometric variations in the body shape of populations can also provide information about the adaptive capacity of natural populations following the introduction of lab-reared population.
Though the efficacy of H. affinis as a mosquito larvae predator has been described previously [3,24], the genetic (e.g., genetic diversity and structure) and phylogenetic characteristics of H. affinis populations in South Korea are yet to be reported. In particular, the genetic characteristics of H. affinis should be investigated to ensure success of the long-term biological control procedure including field application. Accordingly, the aims of the present study were to use mitochondrial Cytochrome c oxidase subunit I (COI) sequences to document the genetic diversity, population structure, and phylogenetic position of natural and lab-reared H. affinis populations in South Korea and use geometric morphometric analysis to investigate the populations' morphological divergence. These results will provide valuable information about the current conservational status of natural H. affinis populations and will elucidate the genetic and adaptive properties of "mass laboratory production population", needed to establish a successful mosquito biological control system.

Sample Collection and COI Sequencing
In total, 82 H. affinis specimens (four larvae and 78 adults) were collected from 21 sampling locations in South Korea ( Figure 1) and from one laboratory rearing (breeding) population at the Deokso Korea University Research farm. Specimens (N = 15) of an endemic and co-occurring congener, Hydrochara libera [25,26], were also collected from one location, to investigate the species' phylogenetic relationship to H. affinis. All specimens were stored in 96% ethanol for genetic analyses.
Deokso Korea University Research farm. Specimens (N = 15) of an endemic and co-occurring congener, Hydrochara libera [25,26], were also collected from one location, to investigate the species' phylogenetic relationship to H. affinis. All specimens were stored in 96% ethanol for genetic analyses.
Genomic DNA was extracted from a femur of the adult specimens and from a whole leg of the larva specimens using the DNeasy Blood and Tissue Kit (Qiagen, USA). A fragment (771 bp) of the mitochondrial Cytochrome c oxidase I gene (COI) was then polymerase chain reaction (PCR)-amplified using specific forward (C1-J-2183: CAA CAT TTA TTT TGA TTT TTT GG) and reverse (TL2-N-3014: TCC AAT GCA CTA ATC TGC CAT ATT A) primers [27] and the following conditions: initial denaturation of 1 min at 94 °C; 35 cycles of 30 s at 94 °C, 30 s at 50 °C, and 1-2 min at 72 °C; and a final extension step of 7 min at 72 °C. After electrophoresis on 1.5% agarose gels, the amplification was verified visually using UV light. The verified PCR products were purified enzymatically using Exonuclease I and Shrimp Alkaline Phosphatase (New England BioLabs, USA) and then sequenced by Macrogen INC Sequencing (Korea) on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems, USA). The COI sequences of H. libera and H. affinis were deposited in GenBank under accession numbers MT506869-MT506883 and MT549695-MT549757, respectively.

Phylogenetic Reconstruction and Haplotype Network Analysis
Phylogenetic tree reconstruction was performed using 110 COI sequences, which included 82 and 15 newly obtained sequences for H. affinis and H. libera, respectively. Additional sequences (N = 13) for Hydrochara species (N = 11), including H. affinis (KF128915) and H. libera (KY554236), and two outgroup species (N = 2), Sternolophus rufipes (KF128913) and S. marginicollis (KC935325), were downloaded from the GenBank database. Sequence alignment was performed using the Clustal-W multiple sequences alignment package [29], implemented in BioEdit 7.1.9 [30], and MEGA7.0.14 [31] was used to estimate P-distance, to evaluate the overall mean distance and the evolutionary divergence between specific sequences. For phylogenetic analysis, Bayesian phylogenetic inference was performed using MrBayes 3.2 [32], with GTR+G as the best-fit substitution model, which was determined using AICc in jModeltest 2.1.7 [33]. The MCMC run involved 5,000,000 generations with Genomic DNA was extracted from a femur of the adult specimens and from a whole leg of the larva specimens using the DNeasy Blood and Tissue Kit (Qiagen, Germanton, NC, USA). A fragment (771 bp) of the mitochondrial Cytochrome c oxidase I gene (COI) was then polymerase chain reaction (PCR)-amplified using specific forward (C1-J-2183: CAA CAT TTA TTT TGA TTT TTT GG) and reverse (TL2-N-3014: TCC AAT GCA CTA ATC TGC CAT ATT A) primers [28] and the following conditions: initial denaturation of 1 min at 94 • C; 35 cycles of 30 s at 94 • C, 30 s at 50 • C, and 1-2 min at 72 • C; and a final extension step of 7 min at 72 • C. After electrophoresis on 1.5% agarose gels, the amplification was verified visually using UV light. The verified PCR products were purified enzymatically using Exonuclease I and Shrimp Alkaline Phosphatase (New England BioLabs, Ipswich, MA, USA) and then sequenced by Macrogen INC Sequencing (Korea) on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). The COI sequences of H. libera and H. affinis were deposited in GenBank under accession numbers MT506869-MT506883 and MT549695-MT549757, respectively.

Phylogenetic Reconstruction and Haplotype Network Analysis
Phylogenetic tree reconstruction was performed using 110 COI sequences, which included 82 and 15 newly obtained sequences for H. affinis and H. libera, respectively. Additional sequences (N = 13) for Hydrochara species (N = 11), including H. affinis (KF128915) and H. libera (KY554236), and two outgroup species (N = 2), Sternolophus rufipes (KF128913) and S. marginicollis (KC935325), were downloaded from the GenBank database. Sequence alignment was performed using the Clustal-W multiple sequences alignment package [29], implemented in BioEdit 7.1.9 [30], and MEGA7.0.14 [31] was used to estimate Sustainability 2020, 12, 5481 4 of 12 P-distance, to evaluate the overall mean distance and the evolutionary divergence between specific sequences. For phylogenetic analysis, Bayesian phylogenetic inference was performed using MrBayes 3.2 [32], with GTR+G as the best-fit substitution model, which was determined using AICc in jModeltest 2.1.7 [33]. The MCMC run involved 5,000,000 generations with four chains, and the first 25% of the samples were discarded as burn-in. Maximum-likelihood (ML) analysis was performed using the same model with 1000 bootstrap replicates and was implemented using PhyML Web-Servers [34].
Haplotypes of the 82 H. affinis individuals were determined using Neighbor-Joining algorithms in DnaSP v5 [35], and a haplotype network was inferred using HAPSTAR v0.7 [36]. Genetic diversity indices, including haplotype diversity (h) and nucleotide diversity (π), were estimated using ARLEQUIN v3.5 [37], and Tajima's D and Fu's Fs statistics were calculated to test for neutrality. Haplotype richness was calculated for the five populations from which at least five specimens (N ≥ 5) had been collected using a rarefaction method in CONTRIB v1.02 [38], which corrects for unequal sample sizes among populations.

Geometric Morphometric Analyses
In total, 86 H. affinis specimens from lab-rearing and natural population, were used for morphometric analyses. The ventral side of each individual was independently photographed three times to reduce imaging errors using a Nikon D750 digital camera (Nikon Co., Japan), and the photographs were converted to .tps files using tpsUtil v.1.74 [39] to generate landmarks. Twenty-four anatomically homologous landmarks that are commonly used in beetles [40,41] were selected for morphometric analysis (Figure 2a), and the specimen data were digitized using tpsDig v. 2.30 [42]. The entire configuration dataset was regarded as object symmetry to remove configuration reflection [43,44]. Full Procrustes fit was implemented in MorphoJ v. 1.06d [45] to obtain shape variables and remove extraneous features, such as size, orientation, and location [43,45,46].  Procrustes ANOVA (analysis of variance) was performed to compute imaging errors from 258 raw coordinate configurations [43,47]; then, the three replicates of each specimen were combined. Canonical variate analysis (CVA) was performed, with 1000 round permutation tests, for the four populations Sustainability 2020, 12, 5481 5 of 12 (ST1, 3, 19, and STL) from which more than five individuals had been collected, and populations ST1 and ST2 were considered one population (N = 14), owing to their geographic proximity (~4 km).

Haplotype Analysis
The 771-bp COI sequences generated from the 82 H. affinis specimens yielded 50 variable sites, including 32 parsimony-informative sites and 18 singletons. The haplotype diversity (h) values ranged from 0.833 to 1.000, and the haplotype richness (hr) of populations ST1, ST2, ST3, ST19, ST20, and STL ranged from 2.770 to 4.000 (Table 1). Among the six locations from which five or more specimens were collected, ST2 (N = 8) yielded the greatest haplotype richness (hr) and nucleotide diversity (π). All the Tajima's D and Fu's Fs values were negative, except the Tajima's D values of the STL population.
In total, 46 distinct COI haplotypes were identified (Table 1, Figure 3), with the greatest number (NH = 12) identified in the laboratory population (STL) ( Table 1). Most (N = 33, 71.7%) of the haplotypes were obtained from single specimens (i.e., singletons), and only three of the distinct haplotypes, namely H24 (14.6%), H5 (13.4%), and H1 (6.1%) accounted for more than 5% of the specimens (Figure 3). Furthermore, the three main haplotypes (H1, H5, and H24), which were separated by less than three mutational steps, were obtained from 28 (34.1%) of the 82 specimens and were detected in 12 (54.5%) out of the 22 populations. In addition, the second most common haplotype (H5) was located at the center of the haplotype network, and the pattern of "star-like" haplotype network was identified from H. affinis in South Korea (Figure 3).

Geometric Morphometric Analysis
According to Procrustes ANOVA, the mean square for the individual level exceeded the imaging error level. Measurement errors only contributed to 2.91% of the individual variation and were, therefore, ignored (Table S1). The CVA revealed differences in the shape of the four included

Geometric Morphometric Analysis
According to Procrustes ANOVA, the mean square for the individual level exceeded the imaging error level. Measurement errors only contributed to 2.91% of the individual variation and were, therefore, ignored (Table S1). The CVA revealed differences in the shape of the four included populations (Figure 2b), but significant differences (p < 0.05) were only observed between the shapes of the ST3 and STL populations and between the ST3 and ST19 populations (Figure 2b, Table 2).

Phylogenetic Status of H. affinis and H. libera in the Genus Hydrochara
The overall mean distance (d) between the 11 Hydrochara species was 0.0377 (SE = 0.0042), and the mean intraspecific genetic diversity for H. affinis and H. libera was 0.0049 (SE = 0.0009) and 0.0051 (SE = 0.0012), respectively.
The topology of the BI tree was congruent with that of the ML tree ( Figure 4). Indeed, the monophyly of genus Hydrochara was strongly supported, with support values of 100 in both trees, and all the H. affinis sequences formed a clear monophyly with the high support values (≥98), as did the sequences from H. libera. No clusters were observed among the sequences from specific localities or regions.

Phylogenetic Status of H. affinis and H. libera in the Genus Hydrochara
The overall mean distance (d) between the 11 Hydrochara species was 0.0377 (SE = 0.0042), and the mean intraspecific genetic diversity for H. affinis and H. libera was 0.0049 (SE = 0.0009) and 0.0051 (SE = 0.0012), respectively.
The topology of the BI tree was congruent with that of the ML tree ( Figure 4). Indeed, the monophyly of genus Hydrochara was strongly supported, with support values of 100 in both trees, and all the H. affinis sequences formed a clear monophyly with the high support values (≥98), as did the sequences from H. libera. No clusters were observed among the sequences from specific localities or regions. Though only some of the relationships with Hydrochara were fully supported, H. affinis was the most closely related to H. flavipes, with support values of ≥97 in both phylogenetic trees and a mean Though only some of the relationships with Hydrochara were fully supported, H. affinis was the most closely related to H. flavipes, with support values of ≥97 in both phylogenetic trees and a mean genetic distance of 0.0508 (SE = 0.0087). Meanwhile, H. affinis and H. libera were more distantly related, with a mean distance of 0.0949 (SE = 0.0115).

Discussion
The high levels of genetic and morphologic variation observed in the H. affinis populations examined in the present study indicate the species' favorable conservation status, genetic diversity, adaptive capacity, and, thus, "suitability" for field application as an effective mosquito control agent. The high level of genetic diversity and richness of the natural H. affinis populations in South Korea (Figure 3, Table 1) indicate the species' large effective population sizes (N e ) and relatively low risk of population reduction or extinction in the sampled areas. Furthermore, the absence of any distinct population structures implies that specimens from any site can be used for direct application in areas where mosquitos and H. affinis already co-occur without disturbing the genetic attributes of the local H. affinis populations. This is important because the invasion of native populations by non-native biological control agents during biological control efforts have been reported [48], and the introduction of new genetic components can disturb the original genetic attributes of the native populations [48]. Indeed, it would be nearly impossible to remove the introduced H. affinis lineages once they have been introduced, especially since the introduced larvae of H. affinis would persist and, after maturing, would interbreed with individuals from the natural populations, as has been reported for larvivorous fish and amphibians that have been introduced for mosquito control [49][50][51].
It is especially interesting that the laboratory rearing population investigated in the present study also maintained a relatively high level of genetic diversity (Table 1), which might indicate the population's high capacity to adapt to various new environments. Indeed, the lab population's genetic diversity (i.e., h and HR) was similar to that of the natural populations (Table 1). This is important because losses of genetic diversity are common among mass laboratory rearing populations of insects, and it is well known that inbreeding depression causes a variety of negative effects on fitness (e.g., low growth rate, reduced reproductive success, high mortality) [11,52,53]. Genetic diversity (i.e., intraspecific variation) is also fundamental to a species' ability to adapt to environmental changes [7,54], and the capacity of biological control agents to adapt to new environmental conditions is critical for successful mosquito control, in particular, since mosquitos occur in a wide variety of habitats, including ponds, marshes, swamps, and variously polluted habitats. Indeed, the adaptive capacity of control agents also influences their survival and longevity, which thereby contributes to their efficacy.
The captive population examined in the present study was established three years prior (in 2016) and had been maintained through the continuous introduction of H. affinis specimens from natural habitats in South Korea. Though the source populations of the laboratory stock were not tracked, the population was observed to share haplotypes with several natural populations. Indeed, the haplotypes that were most common among the natural populations (H1, H5, and H24) also occurred at relatively high frequencies in the laboratory population (52.6% for all three), thereby indicating that the laboratory population was established from natural populations. The laboratory population's relatively high intraspecific genetic variation and frequency of unique haplotypes (47.4%) also suggested that the population's genetic diversity had not decreased significantly since its establishment. Thus, the results of the present study support the notion that the mixture of individuals from genetically distinct populations could be important for establishing and maintaining genetic variation in H. affinis rearing populations. The results also suggest that regular restocking from natural populations could help maintain the genetic diversity of laboratory rearing populations by facilitating outbreeding.
The genetic diversity presented here, for both natural and captive populations of H. affinis (Table 1), could be used as a standard for the level of diversity needed for H. affinis mosquito control systems established in South Korea. The study also provides a basis for evaluating the genetic diversity of rearing populations before their release to mosquito control sites and for detecting reductions in genetic diversity that could inhibit the efficacy of H. affinis.
In addition to measuring genetic diversity and richness, geometric morphometric analysis can also provide insight into the capacity of populations to adapt to new habitats. In the present study, the three natural populations that were examined exhibited clear morphologic differences and were even more distinct from the laboratory population (Figure 2b). Body shape is considered an important ecomorphological phenotype that reflects adaptation to complex aquatic habitats [55], and the shape-environment association with high plasticity has been widely reported [44]. Because morphological variation is affected by both genetic and environmental factors [56] and because the 9 of 12 four H. affinis populations investigated in the present study exhibited low levels of genetic divergence (i.e., non-significant F ST values; Table 2), it is possible that the body-shape variation in H. affinis represents a phenotypically plastic response to local environmental conditions (e.g., food, temperature, water flow, microhabitat), rather than genetic adaptation. The clear difference between the morphology of the natural populations and that of the laboratory population, with more controlled environmental conditions, further supports this conclusion. More importantly, however, these results demonstrate that the laboratory rearing population, which was established using individuals from multiple natural populations, has not lost its capacity to adapt to new environments. Indeed, it is possible that genetic diversity and body-shape variation could be used to evaluate the "suitability" of lab-reared populations as control agents before being released into natural mosquito habitats.
Though H. affinis has been investigated with regard to its use as a mosquito control agent, the present study is the first, to our knowledge, to report the phylogenetic position of the species using molecular phylogenetic analysis. The genus Hydrochara includes 22 species (excluding H. major, which was recently reclassified into the genus Brownephilus) [57][58][59][60], and three species (H. affinis, H. libera, and H. vincina) have been reported to occur on the Korean Peninsula [61][62][63][64][65], with H. affinis and H. libera sharing similar habitats and distributions [25,26]. However, only H. libera has been included in a molecular phylogenetic study [66]. Furthermore, several key morphological characters (e.g., maxillary palpus color, appendage (leg) color, and length of the posterior spine on the sternal keel) have been used to distinguish H. affinis and H. libera [57]. However, the genetic divergence of the two species has yet to be reported. The present study revealed that H. flavipes is more closely related to H. affinis than to H. libera (Figure 4). However, despite the phylogenetic and morphologic differences between H. affinis and H. libera, the two species share similar distributions, ecological niches, habitats, and life cycles [25], which indicates that H. libera could represent another mosquito control agent that is yet to be discovered. The present study represents the first genetic and morphologic assessment of an aquatic insect, H. affinis, as a biological mosquito control agent in South Korea and provides valuable information for establishing efficient systems for the laboratory production of high-quality biological control agents in general.