This study analyzed 1447 species including 1303 lucanid species and 144 outgroup species. All four subfamilies (Aesalinae, Lampriminae, Lucaninae, and Syndesinae), all 105 lucanid genera, and 1303 lucanid species (approximately 99.2% of all described lucanid species) from around the world were included as the inner group, and the outgroups consisted of 4 Diphyllostomatidae species, 16 Hybosoridae species, 16 Geotrupidae (Geotrupinae+Bolboceratinae) species, 11 Passalidae (Passalinae+Aulacocyclinae) species, 6 Glaresidae species, 19 Ochodaeidae species, 43 Scarabaeidae species, 9 Trogidae species, 12 Silphidae species, and 8 Histeridae species. The measurements of all lucanid species and most outgroup species were based on published images (photographs or specimen drawings) (Dallwitz 1980; Dallwitz et al. 1993, 1995, 2000; Sakai and Nagai 1998; Bunalski 1999; Král 2001; Arnett et al. 2002; Ocampo 2006; Mondaca and Smith 2008; Schenk 2008; Nikolajev 2009; Bomans 2010; Fujita 2010; Imura 2010, 2011, 2012; Boilly 2011; Kobayashi and Matsumoto 2011; Okajima and Araya 2012; Okuda 2012a; Palestrini et al. 2012; Paulsen and Ocampo 2012; Ballerio 2013; Král et al. 2013; Paulsen 2013; Bezborodov and Koshkin 2014), and those of the Passalidae species were based on the specimens housed in the Institute of Zoology of the Chinese Academy of Sciences (Suppl. material 1: Table S1). In consideration of sexual dimorphism and male polymorphism in Lucanidae, images of the male specimens were selected in this study, as they are more accessible in publications, and medium-sized specimens were chosen if there were images of different body sizes available.

Geometric morphometric analysis of the variations in the pronotum and elytron were based on one curve for each structure (Fig. 1), and the curves for the pronotum or elytron were resampled by length after 25 and 50 semi-landmarks (SLM), respectively. The curves were digitized with TPS-DIG 2.05 (Rohlf 2006), and the data file was modified as .txt file to convert the semi-landmarks to landmarks (MacLeod 2017). The two lines with the curve number and point number were deleted, and the landmark number was replaced by the point number. This approach was used earlier by Bai et al. (2014a) and Li et al. (2016). The landmark configurations were scaled, translated, and rotated against the consensus configuration by using the Procrustes superimposition method (Bookstein 1991). The differences in the shapes and diversity indexes of the pronotum and elytron were inferred on the basis of principal component analysis (PCA) in MORPHOJ 1.06a (Klingenberg 2011) (Figs 2, 3, 8, 9, Table 1, Suppl. material 1: Table S10). The diversity index was quantified as Procrustes variance, which measures the dispersion of all observations around the mean shape of the respective taxa (Zelditch et al. 2004; Sherratt et al. 2014). The association of morphological diversity with species richness at genus-level was measured using Pearsons correlation coefficient, r in PAST 3.01 software (Sherratt et al. 2014; Hammer et al. 2001) (Suppl. material 1: Table S10, S11).

Figure 1. 

Description of the curves used in geometric morphometric analysis. The positions selected for the pronotum and elytron curves are represented by Prosopocoilus sp. in dorsal view. The curves were resampled in 25 or 50 semi-landmarks (SLM).

A phylogenetic tree was visualized in MESQUITE 2.72 (Maddison and Maddison 2011) on the basis of earlier molecular analysis (Kim and Farrell 2015), and the aligned landmark data were entered into MESQUITE 2.72 as a continuous matrix and linked to the tree (Figs 4, 5). Because the branch lengths (Grafen 1989) were missing, we followed the evaluation proposed by Klingenberg and Marugán-Lobón (Klingenberg and Marugán-Lobón 2013) and assigned an equal length to all branches (i.e., an evolutionary model with the same expected amount of morphological change on every branch was assumed).

The ancestral groundplans of the Lucanidae pronotum and elytron were reconstructed by combining the landmark data with the phylogenetic tree, and the ancestral groundplans of all nodes were reconstructed by using the trace-all-characters and/or landmark-drawing modules of the RHETENOR package in MESQUITE. The ancestral states of all nodes were calculated and exported, and the data computed for the nodes were integrated with the original landmark data for the two characteristics from the 1303 stag beetles in EXCEL and NTSYS-PC (Rohlf 2007), respectively. The thin-plate splines showing the deformation of the landmarks compared with the original computed by MESQUITE were mapped onto the phylogenetic tree (Figs 4, 5).

In this case, the differences in the shapes of the pronotum and elytron among extant and extinct Lucanidae were inferred on the basis of PCA in MORPHOJ 1.06a and PAST 3.01 (Figs 6, 7). The canonical variate analysis (CVA) and discriminant function analysis (DFA) of the landmark data were based on MORPHOJ 1.06a (Suppl. material 1: Tables S2–S9, S12, S13).

The first two principal components of the pronotum and elytron from all 1447 species accounted for 77.37% and 88.40% of the variation among the species, respectively. The first two principal components were plotted to indicate variation along the two axes, which provided 90% equal frequency ellipses containing approximately 90% of the data points of each group (Figs 2, 3). The pronotum morphologies of the outgroups were mostly within the entire morphological variation of the Lucanidae (Fig. 2).

Figure 2. 

Differences in pronotum shape between outgroups and Lucanidae, on the basis of principal component analysis at the species level. The four circles are 90%-equal frequency ellipses of Lucanidae subfamilies.

Figure 3. 

Differences in elytron shape between outgroups and Lucanidae, on the basis of principal component analysis at the species level. The four circles are 90%-equal frequency ellipses of Lucanidae subfamilies.

All the p-values obtained from the permutation tests (10000 permutation rounds) for both Mahalanobis distances and Procrustes distances between any two Lucanidae subfamilies were less than 0.05 for both the pronotum and elytron. Most of the p-values for the pronotum and elytron Mahalanobis or Procrustes distances between the Lucanidae subfamilies and outgroups were less than 0.05, except for some of the distances between the Aesalinae or Lampriminae and the outgroups.

There were significant differences in both the pronotum and elytron between any two of the Lucanidae subfamilies (Suppl. material 1: Table S3, S5, S7, S9). For the pronotum, there was a significant difference between Lampriminae, Lucaninae, or Syndesinae and the outgroups; the Lucanidae and the outgroups partly overlapped, because the pronotum morphology of Aesalinae could not be distinguished from that of Geotrupidae, Glaresidae, Ochodaeidae, and Histeridae (Suppl. material 1: Table S3, S5). However, the differences in the pronotum among these four pairs were not equivalent. On the basis of the Procrustes distance, a measure of the absolute magnitude of the deviation in shape that indicates the extent of the differences between the average group shapes, the differences in the pronotum between Aesalinae and Geotrupidae, Glaresidae, Ochodaeidae, and Histeridae were 0.0512, 0.0637, 0.0488 and 0.0659, respectively (Table S4). For the elytron, there was a significant difference between Lucaninae or Syndesinae and the outgroups, but the morphology of Aesalinae and Lampriminae could not be distinguished from the morphologies of Glaresidae, Scarabaeidae, or Trogidae (Table S7, 9). There was greater overlap in the pronotum morphology of the Lucanidae subfamilies and outgroups than the elytron morphology (Suppl. material 1: Table S2, S6).

In a comparison of the ancestor of Lucanidae and all outgroups (node 2 in Fig. 4), the anterior angle of the pronotum of the ancestor of Lucanidae (node 8 in Fig. 4) is more obtuse (landmarks 4–12 at the splines in Fig. 4). It is shown in Figure 6 that a clear divergence between the two lineages of Lucanidae, the Aesalini lineage and lineage I (node 9 in Fig. 4), is primarily in the direction of the first principal component. Moreover, within lineage I, there is continued diversification in the direction of the first principal component (the horizontal direction in Fig. 6). The pronotum morphology of Lampriminae and Ceruchini showed clear changes in different directions, almost in the inverse direction of the second principal component (the vertical direction in Fig. 6), particularly in the anterior angle (landmarks 4–12 at the splines in Fig. 4).

Figure 4. 

Reconstruction of ancestral groundplans of the pronotum in Lucanidae and the outgroups. The splines indicate deformation of the shapes relative to the reference configuration. The phylogenetic tree was summarized and reconstructed from earlier molecular results (Kim and Farrell 2015).

The elytron of the ancestor of Lucanidae (node 8 in Fig. 5) is more slender than that of the outgroups (node 2 in Fig. 5) (the ratios of the length to the greatest width are 1.571 and 1.235, respectively) and narrower at the end (landmarks 35–50 at the splines in Fig. 5), but it is much wider than that of Diphyllostomatidae (the ratios of the length to the greatest width are 1.571 and 2.170, respectively). The elytron of Aesalini is wider than in lineage I (node 9 in Fig. 5) (the ratios of the length to the greatest width are 1.313 and 1.571, respectively), thus illustrating the major differences in the elytron between the two Lucanidae lineages. The humeral angle is narrower in the ancestor of Sinodendrini and Nicagini (node 15 in Fig. 5) than in the ancestor of Lucaninae, Lampriminae and Ceruchini (node 10 in Fig. 5) (landmarks 10–18 at the splines in Fig. 5).

Figure 5. 

Reconstruction of ancestral groundplans of the elytron in Lucanidae and the outgroups. The splines indicate the deformation of the shapes relative to the reference configuration. The phylogenetic tree was summarized and reconstructed from earlier molecular results (Kim and Farrell 2015).

Figure 6. 

Differences in pronotum shape among each branch and ancestor, on the basis of principal component analysis. Empty dots indicate the number of the node on the phylogenetic tree; solid dots indicate the average shape of the extant subfamily/tribe of each branch.

In terms of the species richness, Lucaninae is the largest subfamily of Lucanidae, comprising more than 90% (1220) of the species of stag beetles, followed by Aesalinae, Syndesinae, and Lampriminae, which have fewer than 50 species each (47, 25, 11). However, morphological diversity of the pronotum and elytra shape does not correspond with species richness (Table 1). Aesalinae, with fairly low species richness, was found to be the most diverse subfamily in both pronotum and elytron morphology, whereas Lucaninae exhibited low morphological diversity (both in the pronotum and elytron). This is an unexpected result considering its extremely high species richness.

In 73 genera (all Lucanidae genera with more than one species), similarly to the subfamily data, the morphological diversity of neither the pronotum nor the elytron was consistent with species richness. There is no significant correlation between morphological diversity and species richness in pronotum at genus-level (Procrustes variance r = 0.15, P =0.21), nor a relationship between morphological diversity and species richness in elytron at genus-level (Procrustes variance r = 0.11, P = 0.33). The total morphological variances in extremely species-rich genera, such as Aegus, Dorcus, Lucanus, and Prosopocoilus, were not predominant.

Figure 7. 

Differences in elytron shape among each branch and ancestor, on the basis of principal component analysis. Empty dots indicate the number of the node on the phylogenetic tree; solid dots indicate the average shape of the extant subfamily/tribe of each branch.

According to the Mahalanobis distances from the DFA (Suppl. material 1: Table S12, S13), all lucanid ancestors (nodes 8–10, 12, and 15) most resemble Nicagini in their pronotum and elytron morphology. Procrustes distances indicated that the pronotum of the ancestor of Lucanidae (node 8) is closest to that of Scarabaeidae, and the elytron is closest to that of Glaresidae. The common ancestor of Lucaninae, Lampriminae, Syndesinae, and Nicagini (node 9) most resembles Hybosoridae in the pronotum and Lucaninae in the elytron, and the ancestor of Lampriminae and Ceruchini (node 12) most resembles Hybosoridae in the pronotum and Sinodendrini in the elytron. The ancestor of Lucaninae, Lampriminae, and Ceruchini (node 10) is most similar to Lucaninae in both the pronotum and elytron, and as shown by both Mahalanobis and Procrustes distances, the ancestor of Sinodendrini and Nicagini (node 15) most resembles Nicagini in the pronotum as well as the elytron.

The ancestral pronotum and elytra were reconstructed, which could be combined with fossil materials to uncover the ancestors’ habitat as well as evolution procedure. There is still a lack of sufficient data as well as studies of the functional morphology of the pronotum and the elytron in Lucanidae, but the functional morphology of other insect clades may allow for certain interpretations. Broad pronotum and elytra of stag beetles may provide advantages during locomotion and hunting prey, like ground beetles (Evans 1977; Forsythe 1981). The diversity of pronotum may reflect occupancy of diverse habitats and niches, as in the case of grylloblattids (Bai et al. 2010).

Figure 8. 

Species richness and morphological diversity of the pronotum at the genus level.

Morphological and taxonomic diversity provide insight into the expansion and contraction of major taxa, and the nature of the relationship between these two aspects of diversity has important implications in evolutionary mechanisms (Foote 1993). Lampriminae has the lowest pronotum and elytron diversity as well as species richness, but species richness and morphological diversity do not always vary consistently in Lucanidae subfamilies. Despite having many more species than the other three subfamilies combined, the subfamily Lucaninae has little morphological diversity. Furthermore, the results from the genus-level data reveal the same pattern in which some of the largest genera are not as diverse as the small groups in terms of either pronotum or elytron morphology.

The theory that species richness generates a variety of forms has been tested and supported in various studies (Williams and Humphries 1996; Roy and Foote 1997; Bell and Barnes 2001), thus suggesting that the relationship between morphological and species diversity should be monotonic or at least positive.

In contrast, numerous studies have indicated that richness and morphology do not always follow a common trend. Foote (1993) has found that morphological variety and taxonomic richness often increase together during the initial diversification of a clade, but two major patterns have been observed as clades decline. In Blastoidea, Trilobita, Libristoma, and Asaphina, morphological diversity continued to increase even in the face of striking decreases in taxonomic richness, but in Phacopida, Scutelluina, and to some extent in Proetida, morphological diversity decreased along with taxonomic diversity. Roy et al. (2001), using data from a large group of IndoPacific gastropods (family Strombidae), have shown that the species richness of a region is a poor predictor of morphological diversity. Areas with only a few species may harbor an impressive array of morphologies, and in contrast, morphological diversity in the most species-rich regions is no higher than that in regions with half the taxonomic diversity. Bell and Barnes (2002) have found no significant correlation between species richness and morphological diversity in cave or boulder habitats, although these variables are significantly correlated in coral reef and soft substratum habitats. In another scenario (Mohedano-Navarrete et al. 2008), the morphological diversity and species richness of Porites corals has been found to vary independently; some regions with few species had remarkably high morphological diversity, including peripheral areas such as Polynesia and East Africa. In a study of 107 families of passerine birds, morphological space was weakly related to the number of species in a family, i.e., the higher species richness of the order Passeriformes in the tropics compared with temperate regions was not matched by increased morphological diversity (Ricklefs 2012).

Figure 9. 

Species richness and morphological diversity of the elytron at the genus level.