Taxonomic reassessment of the genus Dichotomius (Coleoptera: Scarabaeinae) through integrative taxonomy

Sampling

Our total sample set consisted of 304 individuals of Dichotomius (31 species). The morphological analysis of male genitalia included 208 individuals from 28 species (Table S1), whereas the genetic analysis consisted of 145 specimens from 16 species; 52 of these sequences were obtained from GenBank (Table S1). This is representative of 14 species-groups and three subgenera in Dichotomius. All specimens for which we obtained data (DNA or morphology) came from the following biological collections: (i) Colección Alejandro Lopera-Toro (CALT-ECC, Colombian Collection ID 2), (ii) Museo de Historia Natural Universidad Distrital (MUD, Colombian Collection ID 46), and (iii) Colección de Artrópodos de la Universidad del Rosario (CAUR, Colombian Collection ID 229). These individuals were identified by experts or using most recent taxonomical keys (Nunes, 2017; Sarmiento-Garcés & Amat-García, 2014; Vaz-de-Mello et al., 2011).

Morphometric analyses

Because the shape of the male genitalia is considered one of the most informative morphological characters in the classification of Dichotomius species (López-Guerrero, 2005; Sarmiento-Garcés & Amat-García, 2014), we analysed the quantitative variation of the aedeagus in 208 individuals (28 species; Table S1). Male genitalia preparation followed a standard procedure: we detached the last two abdominal segments, soaked them in 10% KOH at 60 °C–70 °C for 12 h and neutralized them in 1% acetic acid to finally store them in glycerine (Sarmiento-Garcés & Amat-García, 2014). Then, we cleaned and dissected the aedeagus. Finally, we photographed the aedeagus in dorsal view and using a Leica DFC320 digital camera coupled to a Leica S6 stereoscope at 4X magnification.

We applied landmark-based geometric morphometrics to these photographs in order to analyse genital shape. We used tpsDig v.2.31 (Rohlf, 2004) to digitise 33 landmarks per individual that describe the outline of the aedeagus, all of them were placed on the parameres (Fig. S2A). This landmark dataset was subjected to superimposition using a Generalized Procrustes Analysis (GPA) in the R package ‘geomorph’ (Adams & Otárola-Castillo, 2013). For this, the software aligns, scales and rotates the configurations to line up the corresponding landmarks as closely as possible, minimizing differences between landmark configurations without altering shape. Then, we obtained partial warps (or shape variables) that indicate partial contributions of hierarchically scaled vectors spanning a linear shaped space. With this information we generated a consensus shape that summarises the aedeagus’ shape variation among all Dichotomius species included (Fig. S3). In this way, each specimen’s shape is quantified by the deviation of its landmark configuration from the average landmark configuration (i.e., consensus shape), which allows to visualise differences between groups. Differences in aedeagus’ shape among species were tested using a Procrustes MANOVA applied to the aligned landmark configurations. This was done using the procD.lm function in the ‘geomorph’ R package (Adams & Otárola-Castillo, 2013).

We implemented a principal component analysis (PCA) on the procrustes aligned data using the plotTangentSpace function in the ‘geomorph’ R package (Adams & Otárola-Castillo, 2013). Of the 66 PCs produced, the first two cumulatively accounted for ∼92% of the total shape variance; therefore, further analyses were performed on these PCs. We used the function plotRefToTarget from the same package to generate the deformation grids representing the extremes (maximum and minimum) of shape variation along the principal components 1 and 2 (PC1 and PC2). We then applied a discriminant analysis of principal components (DAPC) using the R package ‘adegenet’ (Jombart, 2008).

We also applied a model-based hierarchical clustering using the R package ‘mclust’ (Scrucca et al., 2016) in order to identify groups of individuals that resemble each other, independent of other evidence or a priori assignments. This method uses expectation maximization (EM) to estimate the Maximum Likelihood (ML) of alternative multivariate mixture models that describe shape variation in the data and estimates the optimal number of clusters based on the Bayesian Information Criterion (BIC). All models were evaluated for a predefined number of 1 to the maximum number of morphospecies studied (28 in our case, i.e., those for which morphology data was available).

Molecular analyses

We extracted DNA from legs of 95 specimens of Dichotomius using the DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany) following the manufacturer’s instructions with minor modifications: 40 µL of proteinase K were used, protein digestion lasted for at least 2 h and the final elution was made in 100 µL of warm AE buffer. Then, we amplified the 3′ and 5′ ends of the cytochrome c oxidase I gene (COI), and the nuclear gene 28S. All PCR reactions were performed in a final volume of 10 µL containing one µL of 10X Buffer, 0.6 µL of MgCl₂ (25 mM), 0.5 µL of dNTP mix (10 mM), 0.5 µL of each primer (10 µM), 0.05 µL of DNA polymerase (five U/µl; QIAGEN) and 5.85 µL of dH₂O. To amplify the 3′ end of the COI gene we used the primers C1-J-2183 (Jerry: 5′-CAACATTTATTTTGATTTTTTGG-3′) and TL2-N-3014 (Pat: 5′-TCCAATGCACTAATCTGCCATATTA-3′) (Simon et al., 1994). The amplification PCR profile consisted of an initial denaturation step of 94 °C for 5 min, 7 cycles of denaturation at 94 °C for 1 min, annealing at 48 °C for 45 s and extension at 72 °C for 1 min, followed by 33 cycles of denaturation at 94 °C for 45 s, annealing at 52 °C for 45 s and extension at 72 °C for 1.5 min, with a final extension at 72 °C for 10 min. The 5′ end of the COI gene (the barcode) was amplified with the primers LCO1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) (Folmer et al., 1994), using the following PCR conditions: 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 45 °C for 30 s, 72 °C for 1.5 min and a final extension at 72 °C for 7 min. To amplify the 28S gene we used the primers 28SFF (5′-TTACACACTCCTTAGCGGAT-3′) and 28SDD (5′-GGGACCCGTCTTGAAACAC-3′) (Monaghan et al., 2007). PCR cycling was 94 °C for 5 min, 38 cycles of 94 °C for 30 s, 53 °C for 30 s, 72 °C for 45 s and a final extension of 72 °C for 10 min.

All PCR products were purified with ExoSAP and their bidirectional sequencing was carried out by ELIM Biopharmaceuticals Inc. (Hayward, CA, USA). Forward and reverse sequences from each amplicon were verified and assembled into a single consensus contig based on a minimum match of 80% and a minimum overlap of 50 bp using CLC main workbench.

Sequences of each genetic marker were aligned independently using MUSCLE (Edgar, 2004) in MESQUITE v3.04 (Maddison & Maddison, 2011); poorly aligned regions were corrected manually. Protein coding sequences were translated into amino acids to confirm the absence of stop codons and anomalous residues in MESQUITE v3.04 (Maddison & Maddison, 2011). Additional sequences of Dichotomius available in GenBank (Table S1) were downloaded and integrated into the alignments. All sequences generated by us were deposited in GenBank and their accession numbers are listed in Table S1.

We estimated a phylogenetic tree based on the sequence information from the COI, 28S and 16S. All sequences from the latter marker were obtained from GenBank and correspond to the species D. nisus, D. yucatanus, D. parcepunctatus and D. boreus (Table S1). We concatenated all genes into a single alignment (2,546 bp) that included 16 species of Dichotomius and nine outgroups: Deltochilum larseni, Neateuchus proboscideus, Ontherus diabolicus, Pedaria sp., Panelus sp., Australammoecius occidentalis, Euphoresia sp., Brindalus porcicollis, Pleurophorus caesus (Table S1). We calculated a ML tree using IQ-TREE using the entire haplotype set derived from all species and individuals (Nguyen et al., 2015) with 1,000 ultrafast bootstrap replicates. This was done based on the substitution model showing the smallest AIC score for each partition (i.e., COI, 28S and 16S), which was also selected using IQ-TREE (Nguyen et al., 2015; Table S2).

To test whether D. satanas exhibits genetic clustering associated to the Colombian Cordilleras of the Andes as previously suggested (Sarmiento-Garcés & Amat-García, 2014), we also estimated a ML topology using all sequences available for the Colombian specimens of this species (COI and 28S) and using the conditions aforementioned. The sequences were all concatenated into a single alignment of 2,145bp that included one individual of D. boreus, D. quinquelobatus and D. protectus (outgroups) and 60 individuals of D. satanas: 8 from the Central Cordillera of Colombia, 14 from the West Cordillera of Colombia and 38 from the East Cordillera of Colombia.

Finally, we used DnaSP version 6.12.01 (Rozas et al., 2003) to calculate diversity parameters (i.e., number of haplotypes (H), haplotype diversity, genetic diversity (π and θ) and Tajima’s D) for all species and for D. satanas, as well as summary statistics of population differentiation among populations of D. satanas.

Species delimitation analyses

We implemented a joint Bayesian inference based on genetic and phenotypic data to delimit species using iBPP (Solís-Lemus, Knowles & Ané, 2014). This was done using two independent data sets: (i) all species, and (ii) D. satanas from Colombia only. In both cases, we ran the program for three different datasets: (i) morphological and molecular data combined, (ii) morphological data alone, (iii) molecular data alone. In all cases, we used the species-tree topology from IQ-tree as the guide tree. The morphological character matrix used as input included the values of PC1 and PC2 from the geometric morphometric analyses. The molecular matrix included all sequences available for the markers COI, 16S and 28S. We specified nine combinations of the prior distribution for the ancestral population size (θ) and the root age of the tree (τ) ranging from scenarios that represent large population sizes and a deep divergence time (θ = G(1, 10) and τ = G(1, 10)) to those representing small population sizes and a shallow divergence time (θ = G(2, 2000) and τ = G(2, 2000)) as previously used (Eberle, Warnock & Ahrens, 2016; Olave et al., 2017). We used default values of σ² and κ = 0, thus these priors are non-informative and the program estimates them. The MCMC analysis was run over 50,000 generations, sampling every 1,000 steps and using a 10% burn-in. We confirmed the robustness of the results by running the analysis with both the algorithms 0 and 1 for rjMCMC searches. As results were very similar, we present those of algorithm 1. The parameters of the locus-specific rates of evolution were fine-tuned using an auto option.

Morphological analyses

When we tested for aedeagus shape variation in the entire Procrustes shape space, we found differences among all categories tested (i.e., subgenera, species-groups and species; Procrustes MANOVA p < 0.001 in all cases). The PCA of the aedeagus shape revealed that most of its variation is contained in few dimensions. The first two PCs accounting for 91.9% of the total variance. PC1 explained 84.16% of the aedeagus shape and was driven by the width of the lateral outer margins in the apex of the parameres, ranging from being broad to narrow (Fig. 1A and Fig. S2). PC2 explained 7.7% of morphological of the aedeagus shape variation and describes the shape formed by the sides of the parameres (Fig. 1A and Fig. S2). The DAPC suggests the existence of four discrete genitalia morphology groups within Dichotomius (Fig. 1B and Fig. S4). The first group (depicted in red tones) was composed mostly by members of the subgenus Selenocopris sensu (Nunes, 2017) from the species-groups Agenor, Batesi and Inachus (i.e., D. agenor, D. batesi, D. belus, D. deyrollei, D. favi, D. fortestriatus, and D. yucatanus). This group also contained individuals of the subgenus Dichotomius s.s., exclusively those in the species-group Carolinus (i.e., D. amicitiae and D. coenosus). Finally, the species D. fonsecae (subgenus Cephagonus, species group Fissus) also clustered in this first group. The second group (depicted in green tones) was mainly formed by species that belong to the subgenus Dichotomius s.s. from the species-groups Boreus, Buqueti and Mamillatus (i.e., D. boreus, D. compresicollis, D. mamillatus, D. podalirius, D. riberoi and D. robustus); the species D. inachoides (subgenus Selenocopris, species-group Agenor) also grouped here. The third group (yellow) consisted exclusively of individuals from D. nisus (isolated species in the Selenocopris subgenus sensu Nunes, 2017). The fourth group comprised only species from the subgenus Dichotomius s.s., species-group Mormon, namely: D. alyattes, D. andresi, D. ohausi, D. protectus, D. divergens, D. quinquelobatus, D. quinquedens and D. satanas (blue tones). Although the species D. costaricensis and D. worontzowi (both of the Dichotomius s.s. subgenus) appeared well differentiated from any other species or group, we only have one sample for each of them, preventing us from making strong inferences. Consistently, mclust identified four clusters entirely coincident with the groupings obtained above (Fig. 1C). This variation is best explained by a model with ‘diagonal distribution, variable volume and equal shape’ (VEI; BIC = 1,152.184).

In summary, variation in genitalia morphology is not entirely consistent with the current taxonomy of Dichotomius (Nunes, 2017). Specifically, D. (Selenocopris) nisus (yellow) appears as different from other species in the subgenus Selenocopris (red). Also, species in the Carolinus group (D. amicitiae and D. coenosus), currently classified as members of Dichotomius s.s., cluster with species from the subgenus Selenocopris (red). Species in the subgenus Dichotomius s.s. formed two clusters, one that contains lowland species (green) and the other composed only by highland Andean species (blue).

Molecular analyses

We found Dichotomius as a monophyletic genus with two well-supported deep clades (Fig. 2, Fig. S5). The first clade contains D. (Selenocopris) nisus sister to species from the subgenus Dichotomius s.s. The second clade is composed of species from the Selenocopris subgenus, except for D. carolinus, which is currently included within Dichotomius s.s. Within the first clade (1 in Fig. 2), all species that belong to Dichotomius s.s. were grouped by species-group, with the Mormon, Boreus and Mamillatus groups forming each a monophyletic cluster (Fig. 2; Fig. S5). Within each of these species-groups most species appeared as monophyletic, except for D. satanas. This species formed two monophyletic clades, one consisting of Colombian specimens and the other composed by Central American individuals (Fig. 2; Fig. S5). Within the second monophyletic clade (2 in Fig. 2) we could not test the monophyly of the species groups due to low species sampling, yet we observed the Agenor species-group as paraphyletic (Fig. 2; Fig. S5). In general, mtDNA showed higher haplotype diversity than the 28S nuclear gene (Table 1).

When populations of D. satanas from Colombia were analysed separately to evaluate whether this species displays genetic clustering associated with geography or phenotype (Sarmiento-Garcés & Amat-García, 2014), we mainly observed clustering and genetic differentiation associated to the three Cordilleras of the north of the Andes (Fig. 3, Table 2). Individuals from the Central and the Western Cordilleras were reciprocally monophyletic, and both were sister to the Western Cordillera clade. Interestingly, this phylogenetic pattern associates to morphological differences in the females: the Central and Western clusters contain females with only two protuberances in the pronotum, while the cluster of the Eastern Cordillera includes females with two and four protuberances. At the same time, the latter cluster separates into two inner groups, one that contains only females with four protuberances and the second, where females of two and four protuberances are found (Fig. 3).

Species delimitation

The total-evidence (morphology and DNA) approach to Bayesian species delimitation (iBPP) did not support the a priori morphospecies assignment (Fig. 4). In most θ and τ scenarios tested, the posterior probability for the existence of the 16 morphospecies evaluated was lower than 50%. The only a priori defined species that consistently presented high support for all prior combinations were D. belus, D. nisus and D. mamillatus. Other species were supported only when modelling small population sizes (θ = 0.01) and medium to deep divergence time (τ = 0.05 and τ = 0.1), but never when modelling a shallow divergence time (τ = 0.01; Fig. 4).

The existence of two deep clades was strongly supported, regardless of the θ and τ priors used (1 and 2 in Fig. 4). In the first clade, the existence of species groups Nisus, Mamillatus, Boreus and Mormon was also stronly supported (Fig. 4). In the latter group, the separation of D. quinquelobatus from other members of this group showed high posterior probability values in most scenarios, except for those with τ = 0.01. However, the separation of D. protectus from D. andresi, or Colombian D. satanas from D. alyattes was rarely supported (Fig. 4). This was also observed in the Boreus species-group, where the delimitation between D. boreus and D. podalirius always had low posterior probabilities (Fig. 4).

In the second clade the only well supported species across all the parameter combinations (θ and τ) was D. belus. The later suggest that this species may not be part of the Agenor species group. In contrast, species within the clade sister to D. belus showed very low support in most of the parameter space (Fig. 4). The species delimitation based on molecular or morphological data alone were consistent with the total-evidence approach (Fig. S6). However, the results of these independent data types tended to provide stronger supports to species-groups and some species, especially the molecular data.

Finally, the total-evidence analysis of species delimitation done in D. satanas failed to identify any of the phylogenetic clusters associated to geography as separate species (in most θ and τ scenarios tested the support for these clusters was lower than 60%, Fig. S7A). This suggests that D. satanas is likely a single species with phenotypic polymorphism. However, just as before, the analyses with only molecular data presented stronger supports while the analysis based on morphological data provided very poor support Figs. S7B and S7C).