Molecular phylogeny of the genus Dicronocephalus (Coleoptera, Scarabaeidae, Cetoniinae) based on mtCOI and 16S rRNA genes

Abstract The seven species belonging to the genus Dicronocephalus are a very interesting group with a unique appearance and distinct sexual dimorphism. Only one species among them, Dicronocephalus adamsi, has been known in the Korean fauna. This species is recognized as having a wide distribution from Tibet to Korean Peninsula and is currently represented by two subspecies that have separated geographical ranges. The phylogenetic relationships of Dicronocephalus adamsi were still unclear. The phylogeny of Dicronocephalus is reconstructed with a phylogenetic study of five species including four subspecies based on a molecular approach using mitochondrial COI and 16S rRNA genes. Our results are compared with the results obtained by previous authors based on morphological characters. They show that the tested taxa are divided into two major clades. Clade A consists of two species (Dicronocephalus adamsi + Dicranocephalus yui) and Clade B includes the others (Dicronocephalus dabryi + Dicranocephalus uenoi + Dicranocephalus wallichii). This result generally supports Kurosawa’s proposal except that Dicronocephalus dabryi and Dicranocephalus uenoi are newly recognized as members of a monophyletic group. We propose that Dicronocephalus adamsi drumonti is a junior subjective synonym of Dicronocephalus adamsi adamsi. These results show that three members of the Dicranocephalus wallichii group should be treated as species rather than subspecies. However, further research including analyses of different genetic markers is needed to reconfirm our results.

Among the seven species of Dicronocephalus, only D. adamsi is found in the Korean fauna. This species was described from Korea, but it has been known to have a wide range across Korea, China, Tibet, and Vietnam. The range of this species is divided by a wide geographical gap between Liaoning and Shanxi provinces of China (Young 2012). Legrand (2005) divided D. adamsi into two subspecies based on this distribution pattern and morphological differences. He described populations occurring in west China as D. adamsi drumonti. This classification was accepted by Krajcik (2014), but not by Young (2012).
The subspecies of D. wallichii (D. w. wallichii, D. w. bourgoini, and D. w. bowringi) were originally described as valid species (Hope 1831, Pascoe 1863, Pouillaude 1914. While some authors have treated these taxa as subspecies (Paulian 1960, Mikšić 1971, 1977, Krajčík 1998, Sakai and Nagai 1998, Šípek et al. 2008, Young 2012, Krajcik 2014, some others have treated them as species (Kurosawa 1968, Devecis 2008. The controversy over whether they should be dealt with at the species or sub-species level has continued without in-depth analysis. During a review of the genus Dicronocephalus, several issues were encountered, such as validation of species or subspecies rank of taxa composing D. adamsi and D. wallichi (sensu lato) and the lack of phylogenetic analysis of the genus. To resolve these questions, phylogenetic analysis was performed for the genus using cytochrome c oxidase subunit I (COI) and 16S ribosomal RNA (16S rRNA) mitochondrial gene sequences as well as examination of their morphological diagnostic characters.

Specimen sampling and examination
Fifty specimens of Dicronocephalus belonging to five species and seven subspecies from four countries were obtained ( Fig. 1, Table 1), but we were unable to obtain specimens of the remaining two species, D. bieti and D. shimomurai. For examining male genitalia, these were extracted from the abdomens and cleaned by heating with 10% KOH solution in a WiseTherm®HB-48P heating block at 60 °C for 1~2 hours. Male genitalia were preserved in microvials with glycerine after examination. Photographs of external morphology and genitalia were taken with a Canon EOS 10D camera and stacked with a combineZM program (Hadley 2006). Based on previous studies (Pascoe 1863, Pouillaude 1914, Kurosawa 1968, 1986, Young 2012, diagnostic characters were obtained to provide precise criteria for species identification. In this study, the most recent taxonomic scheme by Krajcik (2014) was followed, especially for subspecies treatment of D. wallichii. All examined specimens are stored in the Department of Agricultural Biology, National Academy of Agricultural Biology (NAAS), Jeonju, Korea.

DNA extraction, amplification and sequencing
Genomic DNA (gDNA) was extracted from middle legs removed from dried specimens of all species and accomplished using a QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. Polymerase Chain Reaction (PCR) was performed in order to amplify the cytochrome c oxidase subunit I gene (COI) and 16S ribosomal RNA gene (16S rRNA) using Accupower PCR PreMix (Bioneer, Daejeon, Korea). The universal primer set LCO1490/HCO2198 (Folmer et al. 1994) for amplifying the DNA barcoding region (658bp) of COI sequences was not successful for all samples; this may be caused by the degraded quality of gDNA (Goldstein andDesalle 2003, Hajibabaei et al. 2006;Wandeler et al. 2007). We applied the PCR methodology for retrieving COI sequences from old specimens given in Han et al. (2014) and designed new primer pairs: LCO-Ceto232F (5'-GCHTTYC-CYCGAATAAATAAYATA-3') corresponding to HCO2198 and HCO-Ceto367R (5'-ACDGTYCADCCNGTTCCTGCNCC-3') corresponding to LCO1490. 16S rRNA was targeted in a 600 bp region with two primers, 16SB/16SA, that successfully amplified in Lucanidae and Elateridae (Hosoya et al. 2001, Hosoya and Araya 2005, Han et al. 2009. PCR amplification conditions were as follows: for COI, initial denaturation at 94 °C for 5 min, then 45 cycles at 94 °C for 30 s, 46 °C for 25 s, and 72 °C for 45 s followed by a final extension at 72 °C for 3 min, and for 16S rRNA, initial denaturation at 94 °C for 5 min, then 40 cycles at 94 °C for 1 min, 50 °C for 1 min, and 72 °C for 45 s followed by a final extension at 72 °C for 5min. The amplicons were purified using a QIA quick PCR Purification Kit (Qiagen, Hilden, Germany) after the product yield was monitored by 0.7% agarose gel electrophoresis.   (Table 1).

Phylogenetic analysis
For the phylogenetic analyses, three data sets were used, a 658 bp fragment of COI, 520 bp fragment of 16S rRNA sequences, and the concatenated COI and 16S rRNA sequences. The data sets were aligned using ClustalW in MEGA 5.2 (Tamura et al. 2011), and genetic distances were calculated using Kimura's two-parameter test (Kimura 1980). The phylogenetic analyses were constructed using maximum likelihood (ML), Bayesian inference methods (BI), and maximum parsimony (MP). ML analysis was performed with GARLI 2.0 (Zwickl 2011), and the analysis was initiated at a random start tree using GTR+I+G model parameters selected by Mr-ModelTest (Nylander 2004), with a 10,000 generation search algorithm and 1,000 bootstrap replications. The frequencies with which to log the best score ("logevery") and to save the best tree to file ("saveevery") were set to 10,000 and 10,000 respectively, and the number of generations without topology improvement required for termination ("genthreshfortopoterm") was set to 5,000. At the end of the analysis, there was no improvement in the tree topology by a log likelihood of 0.01 or better. The bootstrap values were calculated using the SumTrees program of the DendroPy package (Sukumaran and Holder 2010).
BI analysis was performed with MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003). Metropolis-coupled Markov chain Monte Carlo (MCMC) analyses were run with one cold and three heated chains (temperature set to 0.2) for 5,000,000 generations and tree sampling every 100 generations. The posterior probabilities were then obtained and a majority-rule consensus tree was generated from the remaining trees after discarding the first 25% of samples.
MP analysis was performed with TNT 1.1 (Goloboff et al. 2008). The analyses, followed by tree bisection reconnection (TBR) branch swapping, used default options that performed 100 random additional sequences and saved up to ten trees per replication. To obtain the strict consensus tree, symmetric resampling (Goloboff et al. 2003) with a 33% change probability and jack-knifing with a 36% removal probability were implemented using a traditional search with 1,000 replications. Each set of results was summarized in terms of absolute frequency, and the group support values were analyzed. For bootstrap value (BP) in ML and MP, and posterior probability value (PP) in BI, supporting values of <70% as "weak", 70-79% as "moderate", 80-89% as "strong", and ≥ 90% as "very strong" support were used.

Nucleotide information for COI and 16S rRNA
The data set of COI, with no evidence of indel (insertion/deletion) events, had 144 (21.9%) variable sites (Vs). Of these, 140 (21.3%) were parsimoniously informative sites (PIs). The data set of 16S rRNA, with indel events at three sites, consisted of 43 (8.3%) Vs, of which 41 (7.9%) were PIs. There was about 2.6 times more variability and the level of PIs was about 2.7 times greater in COI than in that in 16S rRNA.

Phylogenetic analyses of COI
Phylogenetic inferences based on three analyses (ML, BI, and MP) reconstructed the same topologies for COI ( Fig. 2; for BI, ML and MP tree data not shown, see Suppl. material 1 for sequences), and there was separation into two major clades (A and B) with very strong supporting values (100%), except for ML. Eight ingroup taxa representatives including subspecies were clearly clustered into seven monophyletic groups corresponding to nominal species; the two subspecies of D. adamsi formed one cluster. Their terminal nodes were well supported, but the values of ML and BI were very low in D. yui yui (<50% in ML and 53% in BI) and D. wallichii bowringi (<50% in ML and 56% in BI).
Clade A is composed of D. adamsi adamsi, D. a. drumonti, and D. yui yui with strong bootstrap support (>72%). The two subspecies of D. adamsi did not separate into two distinct subgroups. The genetic divergences between the two subspecies were relatively low (0-1.7%); moreover, D. a. drumonti shared haplotypes with D. a. adamsi from Korea and China. D. yui yui was sister to D. adamsi with distinct inter-specific divergences (5.6%-7.3%).
Clade B is composed of D. dabryi, D. uenoi katoi, and three subspecies of D. wallichii with strong bootstrap supports by ML and BI, but relatively low support (56%-62%) by MP. Among the members of Clade B, D. dabryi and D. uenoi katoi formed a monophyletic group with very strong supporting values in all analyses and with distinct inter-specific divergences (5.6%-8.9%). The intra-specific divergences of these two species (0-1.5% in D. dabryi, 0.2%-2.3% in D. u. katoi) were explicitly lower than their inter-specific values. The three subspecies of D. wallichii were clustered as a monophyletic group and clearly subdivided. D. w. bowringi diverged early from an ancestor, and then D. w. wallichii and D. w. bourgoini underwent subsequent separation with strong bootstrap supports by ML (83%) and BI (99%); however, despite low divergences within each subspecies ranging from 0.3%-0.8%, the genetic divergences between these subspecies were unexpectedly variable ranging from 2.7%-8.1%. Genetic divergences were larger between D. w. bowringi and both D. w. wallichii (4.3%-5.0%) and D. w. bourgoini (4.8%-8.1%), than those between D. w. wallichii and D. w. bourgoini (2.7%-5.7%).

Phylogenetic analyses of 16S rRNA
ML, BI, and MP analyses of 16S rRNA resulted in considerably similar topologies to those of COI ( Fig. 3 for BI, ML and MP tree data now shown, see Suppl. material 2 for sequences), but a polytomy was found in D. yui yui and paraphyly in D. w. bowringi with respect to D. w. wallichii.
The intra-specific pairwise distances of 16S rRNA were relatively low, ranging from 0-0.4%. The inter-specific divergences ranged from 0.8%-6.3%. The distances between the ingroup and outgroup taxa ranged from 9.7%-11.8% (Table 3) Figure 2. Phylogenetic relationships among Dicronocephalus species reconstructed with Bayesian inference using COI sequences. Numbers above branches indicate ML bootstrap values and Bayesian posterior probabilities. Numbers below branches are bootstrap, symmetric resampling, and jacknife support from parsimony searches, respectively. Scale bar represents 10% nucleotide mutation rate.    Figure 3. Phylogenetic relationships among Dicronocephalus species reconstructed with Bayesian inference using 16S rRNA sequences. Numbers above branches indicate ML bootstrap values and Bayesian posterior probabilities. Numbers below branches are bootstrap, symmetric resampling, and jacknife support from parsimony searches, respectively. Scale bar represents 10% nucleotide mutation rate. lowest inter-specific divergence range (0.8%-1.2%) was revealed between D. adamsi and D. yui yui, and this is rather similar to the divergence ranges of the D. wallichii subspecies (0.8%-1.6%).
Dicronocephalus adamsi was clustered as a sister to D. yui yui in Clade A with strong bootstrap support (>90%), while the remaining taxa were clustered into Clade B with relatively low supporting values (>76%) in BI and MP. The monophyly of D. adamsi, D. uenoi katoi, D. w. wallichii, and D. w. bourgoini was well supported by bootstrap analyses (>84%). In contrast, in all analyses a polytomy was found in D. yui yui and ML and BI showed paraphyly of D. w. bowringi. We showed that these phenomena were caused by few parsimony-informative nucleotide variations in conserved regions. A comparison of each of those sequences, showed that D. y. yui has different substitutions at 326 nucleotide position. Two samples (7290 and 7291) have "C", while one sample (7292) has "T". On the other hand, D. w. bowringi has a substitution occurred in 196 nucleotide position. The 7693 sample has "G", while the other samples (7692, 7694, and 7695) and two samples (7274 and 7275) of D. wallichii have "A" at this site (Suppl. material 2).  Numbers are indicated as mean (minimum-maximum) of the pairwise distance. *denotes outgroup taxon

Phylogenetic analyses of COI and 16S rRNA
In the combined data set of COI and 16S rRNA, phylogenetic reconstructions produced topologies congruent with the COI analyses. The nodal supporting values were improved compared with the analyses based on each gene (Fig. 4, see Suppl. material 3 for sequences). Monophyly of the seven taxa including subspecies was strongly supported by bootstrap values >90%, except for low support of 53% and 55% in ML and BI, respectively, for the terminal node of D. w. bowringi. D. w. wallichii was grouped as a sister to D. w. bourgoini based on the results of the COI analyses with a high value in BI (94%) and moderate value in ML (74%), but not in MP (Fig. 4).

Re-examination of morphological diagnostic characters
The 19 diagnostic characters used to classify species or subspecies were re-examined in order to determine whether they are suitable for identification (  Figure 4. Phylogenetic relationships among Dicronocephalus species reconstructed with Bayesian inference using COI and 16S rRNA sequences. Numbers above branches indicate ML bootstrap values and Bayesian posterior probabilities. Numbers below branches are bootstrap, symmetric resampling, and jacknife support from parsimony searches, respectively. Scale bar represents 10% nucleotide mutation rate.  (Fig. 1) 0) grayish brown Kurosawa (1968) 1) dark brown 2) yellowish brown 3) dark yellowish brown 4) green-yellowish brown with pale purple on elytra 2. Color in female 0) dark blackish body without marking Kurosawa (1986) 1) not dark blackish body 3. Pronotal and elytral colors (Fig. 1 (Fig. 6) 0) with triangular umbone Pascoe (1866) 1) without triangular umbone 15. Apicosutural angle (Fig. 7) 0) rounded Pouillaude (1914) 1) projected Metasternum 16. Metasternal process 0) obtuse, rather rounded Kurosawa (1968) Young (2012) 1) rectangular or acute, moderately produced 2) triangularly and sharply produced mentioned in previous studies, 13 are clearly suitable for species or subspecies identification; however, we recognized six characters that are ambiguous and not applicable (Table  5). For example, Pouillaude (1914) mentioned three diagnostic characters as follows: 1) D. dabryi has a different color of the pronotum and the elytra compared with D. wallichii subspecies (Fig. 1); 2) D. w. wallichii can be separated from the others (D. adamsi, D. w. bowringi, D. w. bourgoini, D. dabryi, and D. beiti) by having no angular projection at the base of the anterior edge of the clypeus (Fig. 5); and 3) D. w. bourgoini can be distinguished from the others by the projected apicosutural angle of the elytra (Fig. 6). However, none of these characters has proven to be suitable for species identification. We observed that the color of the pronotum and the elytra of D. dabryi was the same with grayish powder in freshly collected specimens, but it has faded gradually in old specimens (Fig. 1D). Also the anterior edge of the clypeus of D. w. wallichii (Fig. 5G) was sinuate in the middle, similar to that of D. w. bourgoini (Fig. 5H), and did not match the description by Pouillaude. We therefore consider that these characters might have been mistakenly described and illustrated by Pouillaude (1914). In addition, the projection of the apicosutural angle of the elytra of D. w. bourgoini was not distinct and could not separate this taxon from the other   Parentheses denote the characteristic represeted by our examination.
Question marks indicate the ambiguous character state to be difficult determination in our examination.
'C' is clear and 'U' is unclear characters resulted in this study.
species and subspecies (Fig. 6H). We consider that using another character such as "the posterior margin of the elytra is round or truncated" may more diagnostic than the former character as shown in Fig. 6. Pascoe (1863) used the triangular umbone on the shoulder of the elytra (Fig. 7) to distinguish D. a. adamsi from D.w. bowringi. But, we consider that the presence of a triangular umbone is as an unsuitable character. We found this state also in some specimens of D. adamsi, although the size of the triangular umbone was small and variable in each specimen. Kurosawa (1986) used the widest portion of the pronotum as a distinguishing character state, but this was variable in all specimens of D. w. bourgoini and not distinct enough to be used in species and subspecies identification. Legrand (2005) used six diagnostic characters to distinguish between the two subspecies, D. a. adamsi and D. a. drumonti. Among them, we found four characters, namely body size, general body shape, longitudinal bands on the pronotum, and the shape of the triangular umbone of the elytra, to be ambiguous. He also illustrated the metasternal process and the parameres and explained in the key to subspecies that the ridge of the metasternal process does not reach the plate, and the process is weakly raised and more rounded anteriorly in D. a. drumonti. Also, the parameres of D. a. drumonti are shorter and with more acute lateral angles than of D. a. adamsi. However, we found that these characters were variable in the specimens from the two geographically isolated populations (Fig. 8). For example, the shape of the lateral angles of the parameres of Tibetan D. a. drumonti (Fig. 8C, D) is similar to that of a D. a. adamsi from South Korea (Fig 8K, L), and another specimen of D. a. drumonti from Sichuan, China (Fig. 8G, H) resembles a D. a. adamsi from Dandong, China (Fig. 8S, T). We did not find any significant diagnostic characters to separate the two subspecies and therefore the new synonymy is here proposed (Dicronocephalus adamsi drumonti Legrand, 2005 = Dicronocephalus adamsi adamsi Pascoe, 1863, syn. nov).

Discussion
From the results inferred from ML, BI, and MP methods using COI and 16S rRNA genes, the genus Dicronocephalus includes two major lineages, one with D. adamsi and D. yui yui and another with D. dabryi, D. uenoi katoi, D. w. bowringi, D. w. wallichii, and D. w. bourgoini (Figs 1-3). The specimens of eight taxa including subspecies clustered into seven groups and their monophyly was strongly supported in all analyses. However, D. w. bowringi was found to be paraphyletic and the monophyly of D. yui yui was not confirmed in the 16S rRNA based analyses. In the same analyses we also failed to identify the monophyly of D. yui yui (Fig. 3). Paraphyly or polytomy of the two species was the result of a few pasimony-informative nucleotide substitutions. This has a significant effect on phylogenetic reconstructions when the genetic divergences within and between species are low.
In all topologies, D. adamsi is sister to D. yui yui; the same was suggested by Kurosawa (1986). He grouped D. adamsi, D. shimomurai, and D. yui as the adamsi speciesgroup and mentioned that the female dark blackish body without markings might be the main characteristic of this group. The abdomen covered with whitish powder is also a trait that is only shared by D. adamsi and D. yui among the examined species (Pouillude 1914, Kurosawa 1986).
In contrast with the molecular data of the adamsi species-group, our results for the other congeners do not support the view of Kurosawa (1986). D. uenoi katoi is treated  as a separate group in his paper, but it appears a sister taxon of D. dabryi in our study, although the general appearance of D. uenoi katoi is rather similar to that of D. yui yui. Especially, these two species share two characters: the pronotal bands reaching the posterior border and the obtuse metasternal process. Pouillaude (1914) also noted that D. dabryi has tawny erect hair on the pronotum and elytra. We could observe that the pronotum and elytra are sparsely pilose and the hairs are much denser and longer on the ventral side compared with the other congeners. Furthermore, in the male genitalia, the parameres of the two species are similar and much shorter than those of other species. In this study, the pilose body, which is represented as a unique character of D. uenoi katoi by Kurosawa (1986), is considered as autapomorphy, which may have been rapidly acquired during allopatric speciation in Taiwan because D. uenoi katoi was isolated from a continental ancestor. This interpretation disagrees with Kurosawa's presumption that D. uenoi katoi is the most primitive in this genus.
Regarding the status of the subspecies of D. adamsi, Legrand (2005) recognized discontinued distribution and morphological differences between two geographically separated populations; however, we consider almost all of the diagnostic characters as being unsuitable for distinguishing these two subspecies. Furthermore, the molecular data indicates that the two subspecies form a monophyletic group with low genetic divergences (0-1.7%) and individuals of the both subspecies share haplotypes. Therefore, our results provide strong evidence that D. a. drumonti should be synonymized with D. a. adamsi.
The three subspecies of D. wallichii were originally described as separate species (Hope 1831, Pascoe 1863, Pouillaude 1914). Subsequently their status was lowered to subspecific (Paulian 1960, Mikšić 1971, 1977, Krajcik 1998, Sakai and Nagai 1998, Šípek et al. 2008, Young 2012, Krajcik 2014. However, Kurosawa (1968) disagreed with Paulian (1960) as he considered that there were significant morphological differences between them such as the characteristics of the antlers, the clypeus, the marginal carinae of the pronotum, and the metasternal process. Devecis (2008) also proposed that the taxa be restored as species based on the morphological differences such as color of the dorsal setation, shape of the antlers, and length of the pronotal bands. Results of our molecular analyses showed that the three subspecies of D. wallichii form a monophyletic group with high supporting values and large genetic distances. The average pairwise distances (4.7%-6.0%) of COI between D. wallichii bowringi + D. wallichii wallichii and D. wallichii bowringi + D. wallichii bourgoini. D. wallichii wallichii + D. wallichii bourgoini were slightly lower than the average inter-specific distances of D. adamsi + D. yui yui (6.2%) and D. dabryi + D. uenoi katoi (6.9%) ( Table  2). Also, in 16S rRNA analysis, the pairwise distances between the three subspecies of D. wallichii were similar to (0.8%-1.6%) the distance between D. adamsi and D. yui yui (0.8%-1.2%) ( Table 3). Our phylogenetic analyses explicitly explain their evolutionary history. D. w. bowringi is the most primitive among this group and D. w. wallichii might be separated by parapatric speciation in the continental region. Also, D. w. bourgoini might have undergone allopatric speciation after colonizing the volcanic island of Taiwan. Our results support specific rather than subspecific rank of the three members of D. wallichii. We revealed them as being in a monophyletic cluster (Mishler andTheriot 2000, Wiens andPenkrot 2002) with each other separated by distinct genetic gaps in the COI and COI+16S analyses, although not in the 16S rRNA analy-sis. Also, our study showed two distinguishable morphological characters, namely the color of the dorsal body side in males and the shape of the metasternal process (Table  5). However, this evidence is not strong enough to propose specific rank for each of them. A recent study showed that the high genetic divergence of COI alone cannot be a reason for species separation in Cetonia aurata aurata (Ahrens et al. 2013). There is a need for additional analyses with representative sample sizes and the use of multiple genetic loci to reconfirm our results.