DNA barcoding of Deltocephalus Burmeister leafhoppers (Cicadellidae, Deltocephalinae, Deltocephalini) in China

Abstract We investigated the feasibility of using the DNA barcode region in identifying Deltocephalus from China. Sequences of the barcode region of the mitochondrial COI gene were obtained for 98 specimens (Deltocephalus vulgaris – 88, Deltocephalus pulicaris – 5, Deltocephalus uncinatus – 5). The average genetic distances among morphological and geographical groups of D. vulgaris ranged from 0.9% to 6.3% and among the three species of Deltocephalus ranged from 16.4% to 21.9% without overlap, which effectively reveals the existence of a “DNA barcoding gap”. It is important to assess the status of these morphological variants and explore the genetic variation among Chinese populations of D. vulgaris because the status of this species has led to taxonomic confusion because specimens representing two distinct morphological variants based on the form of the aedeagus are often encountered at a single locality. Forty-five haplotypes (D. vulgaris – 36, D. pulicaris – 5, D. uncinatus – 4) were defined to perform the phylogenetic analyses; they revealed no distinct lineages corresponding either to the two morphotypes of D. vulgaris or to geographical populations. Thus, there is no evidence that these variants represent genetically distinct species.


Introduction
China contains threatened biodiversity hotspots, including one spanning the Palearctic and Oriental regions and containing a high level of species diversity (Lin et al. 2010). In these regions, accurate identification of extant species is of great significance, although the taxonomic expertise is limited. Traditionally, identification of most species has been based on morphology. However, the availability of inexpensive DNA sequencing technology now provides additional tools not only for routine species identification but also for testing the validity of morphologybased species concepts. DNA barcoding is a simple, effective tool, that can identify and delimit species, including some complex taxa, rapidly and accurately using a standard short DNA sequence of the mitochondrial cytochrome c oxidase I (COI) (Hebert et al. 2003(Hebert et al. , 2004bWard et al. 2005;Hajibabaei et al. 2006). This method has been widely recognized and accepted in molecular phylogenetic studies (Hebert et al. 2003). The COI-based identification system has achieved remarkable success discriminating species across numerous animal groups, including birds (Hebert et al. 2004b), fishes (Hubert et al. 2008), and the insect orders Lepidoptera (Hebert et al. 2004a;Hajibabaei et al. 2006;Yang et al. 2012;Ashfaq et al. 2013), Ephemeroptera (Ball et al. 2005), and Hymenoptera (Smith et al. 2008). But this technology has also failed to identify species accurately under certain circumstances. For example, in a study of 449 species of Diptera and using 1333 COI sequences, Meier et al. (2006) obtained an identification success rate below 70% due to extensive overlap in inter-and intraspecific genetic distances. Within the dipteran family Calliphoridae, Whitworth et al. (2007) found that only 60% of species tested could be identified reliably.
Deltocephalini feed on grasses and sedges and are diverse and abundant in grassland ecosystems. This group contains 73 genera and 613 species around the world. Deltocephalus, type genus of this tribe contains 62 species distributed in the Old World and New World. Some species of this genus can transmit pathogenic diseases to economically important plants and are important economic pests; therefore, tools are needed for their rapid and accurate identification. Four species are described from China, two of them transmit pathogenic diseases. Identification of leafhopper species in most genera now requires dissection and examination of the male genitalia. However, some taxonomically problematic species apparently exhibit substantial intraspecific variation in male genital structures, and this causes confusion among taxonomists. One such practical example is D. vulgaris, which has well-documented morphological differences in the shape of the aedeagus (Figs 2, 3). Dash and Viraktamath (1998) first reported morphological variation in this species when they reviewed the genus Deltocephalus from India. Webb and Viraktamath (2009) also reported two forms of the aedeagus despite many shared morphological features in the species. Zhang and Duan (2011) redescribed D. vulgaris with detailed drawings and photos, illustrating these obvious morphological differences.
Based on DNA barcoding of leafhoppers, 63 barcodes from 45 species in Japan (15 subfamilies and 37 genera without Deltocephalini) were analysed (Kamitani 2011). DNA barcodes from 546 adult specimens of leafhoppers, planthoppers and treehoppers (Hemiptera, Auchenorrhyncha) were obtained from Barrow Island and analysed (Gopurenko et al. 2013). Species determination of members in the genus Aphrodes (Hemiptera, Cicadellidae) based on vibrational signals, mitochondrial DNA and morphology were performed (Bluemel et al. 2014). A total of 1482 specimens based on DNA barcodes of Nearctic Auchenorrhyncha (Insecta, Hemiptera) were studied by Foottit et al. (2014). The boundaries of seven closely related species of the evacanthine leafhopper genus Bundera (Cicadellidae, Evacanthinae) based on DNA barcoding, morphology and hyperspectral reflectance profiling was investigated by Wang et al. (2016), and a revision of the genus Orosius (Cicadellidae, Deltocephalinae, Opsiini) based on morphological and DNA barcoding was undertaken by Fletcher et al. (2017). Although, DNA barcoding research has been applied to these groups of leafhoppers, until now, a few molecular data are available for Deltocephalus. Therefore, a better understanding of Deltocephalus, and particularly the variation of D. vulgaris based on molecular data, is urgently needed.
In this study, we studied 98 COI sequences from three species of Deltocephalus. DNA barcoding data were used to investigate genetic variation of Chinese populations of D. vulgaris and to determine whether the morphological variants previously identified in this species represent distinct lineages. Our specific aims were to test the feasibility of using DNA barcoding data for identification of species of Deltocephalus, to determine the levels of the genetic variation within D. vulgaris, and to preliminarily discuss its possible correlation with morphological variation and biogeographic patterns.

Taxon sampling
A total of 98 specimens of Deltocephalus (D. vulgaris -88, D. pulicaris -5, D. uncinatus -5) were collected with an insect sweep net in the daytime and by a light trap at night. Specimens were all collected directly into 95% or 100% ethanol and stored in -80 °C prior to study. The sample included D. vulgaris, D. uncinatus and D. pulicaris to facilitate comparison of inter-to intraspecific genetic variation in this group. Deltocephalus vulgaris specimens were divided into 11 groups based on their morphological differences and different geographical distributions in China (

Morphology
Morphological observations were made using an Olympus SZX10 stereoscopic microscope (Olympus Corporation, Tokyo, Japan). All photographs and drawings were modified with Adobe Photoshop CS.

DNA extraction, amplification and sequencing
Total genomic DNA was extracted from the whole abdomen of each leafhopper using the EasyPure Genomic DNA Kit (EE101; Transgen, Beijing, China) following the manufacturer's instructions with the following modifications: abdomen incubated at 55 °C for about 20 hours, and with a nondestructive DNA extraction procedure to allow subsequent morphological observation. Genomic DNA extracts were stored in a freezer at -20 °C. The barcode region (630bp) of the COI gene was amplified using primer combination (Folmer et al. 1994), LCO1490 (5'-GGT CAA ATC ATA AAG ATA TTG G-3') and HCO2198 (5'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3') by the standard polymerase chain reaction (PCR) method. Total reaction volume was 25 μl, containing 12.5 μl of 2×Taq MasterMix, 8.5 μl of double distilled water (ddH 2 O), 2 μl of forward and reverse primer (1 μl, respectively), and 2 μl of DNA template solution. The following thermal cycling protocol was used: an initial denaturation step at 94 °C  Table 1. for 3 min, followed by 5 cycles of denaturation at 94 °C for 1 min, annealing at 45 °C for 1.5 min and extension at 72 °C for 1.5 min, followed by 35 cycles of denaturation at 94 °C for 1 min, annealing at 53.5 °C for 1 min and extension at 72 °C for 1 min, with a final extension of at 72 °C for 5 min, and ending with incubation at 12 °C.
The PCR products were examined using 1% agarose gel electrophoresis with ethidium bromide stain to check for successful amplification. The successful PCR products were sent to Beijing Tsingke Biotechnology Co., Ltd (China) for sequencing of both strands using the original PCR primers. All sequences collected in this study have been submitted to GenBank and accession numbers are shown in Table 1.

Molecular data analysis
The forward and reverse chromatograms were proofread and then assembled and edited using DNASTAR software (DNASTAR, Madison, Wisconsin, USA). Multiple sequence alignments were performed by CLUSTAL X 2.0.21 (Thompson et al. 1997;Jeanmougin et al. 1998). Primer sequences were manually deleted with BIOEDIT 7.0.9.0 (Hall 1999). To ensure that the correct target gene fragment was obtained, all sequences were checked in NCBI by Basic Local Alignment Search Tool (BLAST) (Altschul et al. 1990). To ensure nonexistence of stop codons and pseudogenes, the nucleotide sequences were translated to amino acids by MEGA 7 (Kumar et al. 2016). Sequence composition analyses were performed in MEGA 7. Pairwise genetic distances were calculated using the Kimura 2-parameter (K2P) model in MEGA 7 (Kimura 1980). Haplotypes were defined by DNASP 5.0 (Librado and Rozas 2009). The detailed statistics for haplotypes are shown in Table 1. The substitution saturation tests of 45 haplotype sequences segments were conducted in DAMBE 5.3.74 (Xia 2013) by comparing the index of substitution saturation (Iss) with critical values (Iss.c). To construct phylogenetic trees, neighbor joining (NJ), minimum evolution (ME), Bayesian inference (BI) and maximum likeihood (ML) analyses were performed. NJ and ME analyses (Saitou and Nei 1987) were performed in MEGA 7 under K2P substitution model. Branch support was measured using 1000 replicates in each analysis (Felsenstein 1985). Results were summarized as 50% majority consensus trees in MEGA 7. BI analysis was performed in MRBAYES 3.1.2 (Ronquist and Huelsenbeck 2003). The best-fit nucleotide evolution substitution model was selected by JMODELTEST 2.1.7 (Darriba et al. 2012). The Bayesian information criterion (BIC) was used to compare substitution models. The HKY+G model of nucleotide evolution was used. Two replicate runs with four independent Markov chain Monte Carlo (MCMC) chains (one cold chain and three hot chains) to conduct for 2 million generations, with trees sampled every 1000 generations with default parameter values. The average standard deviation of split frequency was lower than 0.01, indicating that the sampling of posterior distribution was adequate. The average standard deviation of split frequencies and Potential Scale Reduction Factor (PSRF) were used for examining convergence.  The stationarity was determined in TRACER 1.5 (Rambaut and Drummond 2009) by plotting the log-likelihood values versus generation number and the effective sample sizes >200 for all parameters. After stationarity had been reached, the first 25% trees were discarded as burn-in and a 50% majority-rule consensus tree with the posterior probability considered as node support values was constructed by summarizing the remaining trees. ML analysis was performed in RAXMLGUI 1.3.1, a graphical frontend for RAXML (Silvestro and Michalak 2012). All ML analyses with thorough bootstrap were run 10 times starting from random seeds under the GTRGAMMA model. The bootstrap support value (BS) was evaluated by analysis with 1000 replicates. All tree topologies were displayed in FIGTREE 1.4 (Rambaut 2009).

Morphological variation of D. vulgaris
Our specimens from China included representatives of both previously reported morphotypes of the aedeagus of D. vulgaris. They also exhibited a range of more subtle variation in the curvature of the aedeagal shaft in lateral view. Under the current morphology-based concept, this species can nevertheless be identified by the colour pattern and the presence of a shallow apical notch on the aedeagus in posterior view.

Sequence composition
The COI sequences are 630bp in length after alignment and trimming. Details of nucleotide composition are listed in Table 2. As is typical for insect mtDNA, the gene is AT-rich (Liu et al. 2012).

Substitution saturation test
The results of haplotype sequences for the substitution saturation test indicate the value of Iss is smaller than Iss.c; namely, little substitutional saturation was detected, which is strongly informative for constructing phylogenetic trees.

Analysis of the genetic distance and phylogenetic trees
The average genetic distances among morphological and geographical groups of D. vulgaris ranged from 0.9% to 6.3% and among species of Deltocephalus ranged from 16.4% to 21.9% without overlap (Table 3). This effectively reveals the existence of "DNA barcoding gap" and indicates the variation among morphological and geographical groups of D. vulgaris have not reached species level. Forty-five haplotypes (D. vulgaris -36, D. pulicaris -5, D. uncinatus -4) were defined to perform phylogenetic analyses. The phylogenetic analyses based on NJ, ME, BI and ML methods nearly yielded identical trees except for the slight change of the position of a few individuals of D. vulgaris and bootsrap values (Figs 4, 5). Deltocephalus vulgaris haplotypes grouped into several distinct clades. However, these groups included individuals of both morphotypes and formed a distinct monophyletic clade with strong support value (BS(NJ) = 100, BS(ME) = 100, PP = 1, BS(ML) = 97) with no obvious biogeographic structure. Furthermore, different morphotypes of D. vulgaris share the same haplotype (Table 1). Thus, the COI sequence data suggest that previous authors were correct in treating the two morphotypes of D. vulgaris as belonging to the same species.

Discussion
DNA barcoding as a standardised method to provide rapid and accurate species demarcation and has been widely applied in identifying and delimiting taxa since it was first reported by Hebert et al. (2003). Two standard criteria have generally been accepted in delimiting species using COI-based DNA barcodes. Based on the existence of a DNA barcoding gap, the feasibility of COI-based DNA barcoding depends on the fact that genetic distances among species are usually much higher than distances within species, without overlap. Different numbers of single species always form an independent clade in a phylogenetic tree (Wiens and Penkrot 2002;Hebert et al. 2003). Our analysis of COI sequences of Deltocephalus suggests a low level of genetic variation among morphotypes and geographical populations of D. vulgaris, and even different morphotypes of D. vulgaris share the same haplotype (e.g., YNA1 and YNB2; FJA1 and HNB1). The intergroup average genetic distances (0.9%-6.3%) of D. vulgaris among morphotypes and geographical populations is distinctly lower than that among species of Deltocephalus (16.4%-21.9%), without overlap. The phylogenetic tree (Figs 4, 5) recovered three independent lineages representing each of the three species with moderate to high support val-ues. The genetic distances among a few morphotypes and geographical populations of D. vulgaris exceeded the 3% standard threshold (e.g., ZJA and HNB; YNA and HNA). The more detailed genetic distances between groups/species of Deltocephalus are summarized in Table 3. However, all individuals of D. vulgaris grouped into a single clade with strong support comprising several subordinate clades but with no obvious correspondence to morphological or geographic groups. Furthermore, different morphotypes from the same and different geographical distributions of D. vulgaris share the same haplotype (Table 1). We consider that the intraspecific genetic distance of a 3% standard threshold can be an inconsistent standard in different groups and maternal inheritance of mitochondrial genes can be affected in the process of evolution by the same mode of inheritance as Wolbachia infection, which also may result in a higher divergence in host mtDNA (Frezal and Leblois 2008;Muñoz et al. 2011). The low level of variation among morphotypes and geographical populations of D. vulgaris supports the notion that they represent a single species. Differences in morphological characteristics, especially in male genitalia, have been the most reliable standard for discriminating among complex groups for many years. However, some cases of intraspecific variation in genital structures have been reported and these have led to uncertainty in the status of species and morphotypes. Mutanen et al. (2007) reported the male genital features that are most accepted and widely used standards to delimit species have been doubted in comparative study on male genital variation in Pammene luedersiana (Lepidoptera, Tortricidae). Yang et al. (2014) found 31 morphological variants in six species of Mogannia (Cicadidae, Cicadinae), but analysis of molecular data revealed low levels of intraspecific variation, although these morphological features have routinely been used to delimit species in this group. On the other hand, Wang et al. (2016) delimited seven different species of Bundera (Cicadellidae, Evacanthinae) based on molecular data, but only very minor morphological differences were found among six of these species. We are gradually becoming aware that similar morphological variation may have a different significance in different groups of leafhoppers and morphology-based species concepts may require confirmation using other kinds of data. DNA barcoding can efficiently complement morphology-based taxonomy and improve accuracy and rapidity in species identification.
Deltocephalus vulgaris, including 88 individuals in this study, and mainly representing two different forms of the aedeagus, were confirmed to be a single species grouped into a single clade with strongly support value in its phylogenetic trees (Figs 4, 5). Individuals collected both at the same place and time and different times and places have the same two forms of the aedeagus (e.g., FJA2 and FJB2; YNB10 and HNA2), which indicates forms are not related to temperature, humidity, precipitation, day length, altitude or latitude.
Our study shows a low intraspecific genetic distance between Guangdong and Hainan populations of D. vulgaris in southern China, suggesting that the Qiongzhou Strait (Fig. 1), a well-known biogeographic barrier has not significantly restricted gene flow for this species and they even share the same haplotype (Table 1). One logical assumption to explain this discovery is that Hainan and Guangdong arose earlier than the Qiongzhou Strait historically. Therefore, D. vulgaris freely exchanged genes when Guangdong and Hainan had been connected.
In the present study, lack of apparent correlation between morphology and COI haplotype is consistent with the hypothesis that the observed morphological variation is intraspecific. Nevertheless, we acknowledge the possibility that two different leafhopper species may share the same, or similar, COI haplotype. Thus, study of other genes may, in the future, reveal higher levels of divergence between the two forms and support recognition of some morphological variants as separate species.