All specimens of Lucanus were netted or light-trap collected for this study and store in ethanol, including 54 samples of the ingroup (21 L. wuyishanensis collected from Zhejiang, Fujian, Jiangxi province; 17 L. liuyei from Guangxi, Guizhou, Hunan province; five L. swinhoei from Taiwan; 11 L. continentalis from Fujian and Zhejiang province), and 10 samples of the outgroup (one each of L. parryi Boileau, 1899, L simithii Parry, 1862, L. fryi Boileau, 1911, L. klapperichi Bomans, 1989, and six L. fujianensis Schenk, 2008). Voucher specimens and their extracted genomic DNA are deposited in the research collection at the Museum of Anhui University, China. (Suppl. material 1).

The map with collection localities was generated using ArcGIS v. 10.3 (http://www.esri.com/sofware/arcgis) based on the geospatial data from the National Geomatics Center of China (Fig. 3). Photographs of the habitus was taken in .jpg format using a Canon 5D Mark IV with Canon 100 mm f/2.8 macro lens and a twin flash (Figs 1, 2).

Figure 3. 

Sample collection sites for this study.

The specimens were preserved in 99.7% ethanol at −20 °C. Total genomic DNA was extracted from a small portion of the muscle using DNeasy Blood and Tissue Extraction kit according to the manufacturer’s recommendations. The primers used to amplify 28S rDNA and Wingless were adapted from Abouheif and Wray (2002), Ward and Downie (2005), Monaghan et al. (2007), and Wild and Maddison (2008). For COI and 16S rDNA, primers were specifically designed in this study (Table 1).

PCR amplification reactions for the three loci (COI, 16S rDNA, 28S rDNA) were performed in a 25 μL volume containing 1 μL of each primer (forward and reverse) at 10 μM, 2 μL of template DNA solution, 12.5 μL of 2× EasyTaq SuperMix (+ dye), and 8.5 μL of sterile double-distilled water to make up the final volume of 25 μL. WG was amplified by nested PCR, first PCR containing 1 μL of each primer (forward and reverse) at 10 μM, 1 μL of template DNA solution, 7.5 μL of 2× EasyTaq SuperMix (+ dye), and 4.5 μL of sterile double-distilled water, and finally use the 1 μL first amplification product as a template, including 1 μL of each primer (forward and reverse) at 10 μM, 12.5 μL of 2× EasyTaq SuperMix (+ dye), and 9.5 μL of sterile double-distilled water. The polymerase chain reaction amplifications were performed under the following conditions: initial denaturation at 94 °C for 2 min, followed by 35–37 cycles of denaturation at 94 °C for 40 seconds, annealing at 52–60 °C for 50 seconds, and elongation at 70 °C for 1 min, and then a final extension step at 72 °C for 7 min, stored at 4 °C at room temperature. Amplifications were purified using Template DNA Amplify Kit (Ensure Biologicals).

Sequencing was performed using the ABI PRISM BigDye Terminator v. 3.1 Cycle Sequencing Kit (Life Technologies, USA), and cycle sequencing reactions were performed on ABI PRISM 3730xl automated sequencers (Life Technologies, USA) at Sangon Biotech Company, China. All sequences generated in this study were submitted to GenBank under accession numbers (Suppl. material 1).

Sequences of forward and reverse strands were assembled using GENEIOUS PRIME 2019.1.1 (https://www.geneious.com) and then aligned using MEGA 11. Genetic divergences among taxa were estimated using MEGA 11 (Kumar et al. 2018) via K2P-distance. The COI gene of Lucanus was assembled for genetic distance analyses. Finally, we concatenated alignments using PHYLOSUITE 1.2.1 (Zhang et al. 2020). The concatenated dataset was partitioned according to the Akaike Information Criterion (AIC) with PARTITIONFINDER 2.1.1 (Lanfear et al. 2017) for phylogenetic analyses.

Phylogenetic inferences were conducted using four gene markers based on maximum likelihood inference (ML), and Bayesian inference (BI). The BI tree was implemented in MRBAYES 3.2.6 (Ronquist et al. 2012). PARTITIONFINDER 2.1.1 was used to determine the best-fit models. Bayesian inference was conducted using MRBAYES 3.2.6 with two simultaneous runs of 5 × 10⁷ generations. Samples were drawn every 1,000 Markov Chain Monte Carlo (MCMC) step. The average standard deviation of split frequencies should be less than 0.01, with the initial 25% of trees discarded as burn-in. ML analyses were performed using IQ-TREE webserver (Trifinopoulos and Lam-Tung 2016). The “Auto” option was set under optimal evolutionary models, and the phylogenetic trees were constructed using an ultrafast bootstrap approximation approach with 10,000 replicates. Phylogenetic trees were visualized and edited in FIGTREE 1.4.3 (http://beast.bio.ed.ac.uk/figtree).

When defining species relationships by using molecular-set data, there are a variety of analytical approaches available. The Automatic Barcode Gap Discovery (ABGD) analysis was performed in this study for COI using a web interface (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html), which detects a gap in divergence distribution, which corresponds to differences between intraspecific and interspecific distances. When the gap exists, the process works well for species delimitation (Puillandre et al. 2012). An unrooted ML tree was generated in IQ-TREE webserver under the auto options. Poisson Tree Processes (PTP) was performed on a single, unrooted tree in the bPTP server (https://species.h-its.org/). A total of 5 × 10⁵ generations were run with the first 10% as burn-in (Zhang et al. 2013). General Mixed Yule Coalescent model (GMYC) analysis was conducted using BEAST 2.6.0 under a relaxed clock Exponential mode and ESS values assessed convergence. A burn-in with 25% was set to obtain an optimal consensus tree. Delimitation approach in the software R with the package ‘splits’ (available at http://r-forge.r-project.org/projects/splits) using the single-threshold method (Pons et al. 2006; Monaghan et al. 2009).

ABGD Automatic Barcode Gap Discovery;

BBP Bayesian posterior probability;

BI Bayesian inference;

GMYC General Mixed Yule Coalescent;

K2P distance, Kimura 2-parameter distance;

ML maximum likelihood;

MBL maximum likelihood bootstrap;

MOTU molecular operational taxonomic unit;

PCR Polymerase Chain Reaction;

PTP Poisson Tree Processes.

External morphological characteristics of two clades are compared in the table (Suppl. material 2). Lucanus swinhoei and L. continentalis differ slightly in body size, maintaining an overall range of 27–58 mm. Major males display the following features: mandibles strongly incurved at basal 1/3 and at apex; apical teeth bifurcated, upper branch teeth usually larger or equal in size with lower branch teeth; major mandibular tooth located at basal 1/3, triangular, proceeded by more than five inner small teeth. Medium-sized males have the following features: mandibles strongly incurved at basal 1/3, usually straight at apex; apical teeth bifurcated, upper branch teeth usually equal to or smaller than lower branch teeth; major mandibular tooth located at basal 1/3, triangular, proceeded by four or five inner small teeth. Minor males indicate the following features: mandibles weakly incurved at basal 1/3, straight at apex; apical teeth bifurcated, upper branch teeth usually equal to or smaller than lower branch teeth; major mandibular tooth located at basal 1/3, triangular, proceeded by fewer than four inner small teeth.

Lucanus wuyishanensis and L. liuyei differ slightly in body size, maintaining an overall range of 28–53 mm. Major males indicate the following features: mandible weekly incurved at basal 1/3, straight at the middle then strongly incurved at apical 1/4; the major inner mandibular tooth located 2/3 from the apical mandibular fork, sharp, triangular protruding forward and inflated on both sides, four separated small inner mandibular teeth attached below the major inner mandibular tooth, four or five unclear, minor inner mandibular teeth continuously located along the midlength of basal mandibles; four small inner mandibular teeth densely distributed between the major tooth and the apical fork. Medium-sized males indicate the following features: major inner mandibular tooth somewhat triangular, weekly inflated on both sides; more than four unclear, minor inner mandibular teeth continuously located along the midlength of basal mandibles, and more than four small inner mandibular teeth densely distributed between the major tooth and the apical fork. Minor males indicate the following features: major inner mandibular tooth weekly developed, single-point and not triangular; less than two separated small inner mandibular teeth attached below the major inner mandibular tooth; more than three unclear, minor inner mandibular teeth continuously located on 1/2 of basal mandible; more than two small inner mandibular teeth densely distributed between the major tooth and the apical fork.

A concentrated matrix with 2489 aligned positions for data was obtained comprising COI, 16S rDNA, 28S rDNA, and Wingless genes. The phylogenetic analyses using both BI and ML inferences recovered overall a consistent topology (Fig. 4). As outgroups, L. smithii, L. fryi, and L. parryi were separated from other species with high support, forming an independent clade. In addition, the clade L. klapperichi was sister to the clade (L. liuyei + L. wuyishanensis) (BPP = 1, MLB = 94). The clade L. fujianensis was sister to the clade (L. swinhoei + L. continentalis) (BPP = 1, MLB = 100). Nested structures of L. liuyei and L. wuyishanensis occur in all geographic populations.

Figure 4. 

Phylogenetic inferences based on four genes (COI, 16S rDNA, Wingless, and 28S rDNA) by maximum-likelihood inference (MLI) and Bayesian inference (BI) with posterior probability. Both posterior probabilities of MLI (above/left of branch) and bootstrapping values of BI (below/right of branch) are shown at nodes.

Genetic distances (K2P-distances) were calculated for all taxa using COI genes (Table 2). The results showed that the average genetic distance between L. wuyishanensis and L. liuyei populations in each collection area was a low mean range (0.0067–0.0110) (Suppl. material 3) and mean genetic distance of 0.0094 in all taxa (Table 2); the mean genetic distance between L. swinhoei and L. continentalis was 0.0072 (Table 2). The numbers were lower than the minimum mean genetic distances of 0.1592 among interspecific taxa and far less than the mean genetic distance of 0.2090 between interspecies of Lucanus (Lin 2017).

Species delimitation is shown in Fig. 5. Analysis of COI gene by all methods (ABGD, PTP, and GMYC) resulted in two molecular MOTUs, L. wuyishanensis + L. liuyei and L. swinhoei + L. continentalis (Fig. 5). For the concatenated dataset, all three methods suggested that L. wuyishanensis + L. liuyei were one MOTU, whereas GMYC divided L. swinhoei and L. continentalis into two MOTUs.

Figure 5. 

Delimitation of the studied Lucanus species based on COI. Columns are taxonomic identification based on three molecular delimitation methods: ABGD, PTP, and GMYC. The phylogenetic tree is based on the GMYC analysis.