Integrative species delimitation based on COI, ITS, and morphological evidence illustrates a unique evolutionary history of the genus Paracercion (Odonata: Coenagrionidae)

Taxon sampling

The geographic distribution of Paracercion is shown in Fig. 1. A total of 348 individuals of the genus Paracercion were collected from China, Japan, and Vietnam. For ITS, 16 sequences of Paracercion and seven sequences of out-groups downloaded from NCBI (National Centre for Biotechnology Information, https://www.ncbi.nlm.nih.gov/) were also included in the analyses. In total, nine nominal species of Paracercion were analyzed with the molecular data. In addition, the sequences of five genera: Coenagrion, Enallagma, Pseudagrion, Ischnura and Ceriagrion were chosen as out-groups. The geographic coordinates of different locations were recorded using a global positioning system (GPS) tool. For ease of DNA amplification, most of the specimens were first deposited in anhydrous alcohol in the field and then stored in a freezer at −20 °C at the College of Life Sciences, Chongqing Normal University. The detailed collecting information of the 348 specimens is shown in Table S1.

DNA extraction, amplification and sequencing

One hind leg of each specimen was used to extract the genomic DNA by the Universal Genomic DNA Kit (Beijing CoWin Biosciences, Changping, China) following manufacturer protocol. The mitochondrial gene cytochrome c oxidase subunit 1 (COI) and a nuclear segment including ITS1-5.8S-ITS2 (ITS) were first amplified following Yu et al. (2015). However, the amplification success rate of the ITS region was relatively low even using the clone method. All PCR products were visualized using 1% agarose gel electrophoresis and purified using a gel extraction kit (Sangon Biotech, Shanghai, China). The PCR products were sequenced by Beijing Genomics Institute or GENEWIZ based on Sanger’s chain termination method.

Sequence analysis

The sequences files were viewed and edited in BioEdit 7.2.0 (Hall, 1999). Alignments of the COI sequences were translated into amino acids using MEGA X (Kumar et al., 2018) to detect the presence of pseudogenes. All of the sequences were saved as FASTA format and aligned with MAFFT version 7 (Katoh & Standley, 2013).

OTU delimitations

The ABGD analysis was performed using the web version (http://wwwabi.snv.jussieu.fr/public/abgd/). The detailed arguments were as follows: Pmin = 0.001, Pmax = 0.1, setp = 10, relative gap width = 1.5, Nb bins (for distance distribution) = 20 and the genetic distance model was set to Kimura (K80).

Distance-based clustering analysis was performed with the “tclust” function of the R package spider (Brown et al., 2012). Fixed empirical thresholds: 1%, 2%, 2.2% and 3% were selected following Zhang et al. (2017). The genetic distances required for clustering analysis were calculated with the “dist.dna” function of the ape package (Paradis, Claude & Strimmer, 2004) in R 3.6.3 (R Core Team, 2020) using the K2P model (Kimura, 1980).

For GMYC analysis, the ultrametric trees were generated in Beast 2.5.2 (Bouckaert et al., 2019). Models of DNA sequence evolution were estimated by MrModeltest v2.3 (Nylander, 2004) based on Akaike Information Criterion (AIC). The nucleotide substitution model selected for COI and ITS was GTR+I+G. We used the relaxed lognormal clock model and the coalescent constant population tree prior. For each molecular marker, the length of Markov chain Monte Carlo (MCMC) chains was set to 20,000,000 (ITS) and 50,000,000 (COI) with trees sampled every 1000 generations. The convergence of each run was checked with Tracer v. 1.7 (Rambaut et al., 2018). The first 10% of trees of each run were discarded, and the maximum clade credibility tree (the posterior probability limit was set to 0.5) was reconstructed using TreeAnnotator (embed in BEAST package) with the mean node heights option. Afterward, the maximum clade credibility tree was used as input for the GMYC analysis with the splits package (Ezard, Fujisawa & Barraclough, 2013) in R. The single-threshold method was used to generate the putative species.

For bPTP analysis, the Maximum-likelihood trees were generated with IQ-TREE Version 1.6.10 (Nguyen et al., 2015). We obtained branch supports with the ultrafast bootstrap (Hoang et al., 2018) implemented in the IQ-TREE software. ModelFinder (Kalyaanamoorthy et al., 2017) was used to compute the log-likelihoods of an initial parsimony tree for different models and to select the best-fitting models of sequence evolution. The out-groups selected for rooting phylogenetic trees were removed in the bPTP analysis. The number of MCMC generations was set to 200,000; the thinning value was set to 100 and the burn in value was set to 0.1. Other parameters were set to default values. The maximum likelihood solution result of bPTP is discussed in this study.

The ultrametric tree and the combined results of different species delimitation methods were visualized by the R package ggtree (Yu et al., 2017).

Network analysis

The sequences files were transformed into haplotype sequences using the R package haplotypes (Aktas, 2020) with gaps and missing data included. Thereafter the haplotype network of each molecular marker was reconstructed with PopART 1.7 (Leigh & Bryant, 2015) using the TCS method (Clement, Posada & Crandall, 2000).

Morphologic checking

All of the living photos were taken in the field with a digital camera (Nikon D3200) (Bangkok, Thailand). Character photos were taken in the laboratory using the Zeiss V20 (Jena, Germany) microphotography system. Specimens were examined and dissected under a Zeiss V8 stereomicroscope. More than 600 specimens (including dried specimens) were identified by Xin Yu and Xin Ning in terms of the corresponding literatures (Dumont, 2004; Yu, 2008; Yu & Bu, 2011; Ning et al., 2016). All the specimens were stored at the College of Life Sciences, Chongqing Normal University. The identification of the Paracercion species is not easy due to the prevalent convergent in appearance and intraspecific variations. The most important diagnostic morphologic characters are the shapes of male caudal appendages and genital ligulae and these were carefully examined. Generally, the variation of male genital ligulae among species is limited. This cannot provide more information on species delimitation other than to show some degree of relationship. The structure of the caudal appendages is the most credible diagnostic character. Based on this character, all the individuals can be divided into nine morphological defined species.

Molecular data

We obtained molecular data from all the nine species of Paracercion. The COI dataset consisted of 348 sequences from Paracercion and 25 sequences from out-groups with a length of 1,066 bp. There were 312 variable sites and 282 parsimony-informative sites in 348 COI sequences and most of the variation occurred in the third codon position. The ITS dataset consists of 66 sequences from Paracercion and seven sequences from out-groups. We obtained both COI and ITS sequences from all the nine species of Paracercion.

Genetic distance (K2P distance)

The COI genetic distances between P. hieroglyphicum and P. melanotum ranged from 0–1.65%. In comparison, the ITS genetic distances between P. hieroglyphicum and P. melanotum ranged from 9.58% to 10.89%. Similarly, the COI genetic distances between P. barbatum and P. v-nigrum ranged from 0.19% to 1.16%; the ITS genetic distances between P. barbatum and P. v-nigrum ranged from 3.03% to 8.44%. The COI genetic distances of P. barbatum ranged from 0 to 0.19% (average value: 0.08%, Table S2). The ITS genetic distances of P. barbatum ranged from 0 to 7.94% (average value: 3.01%, Table S3). Regardless of the specimen, ChH03, the ITS genetic distances of P. barbatum ranged from 0 to 3.19%. The COI genetic distances between P. barbatum and P. melanotum ranged from 9.46% to 10.58%. The ITS genetic distances between ChH03 and 11 specimens of P. melanotum ranged from 0.16% to 2.17%; the ITS genetic distances between the other specimens of P. barbatum and the 11 specimens of P. melanotum ranged from 5.39% to 8.49%. Although the COI genetic distances between P. sieboldii and P. v-nigrum ranged from 3.13% to 4.05%, the ITS genetic distances between P. sieboldii and P. v-nigrum were relatively small (0.17–4.13%). The overlap between intraspecific and interspecific genetic distances indicated that the barcode gap does not exist in the COI and ITS markers (Fig. S1, Fig S2).

ABGD analysis

For COI (P = 0.0129), 348 specimens were divided into six groups. The specimens of P. barbatum, P. sieboldii and P. v-nigrum clustered into one group (Group 1). The specimens of P. hieroglyphicum and P. melanotum clustered into the same group (Group 4). The specimens of the other four morphologically defined species: P. ambiguum, P. calamorum, P. dorothea and P. plagiosum were each divided into one of four different groups (Table S4).

For ITS (P = 0.0129), 66 specimens were divided into seven groups. The specimens of P. barbatum (excluding ChH03), P. sieboldii and P. v-nigrum clustered into one group. One specimen (ChH03) of P. barbatum and the 11 specimens of P. melanotum were placed in the same group. The specimens of the other five morphologically defined species: P. ambiguum, P. calamorum, P. dorothea, P. hieroglyphicum and P. plagiosum were divided into five different groups (Table S4).

Distance-based clustering analysis

If the threshold was set to 1%, 2% and 2.2%, the number of MOTUs/clusters defined by clustering analysis on COI was seven. In comparison, the number of MOTUs for ITS was 17, 11 and 10 separately (Fig. S3). For COI, if the threshold was set to 3%, 348 specimens could be divided into seven MOTUs. The specimens of P. barbatum and P. v-nigrum were grouped into the same MOTU; the specimens of P. hieroglyphicum and P. melanotum were also grouped into the same MOTU. The specimens of the other five morphologically defined species—P. ambiguum, P. calamorum, P. dorothea, P. plagiosum and P. sieboldii—were divided into five different MOTUs (Table S5).

For ITS, if the threshold was set to 3%, 66 specimens could be divided into seven MOTUs. The specimens of P. barbatum (excluding ChH03), P. sieboldii and P. v-nigrum were grouped into the same MOTU. One specimen (ChH03) of P. barbatum and the 11 specimens of P. melanotum were divided into the same MOTU. The specimens of the other five morphologically defined species—P. ambiguum, P. calamorum, P. dorothea, P. hieroglyphicum and P. plagiosum—were divided into five different MOTUs (Table S5).

GMYC analysis

For COI, 145 haplotypes could be divided into seven putative species (Table 1). The 70 haplotypes of P. barbatum and P. v-nigrum were grouped into the same putative species; the 55 haplotypes of P. hieroglyphicum and P. melanotum were also grouped into the same putative species. The haplotypes of the other five morphologically defined species: P. ambiguum, P. calamorum, P. dorothea, P. plagiosum and P. sieboldii were divided into five different putative species (Table S6).

For ITS, 50 haplotypes could be divided into 12 putative species (Table 1). One haplotype (Hap_4, ChH03) of P. barbatum and the nine haplotypes of P. melanotum were included within the same putative species. Six haplotypes of P. barbatum were separated into three putative species. Five haplotypes of P. calamorum were separated into two putative species. Sixteen haplotypes of P. sieboldii and P. v-nigrum were separated into two putative species. The haplotypes of the other four morphologically defined species: P. ambiguum, P. dorothea, P. hieroglyphicum and P. plagiosum were divided into four different putative species (Table S6).

bPTP analysis

A total of 145 COI haplotypes could be divided into seven tentative species using bPTP analysis. All the haplotypes of P. barbatum and P. v-nigrum were grouped into the same tentative species; all the haplotypes of P. hieroglyphicum and P. melanotum were also grouped into the same tentative species. The haplotypes of the other five morphologically defined species—P. ambiguum, P. calamorum, P. dorothea, P. plagiosum and P. sieboldii—were divided into five different tentative species (Table S7). The result of bPTP analysis on COI was consistent with the GMYC analysis.

A total of 50 ITS haplotypes could be divided into 10 tentative species by bPTP analysis. One haplotype (Hap_4, ChH03) of P. barbatum and the nine haplotypes of P. melanotum clustered into the same tentative species. Six haplotypes of P. barbatum were separated into three tentative species. Sixteen haplotypes of P. sieboldii and P. v-nigrum were merged into the same tentative species. The haplotypes of the other five morphologically defined species—P. ambiguum, P calamorum, P. dorothea, P. hieroglyphicum and P. plagiosum—were divided into five different tentative species (Table S7).

The combined results of species delimitation of different methods based on COI and ITS are shown in Fig. 2 and Fig. 3, respectively.

Network analysis

The COI haplotype network showed that the haplotypes of P. barbatum and P. v-nigrum were merged into one large group; however, they did not share any haplotypes. In comparison, the haplotypes of P. sieboldii were placed into a single group (Fig. 4A). All the haplotypes of P. hieroglyphicum and P. melanotum were also grouped into one large group. Two individuals of P. hieroglyphicum and two individuals of P. melanotum shared the same haplotype (Fig. 4B).

The ITS haplotype network showed that the haplotypes of P. sieboldii and P. v-nigrum were merged into one nested group; the haplotypes of P. barbatum were separated into a single group (Fig. 5A). In addition, one haplotype of P. barbatum and the nine haplotypes of P. melanotum were merged into one group (Fig. 5B).

Introgression in Paracercion

For COI, four methods including ABGD, Clustering, GMYC and bPTP all indicated that P. hieroglyphicum and P. melanotum belong to one putative species. For ITS, the results of all four methods support the conclusion that these two species are distinct. In addition, there are significant differences between P. hieroglyphicum and P. melanotum based on morphological characters (Fig. 6), especially the male caudal appendages (Fig. 7). These differences support the conclusion that these are two separate species. P. hieroglyphicum and P. melanotum are sympatric species in many areas (Fig. 1) and share similar shaped genital ligulae (Fig. 8). They have identical COI haplotypes (Fig. 4B). The morphological and molecular evidence argue for the introgression of mtDNA in P. hieroglyphicum and P. melanotum. Our field observations indicate that these two species have evolved reproductive isolation mechanisms. Males sometimes attempt to mate with nonconspecific females, but the latter are not receptive.

In general, introgression in Paracercion was found in COI rather than ITS, which indicates that mtDNA may be more likely to show introgression than nDNA (Ballard & Whitlock, 2004).

Hybridization between P. barbatum and P. melanotum

For ITS, four methods of species delimitation indicated that at least one specimen (ChH03) of P. barbatum together with all specimens of P. melanotum should be the same tentative species. This indicates that P. barbatum and P. melanotum likely hybridize at Changhu, Yunnan. Several possible hybrid individuals from Changhu were found based on morphological characters (Fig. 9). P. barbatum is restricted to Northeast Yunnan while P. melanotum is widely distributed (Fig. 1), covering the range of P. barbatum. It should be extended further south since the two previous southeast species P. malayanum and P. pendulum were confirmed to be the same species. Among them, the caudal appendage of a male specimen of P. malayanum in the RMNH (Natural History Museum, London, UK) collection (kindly provided by Dr. Rory A. Dow) (Fig. 7O) is identical to male P. melanotum. Shapes of the male caudal appendages and genital ligula of P. barbatum are similar to that in P. melanotum (Fig. 9) and that may be why Needham and Gyger confused them. This also suggests the possibility of hybridization. It is worth noting that recent hybridization events can be detected by ITS rather than COI.

Species group of ‘v-nigrum’

Results from both COI and ITS supported a very closely related species group including P. barbatum, P. sieboldii and P. v-nigrum. For COI, four methods—ABGD, Clustering, GMYC and bPTP—all indicated that P. barbatum and P. v-nigrum are a single putative species. The low mean interspecific divergence between P. barbatum and P. v-nigrum (COI: 0.49%) indicates a recent speciation event in Paracercion. For ITS, the results of GMYC and bPTP support that this pair of species were separated into two or more putative species. Furthermore, P. barbatum seems to be closely related to P. melanotum due to the hybridization discussed above. All the specimens of P. sieboldii and P. v-nigrum can be separated into two putative species in the four different methods based on COI. However, they were merged into the same putative species in ABGD, Clustering and bPTP analyses based on ITS. This discordance of species delimitation may be the result of incomplete lineage sorting in the nuclear DNA.

P. barbatum, P. sieboldii and P. v-nigrum are allopatric. However, their distributions are obviously contiguous, namely, P. v-nigrum has the widest distribution range while P. barbatum is restricted to the southwest corner (northwest Yunnan) and P. sieboldii is isolated at the northeast portion (Japan) of this area (Fig. 1). P. barbatum, P. sieboldii and P. v-nigrum share similar morphological characters especially the shape of genital ligulae (Fig. 6). Therefore, these three species may have originated from one widely distributed ancestral species. Dumont (2004) and Ning et al. (2016) also suggested that P. sieboldii shared a common ancestor with the continental P. v-nigrum.

Species delimitation methods

All methods of species delimitation will sometimes fail to delimit species boundaries properly or will give conflicting results. The use of qualitative judgment will then be necessary (Sites & Marshall, 2004). For species delimitation studies, multiple methods based on different principles and multiple loci underlying different evolutionary history should be included. Theoretically, the most realistic strategy is to accept delimitations that are congruent across methods and may represent actual evolutionary lineages. However, morphological characters are also important in species delimitation studies. Inferences drawn from species delimitation studies should be conservative (Carstens et al., 2013). For ITS, the number of putative species generated from bPTP and GMYC exceeded the total number of species determined from morphological characteristics. This was mainly due to the high rates of false positives of PTP and GMYC methods (Luo et al., 2018).