Morphology versus DNA barcoding: two sides of the same coin. A case study of Ceutorhynchus erysimi and C. contractus identification

Genotyping of 2 well‐known weevil species from the genus Ceutorhynchus (Coleoptera: Curculionidae) distributed in west Palearctic, C. erysimi and C. contractus, revealed phenotype versus genotype inconsistencies in a set of 56 specimens (25 C. erysimi and 31 C. contractus) collected from 25 locations in Serbia and Montenegro. An analysis of the mitochondrial cytochrome oxidase subunit I gene (COI), widely used as a barcoding region, and a nuclear gene, elongation factor‐1α (EF‐1α), revealed stable genetic divergence among these species. The average uncorrected pairwise distances for the COI and EF‐1α genes were 3.8%, and 1.3%, respectively, indicating 2 genetically well‐segregated species. However, the genetic data were not congruent with the phenotypic characteristics of the studied specimens. In the first place, C. erysimi genotypes were attached to specimens with phenotypic characteristics of C. contractus. Species‐specific PCR‐RFLP assays for the barcoding gene COI were applied for the molecular identification of 101 additional specimens of both morphospecies (33 C. erysimi and 68 C. contractus) and were found to confirm this incongruity. The discrepancy between the genetic and morphological data raises the question of the accuracy of using a barcoding approach, as it may result in misleading conclusions about the taxonomic position of the studied organism. Additionally, the typological species concept shows considerable weakness when genetic data are not supported with phenotypic characteristics as in case of asymmetric introgression, which may cause certain problems, especially in applied studies such as biological control programs in which the biological properties of the studied organisms are the main focus.


Introduction
Biodiversity protection, preservation and inventorization are leading concepts in modern biology; however, species identification and precise recognition still remain the key problems in basic and applied studies. The rapid development of molecular techniques during the 1990s, accompanied by the DNA barcoding initiative launched by Hebert et al. (2003), lent insight into the potential of molecular methods to provide better taxonomic resolution in species recognition efforts (Hebert et al., 2004), especially as the decline of traditional taxonomic knowledge has become more than obvious (Taylor & Harris, 2012). However, it is indisputable that during the last 250 years, traditional taxonomy accumulated knowledge that has been incorporated into every aspect of current research into the Earth's biodiversity. Without doubt the barcoding initiative became an extraordinary useful tool for fast species identification and inventorying organisms (Miller, 2007), but the epistemological value of this method has been questioned by many authors (DeSalle et al., 2005;Brower, 2006;Meier et al., 2006;Song et al., 2008). In general, from a taxonomical point of view, the criticism of the barcoding approach was rightly directed (Goldstein & DeSalle, 2011). In addition, barcoding practitioners have made no efforts to stabilize any nomenclature issues in inventoried taxa, despite the fact that the current taxonomy is based exclusively on the typological species concept, following the International Codes of Nomenclature (Zoological, Botanical, Bacteria, etc.). The typological species concept is of essential importance for linking traditional taxonomy with novel approaches, such as DNA barcoding. The primary determination of any organism at the species level is always based on traditional taxonomic methods; if this link breaks, than fixation of inadequate DNA barcoding sequences for particular organisms in public databases offers no benefit, only an unnecessary mess (Peterson et al., 2007). Therefore, traditional taxonomy and DNA barcoding should not be treated as antagonistic methods, but rather as complementary approaches for resolving the phylogeny, taxonomy and nomenclature of living organisms (Packer et al., 2009;Taylor & Harris, 2012).
The weevil genus Ceutorhynchus Germar, 1824 (Coleoptera: Curculionidae) is composed of approximately 380 species, widespread all over the Holarctic region (Colonnelli, 2004). Only several species are known to occur outside of this area, in Central and South Africa (Colonnelli, 2006;Korotyaev, 2008). Nearly all known species of Ceutorhynchus utilize members of the family Brassicaceae and the closely related Resedaceae as host plants. Some of the species are known as serious pests on economically important cruciferous crops (Dosdall et al., 2011;Vaitelytė et al., 2013). On the other hand, several Ceutorhynchus species are treated as beneficial organisms and have been considered in various biological control programs as potential agents for the control of invasive weeds from the family Brassicaceae (Fumanal et al., 2004;Hinz et al., 2008;Gerber et al., 2009;Gerber et al., 2012;Hinz & Diaconu, 2015).
In general, the biology of the Ceutorhynchus weevils is poorly understood, with biological data mostly presented as lists of cruciferous plants that have been noted as potential hosts. The phylogenetic relationships within the genus Ceutorhynchus are also poorly studied, without hardly any comprehensive work elaborating this topic in the literature, with the exception of a few contributions (Laffin et al., 2005;Hinz et al., 2008;Rauth et al., 2011). This is most likely a reflection of the obvious fact that taxonomic issues related to many Ceutorhynchus species are still unresolved, given that in many species groups only a few characters can delimit closely related species within this very species-rich genus. Taxonomical and nomenclatural issues thus represent the main constraints in studying this important weevil group. Korotyaev (2008) noted that the comprehensive classification of the genus Ceutorhynchus is a problematic task and proposed, as a more efficient approach, to gradually distinguish subgenera for the most clearly defined groups of species. Because both the biology and genetics/phylogeny of Ceutorhynchus are understudied, work on the molecular identification of the different species is necessary to clarify their taxonomic positions and relationships. This is especially relevant for applied studies where precise species identification is a precondition for any laboratory testing.
Our study on the molecular identification of Ceutorhynchus species (primarily using barcoding) was initiated by the need for more precise determination of their field host associations and level of polyphagy in the context of studying their potential use as biocontrol agents. The concept of the study was to match adult specimens (identified morphologically and molecularly) that was found to be feeding on different plants with larvae (identified molecularly) developing within the suspected plants, host and niche. This approach should have provided a direct link between morphological and molecular identification, host association and behavior related to the niche for oviposition and larval development. The data presented in this paper are related to the evaluation of the taxonomic status of these 2 well-known polyphagous Ceutorhynchus weevil species C. erysimi (Fabricius, 1787) and C. contractus (Marsham, 1802) using both classical and molecular techniques. Both species are easily recognizable morphologically and are ascribed in the literature as polyphagous and very common (Hoffmann, 1954;Dieckmann, 1972). Ceutorhynchus erysimi is mainly associated with Capsella bursa-pastoris for larval development and is less frequently found on other crucifers (Colonnelli, 2004), whereas C. contractus is proposed to utilize a broad plant host range, including different Brassicaceae, Resedaceae, and Capparaceae (Hoffmann, 1954;Colonnelli, 2004). Both species possess very broad and overlapping distributions in the west Palearctic region, with the northernmost distribution in the subarctic zone and a southeastern distribution expanding to the Middle East and northern Africa. Because of their wide distribution and abundance in the western Palearctic, both species are likely to be the most common species in private and public weevil collections.

Insect sampling
This study included weevils collected in the territory of Serbia and one location in Montenegro. The weevils were collected using 2 methods: (a) entomological net sweeping and (b) collecting mature L3 larvae from their host plant and subsequently rearing them to the adult stage. For sweeping, we targeted particular plant species of interest, avoiding mixed plant populations with 2 or more cruciferous plants. The collected weevils were stored just after collection in 96% ethanol at 5-7°C until DNA extraction. The larvae were collected from various cruciferous plants, but mainly from the root-crown for larvae of C. erysimi, and from basal parts of the stem leaf petioles for C. contractus. Collected L3 larvae were released into plastic containers (6 × 5 cm) filled with 2-cm of sifted loamy soil and closed with gauze leads. The containers were sprayed daily with water to provide humidity to the soil. After emergence, adults were fed with leaves of various crucifers for the next 5 d until they were fully sclerotized and then stored in 96% ethanol until DNA extraction. After DNA extraction, weevils were prepared as dry voucher specimens for taxonomic studies. Several surveys were performed to collect insect material between April and June of 2010-2014 ( Fig. 1, Table S1).

Morphological analysis
Following the descriptions in Hoffmann (1954) and Dieckmann (1972), we used 2 phenotypic characters to clearly differentiate C. erysimi from C. contractus: the body length (rostrum excluded) and the structure and color of elytra.
Polymerase chain reactions (PCR) for both amplifications were performed in a 20-μL final volume containing Kapa Biosystems High Yield Reaction Buffer A with 1.5 mmol/L MgCl 2 (1×), an additional 2.25 mmol/L of MgCl 2 , 0.6 mmol/L of each dNTP, 0.5 μm of each primer, 1 U of KAPATaq DNA polymerase (Kapa Biosystems, Inc., Woburn, MA, USA) and 1 μL of DNA extract. PCR cycles were carried out in a Mastercycler ep gradient S (Eppendorf, Hamburg, Germany) applying the following thermal steps for amplification of the COI barcoding region: 95°C for 5 min (initial denaturation), 35 cycles at 94°C for 1 min, 54°C for 1 min (annealing), 72°C for 1.5 min and a final extension at 72°C for 7 min. The thermal protocol for amplification of the EF-1α gene was similar, with the exception of an annealing temperature of 50°C and the use of 40 cycles.

Phylogenetic analysis of the barcoding mtCOI gene and the nuclear EF-1α gene
Relationships resulting from intrinsic population-level processes (e.g., persistence of ancestral haplotypes, multifurcations, or recombination) are better visualized in reticulated graphs or networks (Posada & Crandall, 2001) than on evolutionary gene trees, although both may be informative at intraspecific levels in closely related species. Thus, gene genealogies for mtCOI and EF-1α were inferred using 2 approaches for haplotype network construction. Median-joining networks (Bandelt et al., 1999) were calculated with the program NETWORK version 4.6.1.2 (http://www.fluxus-engineering.com), keeping the parameter ε = 0 and the postprocessing option for maximum parsimony (MP) calculation. This method starts with minimum spanning trees combined within a single network, and then median vectors (consensus sequences) are added to reduce tree length. Such vectors can be interpreted as possibly extant unsampled sequences or extinct  Table S1. ancestral sequences (Bandelt et al., 1999). In addition, gene genealogies of the mtCOI and EF-1α genes were also inferred using TCS, version 1.21 (Clement et al., 2000), to infer haplotype networks using statistical parsimony (Templeton et al., 1992) with a confidence limit of 95%.

Virtual restriction analysis and gel plotting of the barcoding region of the mtCOI gene
The obtained sequences of the COI barcoding region were edited using FinchTV v.1.4.0 (http://www. geospiza.com) and aligned using the ClustalW program integrated into the MEGA5 software (Tamura et al., 2011). In silico restriction analysis was performed for all registered COI haplotypes to identify suitable restriction enzymes that could be used for identification of C. erysimi and C. contractus. Sequences of all haplotypes of both species were aligned, and the presence of speciesspecific SNPs (single nucleotide polymorphisms) in the recognition sites for restriction endonucleases was determined using the pDRAW32 software (AcaClone Software, http://www.acaclone.com). The SspI endonuclease was selected based on its recognition sequence, which occurred twice within the sequences of the C. erysimi COI haplotypes and once within the C. contractus COI haplotypes, tentatively generating 3-fragment and 2-fragment profiles, respectively. Virtual SspI restriction patterns were separated on a 1.5% agarose gel generated in the pDRAW32 program and compared with the actual enzymatic RFLP (restriction fragment length polymorphism) patterns of amplicons obtained from specimens representing each COI haplotype.

Species-specific PCR-RFLP analysis of the barcoding region of the mtCOI gene
For the species-specific RFLP identification of C. erysimi and C. contractus via the barcoding region of the mtCOI gene, an additional set of 101 specimens, 33 C. erysimi, and 68 C. contractus, were collected from 5 and 7 locations, respectively (Table S2). Some of the specimens for this analysis were collected as L3 larva instars and then reared to adults under laboratory condition, whereas others were collected from particular plant populations by the sweeping net method. Taken together, the material for the morphological and genetic study was collected from 25 localities, of which 6 sites were syntopic for the collection of both species (Fig. 1).
The COI amplicons amplified with the LCO1490/ HCO2198 primer pair were 709 bp in length and were subjected to restriction digestion with the SspI endonuclease (Fermentas, Vilnius, Lithuania) to obtain species-specific patterns typical for the C. erysimi and C. contractus genotypes. All restriction analyses were performed according to the manufacturer's instructions. The restriction products were separated via the automated capillary electrophoresis QIAxcel Advanced System (Qiagen) using the Screening Gel Cartridge (Qiagen) with the following parameters: sample injection voltage 5 kV, sample injection time 8 sec, separation voltage 6 kV and separation time 320 sec. The QX alignment marker 15 bp/5 kb (Qiagen) was used for alignment of the obtained restriction fragments, and the QX DNA size marker FX174/HaeIII (Qiagen) was used for fragment size comparison.

Morphological study
The primary differences that clearly define the phenotypes of C. contractus and C. erysimi are the ground color of the vestiture from the dorsal side and the body size of the adult weevils. All examined C. erysimi specimens fit the description of Hoffmann (1954) and Deickmann (1972); from a dorsal perspective, the elytra were always shining metallic blue or metallic green to different extents, and the average body length (rostrum excluded) was 2.3 ± 0.2 mm (range 1.8-2.7, n = 58). Pronotum dorsally black with more or less bronze shine, densely punctured, 0.6× as long as wide, with median furrow more or less expressed. Elytra with intervals 2-2.4 times as wide as striae; striae slightly convex with more or less expressed transverse wrinkles; preapical tubercles on intervals 5-8, pronounced, shining black. Aedeagus rectangular, truncate apically. In all C. contractus specimens, the elytra shined black dorsally, with some specimens showing very discrete dark black/bluish reflections. Pronotum dorsally black sometimes with discrete bronze shine, densely punctured, 0.9× as long as wide, without median furrow. Elytra with intervals nearly twice times as wide as striae; striae slightly convex with less expressed transverse wrinkles; preapical tubercles on intervals 5-8, less pronounced, shining black. Aedeagus rectangular, truncate, sometimes slightly depressed along the apical margin. The average body length for C. contractus was 1.7 ± 0.2 (range 1.2-2.0, n = 99). Body length was significantly different between C. erysimi and C. contractus (ANOVA, F 1, 56 = 20.6, P < 0.01). Thus, following traditional taxonomical methods we considered 2 possible phenotypes of the specimens in our study: the black and small specimens were considered C. contractus, whereas larger and green/blue-colored specimens were considered C. erysimi.

Mitochondrial COI analyses
The final 633-bp trimmed alignment of mtCOI yielded a total of 46 (7.2%) polymorphic nucleotide sites, of which 29 were parsimony-informative. The maximum ingroup genetic distance between sequences obtained from C. erysimi and C. contractus was 4.3% (uncorrected), whit an average of 3.8% (Table S3). The median-joining network analysis revealed a total of 22 haplotypes, which, according to their substantial genetic divergence, formed 2 distinct clusters (Fig. 2). The C. contractus cluster contained a total of 6 haplotypes and all specimens possessed the typical phenotype characteristic for this species, that is, small in size and black in color. In the second cluster, composed of 16 haplotypes (hereafter referred to as the C. erysimi genotype group), 7 specimens expressed the typical C. contractus phenotype, whereas 25 other specimens were phenotypically C. erysimi. Among the sequences within the C. erysimi genotype group, an intraspecific genetic divergence of 0.4% (range: 0.2%-0.8%) was recorded, whereas the mean divergence among the 6 haplotypes of the C. contractus genotype group was 0.7% (range: 0.2%-1.1%) (Table S3). In addition, the statistical parsimony network did not connect the C. erysimi and C. contractus genotype groups at the 95% parsimony connection limits (data not shown).

Nuclear EF-1α analyses
We reduced EF-1α sequences to unique genetic variants to perform the haplotype network analysis. From the 667-bp alignment of EF-1α sequences (the intron was removed from analysis), there were a total of 25 polymorphic nucleotides (3.7%), of which 10 were parsimony-informative, with a maximum ingroup genetic distance of 1.6% (uncorrected), and an average of 1.3% (Table S4). A total of 21 EF-1α haplotypes were revealed using median-joining network analysis; these were clustered in 2 haplotype groups, referred to as the EF-e genotypes for C. erysimi (8 haplotypes) and the EF-c genotypes for C. contractus (13 haplotypes). A mean divergence of 0.8% was estimated from all obtained haplotypes from both species. Among sequences within each genotype groups, an intraspecific genetic divergence of 0.4% (range 0.1%-0.6%) was recorded (Table S4). With regard to the phenotype of the analyzed specimens, 1 specimen (voucher code L333, Table S1) that was phenotypically identified as C. contractus carried an EF-1α haplotype associated with the C. erysimi EF-e4 genotype group as well as the C. erysimi COI haplotype ery9 (Fig. 2, Table S1).

PCR-RFLP analysis of the barcoding region of the mtCOI gene
A total of 101 individuals were assessed by a PCR-RFLP assay to determine the RFLP profile of the tested specimens. We tested 33 specimens that exhibited the phenotype of C. erysimi and 68 specimens with a typical C. contractus phenotype. The COI amplicons obtained with the LCO1490/HCO2198 primers (709 bp) were subjected to restriction digestion with the SspI endonuclease to obtain information on the restriction profiles that should be characteristic for each species. The SspI endonuclease cuts amplicons of C. erysimi at 2 positions, giving fragments of 385, 235, and 89 bp, whereas the C. contractus genotype possesses 1 restriction site, giving fragments of 474 and 235 bp (Fig. 3G). The restriction profiles that delimited the C. erysimi and C. contractus genotypes were confirmed via virtual digestion using pDRAW32 software on all 22 haplotypes obtained from both genotype groups. Our RFLP analysis confirmed the presence of C. erysimi mitochondrial genotypes in 10 of 68 (14.7%) specimens that exhibited the typical phenotype of C. contractus (Fig. 3), whereas no C. contractus mitochondrial genotypes were found in 33 analyzed specimens with the C. erysimi phenotype.

Discussion
The haplotype network of mitochondrial COI and nuclear EF-1α sequences constructed with the phenotypically well-separated species C. erysimi and C. contractus from the territory of Serbia and Montenegro demonstrated an unusual genetic richness in both of the studied species. In our study, which covered a relatively small geographic area, the observed haplotype richness values for COI and EF-1α were 22 and 21 haplotypes, respectively. We recorded 6 haplotypes for the with mitochondrial COI gene of C. contractus of 31 analyzed specimens from 19 locations, whereas almost 3 times as many haplotypes were recorded for C. erysimi of 25 specimens collected from 12 locations. Surprisingly, more haplotypes were identified in the EF-1α gene for C. contractus (= 13) than C. erysimi (= 8), although the diversity of the EF-1α gene was expected to be lower for both species, given that nuclear EF-1α is a more conservative gene. Nevertheless, the main finding of our study is the discrepancy between phenotypic and genotypic content in a subset of the analyzed specimens. This discrepancy was recorded primarily for specimens with phenotypic characteristics typical of C. contractus; mitochondrial genotypes typical for C. contractus were not recorded in specimens with the typical phenotype of C. erysimi (Fig. 2B). Observed pattern of asymmetric introgression is very common especially in the events of local and introduced species hybridization (Currat et al., 2008;Excoffier et al., 2009) as well as behavioral and mating differences between 2 sympatric species (Patton & Smith, 1993;Gomes et al., 2009). In the case of these 2 well-studied Ceutorhynchus species, the latter is more plausible, having in mind their sympatric occurrence. However, further research is needed to determine exact behavioral and mating differences and as a consequence, the levels of asymmetric introgression. This discrepancy applies to the nuclear EF-1α gene, because we recorded 1 specimen with the C. contractus phenotype but an EF-1α haplotype clustering within the C. erysimi genotype group (Fig. 2C). Accordingly, ambiguity in the mitochondrial data strongly suggests more or less episodic hybridization events between C. erysimi and C. contractus. Genetic diversity in regard to the conservative EF-1α gene also supports the existence of hybridization events between the 2 weevil species, although it is likely that during such events, the surviving offspring subsequently underwent repeated assortative mating with one of the parental lineages. In this scenario, the accumulation of C. erysimi mitochondrial genotypes inside populations of C. contractus weevils should represent a consequence of recent or historical hybridization events with C. erysimi followed by an asymmetric introgression. A plausible explanation for the subsequent sorting of nuclear alleles according to morphotypes of the nuclear genetic content in the C. erysimi and C. contractus populations could be the gregarious behavior of the adult populations, especially those with less expressed polyphagous behavior. Thus, our study supports the idea that nuclear genotyping is more reliable method for defining the real taxonomic positions of specimens with ambiguous mitochondrial contents.
Initially, our intent was to understand the speciation process and enable precise species identification within the weevil genus Ceutorhynchus, keeping in mind that some of the species from this genus have been proposed as biological control agents for the control of invasive cruciferous weeds. The polyphagous behavior of C. erysimi and C. contractus deserves special attention because both Fig. 3 The mtCOI PCR-RFLP species-specific identification of Ceutorhynchus erysimi and C. contractus. (A-F) SspI RFLP profiles of the COI barcoding region (LCO1490/HCO2198) of 101 specimens of C. erysimi (33) and C. contractus (68) identified by phenotype. Restriction fragments were separated using the capillary electrophoresis system QIAxcel advanced (Qiagen). The reference controls for the restriction pattern comparison are C. erysimi (3280) and C. contractus (3299). M: QX DNA size marker FX174/HaeIII (Qiagen). Fragment sizes (bp) of the marker (1353, 1078, 872, 603, 310, 281, 271, 234, 194, 118 and 72) and the alignment marker QX 15 bp/5 kb (15 and 5000) are designated. Arrows indicate specimens with mismatched species designation according to phenotype and COI-RFLP genotype. (G) In silico gel electrophoresis of the COI sequence (delimited by the LCO1490 and HCO2198 primer annealing positions) of control specimens of C. erysimi (3280) and C. contractus (3299) virtually digested with SspI. The abbreviations "ery." and "con." are used for C. erysimi and C. contractus, respectively. species possess a characteristic niche for development: C. erysimi occupies the root crown for larval development, whereas C. contractus uses the basal part of steam leaf petioles. In addition, both species were easy to rear and adults were easily distinguishable. The first ambiguity recorded in our study occurred during our attempt to distinguish larvae of C. contractus collected from diverse cruciferous plants, sampled from the typical host niche for this species, that is, the basal part of leaf petioles. For this purpose, we designed a fast identification protocol for C. contractus larvae according to the RFLP profiles of the barcoding region using the SspI endonuclease for the restriction digest. The observed ambiguity in the RFLP profiles forced us to sequence the barcoding region of these specimens, which revealed COI sequences typical for C. erysimi. Furthermore, our study of adult specimens of both species confirmed the existence of genetic ambiguity between these phenotypically well-segregated species.
Even though the situation of C. erysimi and C. contractus could be treated as an isolated case, the results presented in this study open several important questions related to the "undeclared war" between researchers advocating traditional or barcoding approaches in taxonomy (Brower, 2006;Miller, 2007;Song et al., 2008;Boero, 2010;Goldstein & DeSalle, 2011;Boero & Bernardi, 2014). In any case, for specimens that exhibit the phenotype of C. contractus but share the mitochondrial genotype of C. erysimi, it is not possible to determine their taxonomic position because if the traditional approach is applied, these specimens would be treated as C. contractus sensu stricto. However, this would not be the case if the barcoding approach were applied, especially if species determination were performed at the larval stage, which is advocated as one of the best reasons for applying the barcoding approach for species recognition (Köhler, 2007;Wong & Hanner, 2008). Thus, to stabilize taxonomic status of species with documented introgression, the process should consider applying other relevant methods, such as analysis of additional genetic markers, which would provide more reliable resolution in defining species limits.
In the end, it is not a question of which approach is preferable or which method better fits the needs of modern science, or, more importantly, those of practitioners involved in applied research. In this study, we present a very specific situation where 2 different approaches, traditional taxonomy and barcoding, failed to agree. We believe that such situations are not rare, particularly not in groups with pronounced species radiation or cryptic speciation events.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site:

Table S1
List of Ceutorhynchus erysimi and C. contractus specimens collected for morphological and molecular analyses, sorted by DNA voucher code, species name according to phenotypic and genetic content, host plant affiliation, locality, and haplotype name. Table S2 List of Ceutorhynchus erysimi and C. contractus specimens used for PCR-RFLP genotyping, sorted by DNA voucher code, phenotypic identification, host plant affiliation, locality and species designation according to PCR-RFLP analysis of the COI gene. Table S3 Average mitochondrial DNA cytochrome oxidase subunit I (COI) divergence based on a pairwise analysis (p-distance method) of haplotypes of the Ceutorhynchus erysimi and C. contractus genotype groups. Table S4 Average nuclear DNA elongation factor-1α gene (EF-1α) divergence based on a pairwise analysis (p-distance method) of haplotypes of the Ceutorhynchus erysimi and C. contractus genotype groups.