Relationships between species richness and morphological diversity: insight from a geometric morphometric analysis of scarabs

Yijie Tong University of Chinese Academy of Sciences; Institute of Zoology, Chinese Academy of Sciences Haidong Yang Institute of Zoology, Chinese Academy of Sciences; Guangdong Institute of Applied Biological Resources Xingke Yang Institute of Zoology, Chinese Academy of Sciences; Guangdong Institute of Applied Biological Resources Ming Bai (  baim@ioz.ac.cn ) https://orcid.org/0000-0001-9197-5900


Results
In this study, 5148 scarabs covering 92% families and 70% scarabaeoid subfamilies of Scarabaeoidea were selected as a test group. The taxonomic categories of family and subfamily were selected for phylogenetic representation and four feeding types (omnivory, phytophagy, coprophagy and non-feeding) were selected as ecological factors. A geometric morphometric analysis on the pronotum and elytron, which represent the major aspect of the morphology in dorsal view, were conducted. Correlations between MD and SR among sub-/families were always found to be consistent, except a genus-level test of the pronotum on the subfamily category and a species-level test of elytron on the family category. These phenomena re ected a strong phylogenetic in uence. However, morphological diversities were not correlated signi cantly with species richness for any feeding type (ecological factors). In addition to this hypothesis test, the morphological diversity and known world species richness (~ 32000 species and 2453 genera of Scarabaeoidea) were estimated as well. Our global results were very similar to our sampling results.

Conclusions
Based on our analyses the hypothesis of ecological factor (represented by feeding habit) correlations was rejected in scarabs. On the contrary, phylogenetic factor correlations were found to be statistically signi cant. The morphological diversity of scarabs in most sub-/families tted their species richness. At the same time, the correlation our test parameters revealed variation in both groups and characters, which might be caused by morphological changes under coevolution with different ecological factors.

Background
Biodiversity plays a host of important roles in biosphere operations. This factor is usually quanti ed via proxy indicators, such as morphological forms and genes, and used to study changes in ecosystems. Use of different biodiversity indicators highlights the disparate nature of biodiversity investigations [1,2]. Among these, morphological diversity and species richness have often been used to explore patterns of biodiversity in different ecosystems [3,4]. Morphological diversity re ects biological feedback and ecosystem functionality [5]. Species richness reveals the number of species in a community or area.
Because both morphological diversity and species richness contain certain meaningful information, analysis of the relation between these two indicators has been regarded widely as a matter of importance [6][7][8].
Most studies have revealed that morphological diversity is positively correlated with species richness in many groups (from unicellular organisms to most animal groups, plant, etc.), which is presumed to re ect the in uence of phylogenetic factor [9][10][11][12]. However, if ecological factors are considered tests of the relation between morphological diversity and species richness have not always yielded consistent results [13][14][15]. The inconsistent conclusions of previous studies remind us that the relation of these two indicators is not simple. New approaches and larger datasets might provide greater insight into this problem.
In this study, Scarabaeoidea was selected as the test group owing to its worldwide distribution and complex habits, which set up their diverse external forms and behaviors [16][17][18]. Two samples were selected to represent the species/genus richness of Scarabaeoidea. For the larger sample, all known species/genera of Scarabaeoidea were included. However, to date it has not been possible to collect morphological diversity data for all ~ 32 000 species known globally. To facilitate a better understanding the problem, a large number of scarab species (5148, c. 1/6 of known species) were selected. This sample included 11 scarabaeoid families (92% of living families) and 26 subfamilies (70% of living subfamilies). All specimens were identi ed and analyzed in two ways, re ect patterns of phylogenetic relation (family and subfamily) and ecology (four feeding types: omnivory, phytophagy, coprophagy and non-feeding). Thus partitioned, these data were used to explore the correlation between morphological diversity and species richness.
Geometric morphometrics, a useful approach on the quantitative morphology analysis, was applied in this study. Two character complexes -the pronotum and elytron -which represent major parts of the body structure in dorsal view, were employed as the foci of character tests. The former supports prothorax's muscle system and prothoracic legs' locomotion [19] and the later represents an autapomorphy of Coleoptera [20]. These characters differ from the mouthparts with their strong selective pressures from food types, the pronotum and elytron carry a large amount of more generalized morphological information (e.g. their outer contours) that are unde ned from the standpoint of ecological selection pressure(s).

Results
Species/genus richness among groups Species/genus richness values from the four feeding types and 37 taxonomic groups (11 families and 26 subfamilies) are illustrated in Fig. 1 (see also Suppl. 1). Phytophagous species exhibited both the highest species richness (3651) and genus richness (575) values while non-feeders were found to exhibit the lowest species richness (29) with only 2 genera. The omnivorous and coprophagous groups exhibited similar species richnesses (370 and 369, respectively), while the genus richness of the former was higher than the later (62/36). In the family category, for our sample Scarabaeidae and Lucanidae both exhibited very high species richness (2776/1305) and genus richness (534/81) values primarily because of their large and diverse sub-groups; Glaphyridae, Glaresidae and Pleocomidae exhibited very low species (29/29/26) and genus (2/1/1) richnesses (with the species richnesses of the Ochodaeiade and Passalidae being higher than those of all three families mentioned above; 36/33), as well as higher per capita genus richness values (9/14). In addition, the Diphyllostomatidae was found to exhibit the lowest species richness value (3) with only one genus included in the test sample. Similar patterns were also found in subfamily category test, with most subfamilies exhibiting high species richness and relatively high genus richness values (Fig. 1).

Morphological Variations Of Pronotum And Elytron
The rst three principal components (PCs) of the pronotum and elytron from 5148 species accounted for 87.90% and 93.36%, of observed shape variation respectively (Fig. 2, 3). Shape models were calculated at equally spaced intervals along these PC axes to document the deformation of test characters, with the overlaying of these shape models used to reveal the comparative along-axis shape-change trends (Fig. 2, 3). Variation summaries for the other principal components from this analysis are provided in Supplement 2.
Along the positive direction of the rst PC axis, pronotum shape changed from a trapezoidal to inverted trapezoidal outline with the producing anterior angle, and retracting anterior margin, making the anterior angle sharper, while the posterior angle diminished and became less distinct. Along the positive direction of the second PC axis, the pronotum became stretched longitudinally and diminished horizontally with both anterior and posterior angles diminishing and becoming less distinct and the entire pronotum becoming markedly more oval. Along the third PC axis, the entire pronotum tends toward an inverted triangle with the anterior angle extending outward horizontally, the posterior angle shrinking inwardly, and the posterior margin stretching longitudinally.
Along positive direction of its rst PC axis, the elytron becomes laterally narrower with the anterior angle contracting inward. The shrinking smaller scutel causes the anterior margin to become shorter with the lateral margin, that protruded outward originally, contracting inward. Along the positive direction of the second elytron PC axis the anterior angle becomes larger and the posterior margin protrudes. The entire elytron stretches longitudinally and shrinks horizontally. Along the positive direction of third elytron PC axis, the smaller scutel attens the front margin. At the same time, the posterior margin shrinks inward to make its trace a smooth curve.

Morphological Diversity Among Groups
Morphological diversity of both the pronotum and elytron among the 11 families and 26 subfamilies re ects the in uence of phylogenetic factors. The Scarabaeidae exhibited both the highest morphological diversity ( Fig. 1) and species richness (Suppl. 1). The Diphyllostomatidae, a small group for which only 3 species were included in this study, was found to have the lowest relative morphological diversity (Fig. 1). However, the relationship between morphological diversity and species richness was not always consistent. The Lucanidae exhibited high species-richness values despite its fairly low morphological diversity (Fig. 1). The Hybosoridae, with a comparatively small number of test species (57), was found to have the second highest pronotum diversity.
In the test of 30 subfamilies ( Fig. 1) (Suppl. 1), results similar to the family test were obtained. Morphological diversities of certain subfamilies were also found to lack a strong association with species richness. The Aesalinae were found to have the highest pronotum diversity, despite their low species richness (47) in our sample. Similarly, with 263 species included in our sample the, Melolonthinae was found to have fairly low elytron diversity.
In addition, morphological diversity was also found to vary with the character of features tested. The changing trends of characters and species richness values were not always highly associated (Fig. 1). For example, elytron diversity of the Lucanidae was low despite its high pronotum diversity. In the subfamily test, the Aesalinae were found to exhibit a high pronotum diversity, but a fairly low elytron diversity.
Moving on to the ecological factors ( Fig. 1), phytophagous species were found to exhibit both the highest morphological diversities and guild richness to comprise more than 70.92% of all included species (Suppl. 1). The non-feeding group, containing two families (only 29 species and 2 genera), was found to exhibit the second highest pronotum diversity. However, as with the previous analysis, the relation between morphological diversity and feeding-guild richness was not always consistent. Coprophagous taxa exhibit high guild richness values despite their fairly low morphological diversity. Relative to be far more numerous phytophagous species, omnivorous species were found to exhibit the highest elytron diversity. At the same time, we found the changing trends of characters tested and species richness values also displayed inconsistent associations. The phytophagous group was found to have the highest pronotum diversity, however, its elytron diversity was not unusually high.

Relationship Between Morphological Diversity And Species Richness
Based on ordinary least squares (OLS) regression analysis, an ANOVA (F-test) was used to test the signi cance of the regression between morphological diversity and species richness quantitatively, and the equations of the best-t lines estimated to verify the linear relations (Fig. 4). Test results varied both by groups and traits (Suppl. 3).
Three level of relations were revealed by the phylogenetic factor examined in this study. Firstly, the correlation between morphological diversity and species richness among sub-/families was found always to be consistent, with two exceptions: the genus level test of pronotum on subfamily (p-value = 0.188) and the species level test of elytron on the family category (p-value = 0.089), the regressions of which were non-signi cant statistically. These results re ect the in uence of phylogeny on the mediation of morphological shape in these two broadly construed regions of the phenotype. Correlation of the family category with pronotum shape was signi cant statistically (p-value = 0.005 for species; p-value = 0.029 for genera) with the subfamily category test being marginally non-signi cant (p-value = 0.049 for species). With regard to the elytra, the comparative result was the opposite in all tests returning highly signi cant statistical results (p-value = 0.043 for genera in family category; 0.006 in species and 0.020 for genera in subfamilies).
Application of the non-parametric Spearman correlation coe cient to our test parameters revealed variation in both groups and characters. For the family test, the correlation between test parameters was found to be higher than it was for the subfamily both for the pronotum test (r = 0.5321/0.7335 for genus/species level test in family category and 0.1755/0.3397 for genus/species level test in subfamily category test respectively) and for the elytron test (r = 0.4587/0.5968 in genus/species level test in family category and 0.3898/0.5565 in subfamily category). At the same time, linear modeling in the specieslevel test was found to be lower than it in the genus-level test for both character complexes.
In addition, morphological diversity and known global species richness values for these same categories (11 families and 26 subfamilies include 32000 species and 2453 genera), were calculated. Based on a covariance analysis no signi cant difference was found in the slopes and intercepts of tting linear models between the global and sampling (p > > 0.05) ( Table 1), and the correlation between morphological diversity and species/genus richness was same as in the sampling test. This suggests our analysis represented the correlation between species/genus richness and morphological diversity objectively. For the family test, this correlation was found to be higher than it in the subfamily in the pronotum (r = 0.6800/0.7130 for genus/species level test in family category and 0.04260/0.2460 for genus/species level test in subfamily category test respectively) and in the elytron (r = 0.7130/0.7710 in genus/species level test in family category and 0.6200/0.4020 in subfamily category). The species-level test was found to exhibit a higher correlation than in the genus level test in family category, on the contrary the species-level test was found to exhibit a lower correlation than in the genus level test in subfamily category (r = 0.4020/0.2460 in species level test and 0.6200/0.4260 in genus level test).
However, morphological diversities were found to be insigni cantly correlated with species richness in all feeding types based on linear regression analyses (ecological factor) (p > > 0.05) (Suppl. 3, Fig. 4).

Discussion
In this study, we compared the species richness and morphological diversity in the pronotum and elytron of ~ 32000 scarab species using geometric morphometric and regression analyses. The in uence of phylogenetic and ecological factors on the correlation between morphological diversity and species richness was examined quantitatively. Results showed that the phylogenetic factor correlations were found to be statistically signi cant. On the contrary, the hypothesis of ecological factor (represented by feeding habit) correlations was refuted by the lack of statistical signi cance revealed by our tests.
This is the rst study to examine the relationship between morphological diversity and species richness in all known species of Scarabaeodiea. The role of phylogenetic and ecological factors on the associations between pronotom and elytron shape variation with taxonomy and an important aspect of ecology were revealed in the context of quantitative analyses. Only two characters that are unde ned from the standpoint of ecological selection pressure(s) were selected in the study. However, there are many more potential candidate characters available for study not in scarabs and in other beetles. In future studies more traits and samples from different groups of organism should be added to improve our understanding of the relation between morphological diversity and species/guild richness.
Different morphologies are associated with diverse functional aspects of niches [21,22] and both the evolution and development of taxa. These account for the variable morphological characteristics of scarabs [18,23]. When compared to feeding habit, phylogenetic factors were found to explain observed correlation between morphological diversity of the pronotum and elytron and the species richness of scarabs in a certain extent. In this study, complex and irregular relations existed in different groups on the entire family-category and subfamily-category by scarabs' diverse habits and the quite different coevolution of test characters. Scarabaeidae was found to exhibit the highest species richness and morphological diversity, a result that can be explained by the variety of its habitats -from rainforest to desert -and its taxonomic richness [24]. Low elytron shape diversity was found in the test of Passalidae, which are adapted to live in the tunnels of rotten wood. A total of 3 species of Diphyllostomatidae have been recorded worldwide, its small population and degraded characters determined that this group was relatively conservative [16], leading to its very low external morphological diversity. The Hybosoridae is a large, diverse and globally distributed family of c. 230 species and 30 genera. Its diversi ed food habits, which include fresh corpses, small insects or fungi [25], is associated positively with its high pronotum diversity of morphology despite only 57 species being included in this study. Meanwhile, based on our regression analysis results, correlation of morphological diversity to species richness in subfamilycategory test was found not as good as that exhibited by the family category. In this study, our inconsistent results showed that the factors affecting these correlations are complex and diverse with phylogenetic factor exhibiting high correlations between species/genus richness and morphological diversity to some extent. Additionally, some of the high ranking traditional taxonomic categories (e.g. sub-/families) are arti cial phylogenetically and so were not formally tested in the phylogenetic analyses [20,24,26]. In this case, the results of this study could be slightly changed if the classi cation system of Scarabaeoidea be modi ed in future. It would be interesting to make additional test following the stable classi cation of Scarabaeoidea in future based on the comprehensive sampling and phylogenetic analyses.
Biological diversity can be quanti ed in terms of species richness and morphological diversity [27], both of which are dependent on ecological variability and species interactions [28,29]. As a result of selectionbased feedback due to the variability of ecological factors, new forms emerge as a result of biological evolution [30]. In our study, feeding habit was found not to be a major factor affecting the correlation between scarab pronotum and elytron morphological diversity relative to this group's variation in species richness, complex habits and coevolution. This phenomenon can be found in some other groups as well. For example, the special galeas of stag beetles adapt to the relatively xed food, but their conservative diets don't affect the diverse morphological changes of their pronotums which play an important role in their ability to raise their heads while raising opponents with their mandibles [31]. More conservative external morphologies and decreased morphological diversities have been found in coprophagous and non-feeding scarab groups [32], presumably because of their specialization and lack of diet diversity, while the variety of the external forms that characterize phytophagous and omnivorous groups arise as a result of coevolution with their diverse diets. These well-understood principles account for the shape diversity exhibited by scarabs' pronotum and elytron in coprophagous, non-feeding, phytophagous and omnivorous groups. In addition, our sample of phytophagous scarabs was found to exhibit a very high pronotum diversity, but fairly low elytron diversity, which implies that the pronotum is under more direct selection pressure from diverse environment and surrounding ecological factors, presumably as a consequence of playing an important role the structuring of the muscle systems responsible for leg locomotion. The elytron's role is mainly one of protecting the beetles' abdomen, and so its requirement for environmentally mediated shape adjustment is far less pronounced than that of pronotum.

Conclusion
As a large and diverse group, the Scarabaeoidea was found to have a varied morphological appearance of pronotum and elytron. Procrustes PCA revealed that this variation re ects associations with various constraining factors. Coprophagous and non-feeding groups exhibit specialized external forms despite their fairly stable habits. Feeding habit was found not to be associated with morphological diversity or species richness, while phylogenetic factors were found to reveal strong associations with species/genus richness and with morphological diversity. In uenced by niche and coevolution, the species/genus richness and morphological diversity values did not always display consistent family-category and subfamily-category relations.
Nonetheless, there remain many limitations in our results. Our phylogenetic factor and feeding habit analyses revealed only a small fraction of the impact these factors may have on the correlation between species/genus richness and morphological diversity. In addition, present limitations on the composition of our test sample for dimorphic groups makes our results less generalized, and perhaps less stable, than they might have been otherwise. These issues can be addressed when more complete samples become available.
In future studies more traits should also be considered to improve our understanding of morphological diversity and richness, and more parts or ecological factors included to explore which aspects of the morphology are more meaningful for explaining the correlation between morphological diversity and species richness.
Data analysis Two curves were extracted from pronotum and elytron contours to represent their external forms. Each curve was resampled into 25 or 50 semi-landmarks (SLM) by length, respectively (Fig. 1). The curves and semi-landmarks were digitized with TPS-DIG 2.05 [51]. Data les used for morphological analysis were completed by including semi-landmarks and landmarks [52] into the same data les [23,53,54]. Principal component analysis (PCA) and geometric modelling of mathematical spaces formed by the PC axes were used to ordinate shape variation of the whole dataset [55]. Patterns of test shape (pronotum and elytron) variation represented by the rst three principal components were included in this study. Each shape

Declarations
Ethics approval Not applicable. Consent for publication model shown in the shape-model matrix was calculated at an equally spaced interval along each PC axis (Fig. 2, 3). It should be remembered that the apparent deformations indicated different shape variations in the two-dimensional plot constructed by different principal components. This pattern was formed by the rst three principal components with the largest proportion was shown (Fig. 2, 3), deformation of test characters changed to some extent in the plots which were structured by other principal components (Suppl. 2). For the correlation in feeding habits test, test samples of 26 subfamilies and 2 families which had no subfamilies were divided into omnivorous, phytophagous, coprophagous and non-feeding categories based on their food habits [18]. The Scarabarinae were described as coprophagous as a result of their dung-type soft food, some groups of scarabs (Amphicominae, Glaphyrinae, Asealinae, Lampriminae, Lucaninae, Syndesinae, Passalinae, Aulacocyclinae, and the subfamilies of Scarabaeidae: Cetoniinae, Dynastinae, Melolonthinae, Orphninae, Rutelinae) were treated as phytophagous as a result of their planteating habits [32,[56][57][58][59][60][61][62][63]. Some groups with complex and diverse food habits (e.g., mycophagy, coprophagy) were assigned to the omnivourous groups as omnivorous (Aphodiinae, Anaidinae, Ceratocanthinae, Hybosorinae, Bolboceratinae, Geotrupinae, Ochodaeinae, Troginae, Omorginae, Liparochrinae and Chaetocanthinae) [64,65]. In addition, Pleocominae and Diphyllostomatidae were assigned to the non-feeding groups because of their vestigial mouthparts and unclear ecology [26,32].
For the correlation test based on phylogenetic factor, the whole dataset was separated into 11 families and 26 subfamilies (Suppl. 1). Correlations between pronotum and elytron morphological diversity and species/genera richness values in family category and subfamily category were analyzed via regression analysis ( Fig. 1) [66,67]. Also, in the subfamily category test, Diphyllostomatidae and Glaresidae were treated as a subfamily-group for analyzing because that these two families had no subfamilies. Exploration of correlations between morphological diversity and species richness was achieved through regression analysis [66]. Morphological diversity was quanti ed as the total Procrustes variance [68] in MORPHO J 1.06a [69]; species richness value was obtained by counting the number of genera and species in different test groups (Fig. 1). A best-t line estimated from morphological diversity and species/genus richness values was used to test the linear relationship of these parameters. The correlation between morphological diversity and species richness was measured using Spearman correlation coe cient [54,66,68], where necessary the species richness was transformed to the Log 10 scale (Fig. 4) [70]. In addition, the worldwide species/genus richness was counted [20,[71][72][73], and the regression analysis used to document the correlation between the worldwide species/genus richness and the test groups' morphological diversity based on test species/genus richness and morphological diversity via covariance analysis [74] to verify our results (Fig. 4, also see Suppl. 3).

Figure 1
Species/genus richness and morphological diversity among groups. Four parameters were shown based on phylogenetic factor and feeding types. The curves on the pronotum and elytron of the beetles (Callistethus plagiicollis plagiicollis) demonstrated the morphological representation used in geometric morphometric analysis. Curve 1 (orange) was resampled in 25 semi-landmarks (SLM); curve 2 (green) was resampled in 50 semi-landmarks (SLM). Samples grouped by phylogenetic factor. c, d: Samples grouped by ecological factor. Each shape model, shown in the lower matrix was calculated at equally spaced intervals along PC axes with coordinates for individual shape model calculations indicated by marks and numbers in the two-dimensional PC shape spaces. Overlaying shapes got from the deformations along each axis were shown at the end of the matrix.

Figure 3
Morphological variations of elytron in Scarabaeoidea. 93.362% of observed shape variation was represented by rst three PCs. a, c: Variations of PC1 and PC2. b, d: Variations of PC 2 and PC3. a, b: Samples grouped by phylogenetic factor. c, d: Samples grouped by ecological factor. Each shape model, Page 22/24 shown in the lower matrix was calculated at equally spaced intervals along PC axes with coordinates for individual shape model calculations indicated by marks and numbers in the two-dimensional PC shape spaces. Overlaying shapes got from the deformations along each axis were shown at the end of the matrix.

Figure 4
Relationship between morphological diversity and species richness. a-d: Samples grouped by phylogenetic factor. e, f: Samples grouped by ecological factor. Correlation between pronotum diversity and species/genus richness in family category (a), subfamily category (b) and feeding types (e).
Correlation between elytron diversity and species/genus richness in family category (c), subfamily category (d) and feeding types (f). Values of species/genus richness were transformed to their log10 values. TGR= tested genus richness, TSR= tested species richness, WGR= worldwide genus richness, WSR= worldwide species richness.
Supplementary Files