Abstract

Pangolins are among the most endangered groups of mammals, comprising eight extant species delineated into three genera. Despite several studies dedicated to their skeletal anatomy, the potential taxonomic insight from cranial morphological variation in extant Pholidota is yet to be assessed with modern geometric morphometric methods. We present the first comprehensive study on the cranial morphology of extant pangolins and discuss its implications for the taxonomy and evolution of the group. We performed landmark-based morphometric analyses on 241 museum specimens to describe the variation in skull shape in seven of the eight extant species. Our analyses revealed genus- and species-level morphological discrimination, with Asian species (Manis spp.) being grouped together, whereas African pangolins present distinct skull shapes between small (Phataginus spp.) and large (Smutsia spp.) species. Analyses of allometry also identified a set of traits whose allometric trajectories distinguish Asian from African specimens. Finally, we uncovered intraspecific variation in skull shape in white-bellied pangolins (Phataginus tricuspis) that partly corroborates recent DNA-based differentiation among biogeographically distinct populations. Overall, our results shed light on the morphological diversity of the skull of these enigmatic myrmecophagous mammals and confirm the genus-level classification and cryptic diversity within the white-bellied pangolin revealed by molecular phylogenetics.

INTRODUCTION

Extant pangolins [Pholidota (Weber, 1904) Manidae (Gray, 1821)] are currently the most heavily poached mammals in the world (Challender et al., 2014; Zhou et al., 2014). The increasing demand for pangolin scales in Chinese traditional medicine is driving this entire group to the brink of extinction. In that context, detailed morphological and genetic studies constitute a prerequisite to trace the geographical origin of seized specimens and may prove fruitful to delimitate new and more effective conservation management units (von Helversen et al., 2001; Hebert et al., 2004; Bickford et al., 2007; Moraes-Barros et al., 2007; Galatius et al., 2012; Hautier et al., 2014; Sveegaard et al., 2015; Gaubert et al., 2016). The eight extant pangolin species are restricted to tropical and intertropical regions of sub-Saharan Africa and Southeast Asia (Gaubert, 2011; Gaudin et al., 2019). They are classified in a single family (Manidae) within the order Pholidota, which is most closely related to carnivores within placental mammals (Murphy et al., 2001a, b). Pangolins have evolved a set of highly distinctive morphological traits, such as toothless jaws and an elongated rostrum, linked to their highly specialized diet of ants and termites (Ferreira-Cardoso et al., 2019; Gaudin et al., 2019). Despite this unique evolutionary history and their current protection status, pangolins are among the least studied placental mammals, with some aspects of their phenotypic variation (i.e. cranial shape) still being completely unexplored (Gaubert et al., 2018).

The taxonomy of extant pangolins has been relatively unstable, varying from a single genus (Manis Linnaeus, 1758; Jentink, 1882; Emry, 1970) to six genera (Pocock, 1924). Based on both morphological and molecular phylogenetic studies, the four African species have been split into two genera (Phataginus Rafinesque, 1821 and Smutsia Gray, 1865), whereas the four Asian species remained in a single genus (Manis) (Gaudin et al., 2009; Gaubert et al., 2018). A recent molecular phylogeographical study also identified six cryptic lineages within the widely distributed white-belied pangolin (Phataginus tricuspis Rafinesque, 1821), showing an unexpected intraspecific molecular divergence (Gaubert et al., 2016). Likewise, several genetic clusters have been identified in the range of the Sunda pangolin (Manis javanica Desmarest, 1822), although their geographical delineation remains unclear (Zhang et al., 2015; Nash et al., 2018). Many aspects of inter- and intraspecific diversity of pangolins are still unexplored, and this newly described cryptic diversity has great potential to influence regional conservation policies and to identify the geographical origin of illegally trafficked specimens.

Geometric morphometrics has proved to be an efficient way to define the patterns of shape variation associated with species delimitation (interspecific taxonomy; e.g. Cardini & O’Higgins, 2004; Villemant et al., 2007) and to unveil hidden morphological variation (intraspecific taxonomy; e.g. Hautier et al., 2014, 2017; Miranda et al., 2018). Such methods also enable us to understand the role of interactions between size and shape, with the variation in shape correlated with size (allometry) being one of the main factors contributing to the integrated evolution of cranial shape (Cardini & Polly, 2013; Klingenberg, 2013; Cardini, 2019). In mammals, ontogenetic and evolutionary allometry is mainly associated with elongation of the rostrum (Cardini & Polly, 2013; Cardini et al., 2015). However, the variation in shape and size of the pangolin skull has never been described in detail or quantified formally, although the elongated, toothless snout constitutes one of the main characteristics of their skull. A precise characterization of the changes in shape associated with growth in pangolins should enable us to understand the extent to which allometry has contributed to the evolution of the skull.

Here, we applied geometric morphometric methods to study the variation in shape of the skull within and among the eight extant pangolin species, with the aim of assessing their current taxonomy and the recently identified molecular lineages. First, we examined the patterns of ontogenetic allometry in extant Pholidota using regression and trajectory analyses. Second, we explored the variation in cranial morphological among and within the three extant genera by performing regression, ordination and discriminant analyses. Finally, we investigated the variation in the shape of the skull in two wide-ranging species (P. tricuspis and M. javanica) in order to assess whether molecular-based cryptic lineages (Gaubert et al., 2018; Nash et al., 2018) differ in skull shape.

MATERIAL AND METHODS

Biological sampling and data collection

The material used in this study belongs to the collections of the Natural History Museum (BMNH) in London (UK), the Museum für Naturkunde (MfN) in Berlin (Germany), the Muséum National d’Histoire Naturelle (MNHN) in Paris (France), the Royal Museum for Central Africa (KMMA/RMAC) in Tervuren (Belgium), the American Museum of Natural History (AMNH) in New York (NY, USA) and the National Museum of Natural History (USNM) in Washington (DC, USA).

Our dataset is the result of the landmarking of 243 specimens from the eight extant species of pangolins (Fig. 1): P. tricuspis (Fig. 2; white-bellied pangolin, N = 97), Phataginus tetradactyla (Linnaeus, 1766) (black-bellied pangolin, N = 23), Smutsia gigantea (Illinger, 1815) (giant ground pangolin, N = 16), Smutsia temminckii (Smuts, 1832) (Cape pangolin, N = 21), Manis crassicaudata Geoffroy, 1803 (Indian pangolin, N = 10), M. javanica (Sunda pangolin, N = 38), Manis pentadactyla Linnaeus, 1758 (Chinese pangolin, N = 36) and Manis culionensis (de Elera, 1915) (Palawan pangolin, N = 2; Supporting Information, Appendix S1). The skull shape of the Palawan pangolin could not be assessed quantitatively owing to the low number of specimens available and the relatively large number of missing landmarks [a principal components analysis (PCA) including this species is provided in the Supporting Information, Fig. S1]. Taxa were identified based on a list of morphological criteria identified in the present study and in previously published works (Supporting Information, Appendix S2; Hatt et al., 1934; Gaudin et al., 2009, 2019). Different sizes of skulls were included to account for the change in shape during ontogeny. Given the absence of teeth in pangolins, the determination of age is not straightforward. For each species, juveniles were defined arbitrarily as those for which size (estimated by the centroid size; see below) fell below the first quartile value (25% smallest specimens) for each species.

Figure 1.

Geographical distribution of the sampled specimens (N = 243) belonging to eight species of pangolins. Mean shapes are illustrated for each species, except for Manis culionensis. The phylogenetic relationships used the branch lengths from Gaubert et al. (2018). 1, Phatagininae; 2, Smutsiinae; 3, Maninae.

Figure 1.

Geographical distribution of the sampled specimens (N = 243) belonging to eight species of pangolins. Mean shapes are illustrated for each species, except for Manis culionensis. The phylogenetic relationships used the branch lengths from Gaubert et al. (2018). 1, Phatagininae; 2, Smutsiinae; 3, Maninae.

Figure 2.

Map of Africa, with the locations of the sampled Phataginus tricuspis (N = 96) coloured according to the cryptic lineage (Gaubert et al., 2016) to which they were attributed (right). The tree topology of the intraspecific affinities is based on Gaubert et al. (2018).

Figure 2.

Map of Africa, with the locations of the sampled Phataginus tricuspis (N = 96) coloured according to the cryptic lineage (Gaubert et al., 2016) to which they were attributed (right). The tree topology of the intraspecific affinities is based on Gaubert et al. (2018).

We placed 75 three-dimensional landmarks on pangolin skulls using a Revware MicroScribe M 3D digitizer (Fig. 3; Supporting Information, Table S1). Our selection of landmarks was based on previous works focused on mammalian taxa (e.g. Goswami, 2006; Hautier et al., 2017). In a significant number of specimens, the premaxillae were absent, loosely attached or broken and could not be landmarked. In pangolins, the jugal bone is often absent. However, some M. pentadactyla specimens presented a complete zygomatic arch (see also Emry, 2004). In such cases, the landmarks 58/72 (zygomatic process of the maxillae) were hard to define accurately. They were therefore considered as missing and estimated a posteriori.

Figure 3.

Landmarks digitized on the skull of Phataginus tricuspis (BMNH 34.6.2.92) in lateral (A), ventral (B) and dorsal (C) views. Red and blue numbers represent landmarks placed ventrally and dorsally, respectively.

Figure 3.

Landmarks digitized on the skull of Phataginus tricuspis (BMNH 34.6.2.92) in lateral (A), ventral (B) and dorsal (C) views. Red and blue numbers represent landmarks placed ventrally and dorsally, respectively.

Thin plate spline interpolations (Gunz et al., 2009) were computed to estimate missing landmarks for each group (e.g. species, cryptic lineages). This approach was implemented in the software package ‘geomorph’ v.3.5.0 (Adams et al., 2017) in R (R Development Core Team, 2013). A generalized Procrustes analysis (Rohlf & Slice, 1999) was performed on all sets of landmarks. All specimens were scaled to centroid size, optimally translated and rotated using a least-squares criterion. The coordinates retrieved by the generalized Procrustes analysis represented the shape of the skull of each specimen. An ANOVA was performed on a subset of our data for which sex determination was available (N = 120), in order to assess sexual dimorphism in skull shape.

Allometry in extant pangolins

The study of allometry can focus on three different levels (Cheverud, 1982; Klingenberg, 2016). Morphological changes can be associated with phylogenetic differences in size (evolutionary allometry), variation in size within a single ontogenetic stage (static allometry) or variation in size attributable to individual growth within a single species (ontogenetic allometry; e.g. Foth et al., 2016; Esquerré et al., 2017; Gray et al., 2019). Here, we quantified evolutionary allometry and considered three traits of ontogenetic allometric trajectories: direction (slope), magnitude (length) and intercept. These aspects were controlled/investigated with three different analyses.

First, a phylogenetic multivariate regression of the mean shape of adult specimens (N = 173) against log-centroid size was performed to assess evolutionary allometry. This analysis was performed with procD.pgls from ‘geomorph’ using a consensus phylogeny from Gaubert et al., (2018). Then, we focused on ontogenetic allometry using interspecific (N = 241) and intraspecific datasets. The cryptic lineages dataset including P. tricuspis included only 96 specimens (Fig. 2; KMMA 30808 was discarded owing to imprecise information on location). Details about the analysis of the M. javanica dataset can be found in the Supporting Information (Appendix S3). First, a multivariate regression was performed to assess covariation patterns between the logarithm of the centroid size and Procrustes-aligned coordinates (raw shapes) using procD.lm from ‘geomorph’. The hypothesis of parallel group slopes was assessed with a homogeneity of slopes (HOS) test. This test consists of an analysis of covariance (ANCOVA) to assess the influence of size, groups and the interaction of size and groups on skull shape. The HOS test includes pairwise comparisons between groups (species/cryptic lineages) to assess significant differences of both the direction (angles) and magnitude (amount of change in shape with size) of allometric trajectories. Significance was assessed with a residual randomization permutation procedure with 10 000 iterations. The HOS tests were performed with advanced.procD.lm from the package ‘geomorph’. The HOS tests were complemented by graphical representations of allometric trajectories and skull shape deformations (Adams & Nistri, 2010; Esquerré et al., 2017). We plotted the first principal component (PC1) of the predicted values of multivariate regression of shape on log-centroid size vs. log-centroid size regressions for each species. We then assessed the significance of differences between the intercepts using a Tukey means comparison, to test for changes in shape explained by species differences in the resulting morphospace. Thin plate splines were generated using the function tps3D from the R package ‘Morpho’ v.2.5.1 (Schlager, 2017) in order to characterize differences in shape between the smallest and the largest specimens for each species (ontogenetic allometry). We then visualized landmark displacement during ontogeny using the function deformGrid3d from the same package. A phenotypic trajectory analysis (PTA) was also performed as a complementary analysis to the HOS test (see Supporting Information, Appendix S3; Adams & Collyer, 2009; Collyer & Adams, 2013).

Second, when the HOS test revealed parallel trajectories, we performed an additional analysis (overlap test) on multivariate shape data to test for overlap in ontogenetic allometric trajectories by comparing their differences with a set of 10 000 permutations (Piras et al., 2011; Esquerré et al., 2017). Intercepts were tested at x = 5.125, because we lacked fetuses and neonates of small size (close to x = 0), which could result in an incorrect estimate of minimum size (x = 0). The overlap test was performed only on the interspecific dataset.

Finally, if slopes were overlapping between species, a third analysis was performed to identify peramorphosis/paedomorphosis. This ‘heterochrony test’ enables the characterization of differences in skull shape at maximum size (Piras et al., 2011; Esquerré et al., 2017). Significance was assessed by comparing these differences with a set of 10 000 permutations. This test was performed only on the interspecific dataset.

Visualization of shape variation and statistical analysis

The variation of skull shapes was visualized using a PCA (Dryden & Mardia, 1993). Data were analysed without allometric correction in order to allow for a comparison of the results with morphological discrete characters used in phylogenetic analyses. Analyses with the allometry-corrected shapes for the interspecific and P. tricuspis datasets can be found in the Supporting Information (Appendix S3). Allometry-corrected shapes were obtained as the residuals of pooled within-group regressions. This method is used to obtain a common estimate for the allometry when comparing several groups (Sidlauskas et al., 2011; Benítez et al., 2013; Klingenberg, 2016). Pooled within-group regressions were performed in MorphoJ v1.06d (Klingenberg, 2011). For simplicity, when PCAs were performed on allometry-corrected shapes, axes were designated as PCres.

Interspecific variation in extant pangolins

A PCA was performed excluding specimens considered as juveniles (N = 173). A second PCA was performed on a dataset including juveniles and two specimens of M. culionensis (N = 243). We used a .ply surface of a micro-computed tomography-scanned P. tricuspis skull (BMNH 34.6.2.92), the species most closely resembling the mean shape. The surface was then deformed to the mean skull shape of Pholidota and used to visualize the variation in skull shape along the first three principal components (PCs). Triangular mesh warping via thin plate spline was performed with the package ‘Morpho’.

A multivariate ANOVA was performed to assess whether skull shape differed between taxonomic groups. Pairwise comparisons between least-squares means were performed using advanced.procD.lm. If taxonomy had a significant effect on skull shape, a leave-one-out cross-validated linear discriminant analysis (LDA) was performed on a set of PCs explaining 90% of variance. The leave-one-out procedure allows evaluation of the accuracy with which unknown specimens can be identified (e.g. Evin et al., 2013). Linear discriminant analyses and post hoc classification methods were performed with the ‘MASS’ package (Venables & Ripley, 2002) in R.

Intraspecific variation in P. tricuspis and M. javanica

A second subset including only adult P. tricuspis specimens was analysed separately in order to describe the intraspecific variation in skull shape (N = 71). We used the cryptic genetic lineages defined by Gaubert et al. (2016) to test the variation in skull shape linked to geographical distribution within P. tricuspis. Specimens were sorted according to six regions (Fig. 2): Western Africa (WAF), Ghana (GHA), Dahomey Gap (DHG), Western Central Africa (WCA), Central Africa (CAF) and Gabon (GAB). Skull shape difference tests and cross-validated LDA described in the previous paragraph were repeated on this dataset. If a posteriori attribution errors were consistently detected between two regions, these were merged, and the protocol was repeated with the new specimen sorting in order to test for its potential in assigning specimens with unknown geographical origin.

An equivalent protocol was applied to assess intraspecific variation in the skull of M. javanica. The geographical delimitation of cryptic lineages within this species is still uncertain (Zhang et al., 2015; Nash et al., 2018). These preliminary analyses can be found in the Supporting Information (Appendix S3).

The original landmark coordinates used in this study are provided in the Supporting Information (Appendix S1).

RESULTS

Allometry in extant pangolins

A multivariate regression (Table 1) performed on raw shape variables revealed significant effects of log-transformed centroid size (F1,227 = 103.59, P < 0.001, R2 = 0.16), species grouping (F6,227 = 49.63, P < 0.001) and of an interaction between centroid size and species grouping (F6,227 = 1.35, P < 0.001) on shape. When we accounted for phylogeny, the effect of size was marginally non-significant (P = 0.09), with evolutionary allometry explaining roughly one-third of skull shape variance (R2 = 0.27; Supporting Information, Table S2). The non-significant P-value was probably attributable to the low number of species (N = 7). The HOS test pairwise comparisons of ontogenetic allometric trajectories revealed small, significant differences (low z-values) between the slopes of S. gigantea compared with both M. javanica and M. pentadactyla (Supporting Information, Table S3). The remaining species did not present significantly different slope angles, which implied that allometric trajectories were parallel within each genus. The results for the PTA are presented in the Supporting Information (Appendix S3; Fig. S2; Tables S4 and S5).

Table 1.

ANOVA of shape (Procrustes coordinates) ~ log(centroid size)*species

 d.f. R2 F P-value 
Log(centroid size) 0.16 103.59 < 0.001* 
Species 0.47 49.63 < 0.001* 
Log(centroid size):species 0.01 1.35 < 0.001* 
Residuals 227 0.37 – – 
Total 240 – – – 
 d.f. R2 F P-value 
Log(centroid size) 0.16 103.59 < 0.001* 
Species 0.47 49.63 < 0.001* 
Log(centroid size):species 0.01 1.35 < 0.001* 
Residuals 227 0.37 – – 
Total 240 – – – 

The randomized residual permutation procedure used 10 000 permutations.

*Significant P-value.

Table 1.

ANOVA of shape (Procrustes coordinates) ~ log(centroid size)*species

 d.f. R2 F P-value 
Log(centroid size) 0.16 103.59 < 0.001* 
Species 0.47 49.63 < 0.001* 
Log(centroid size):species 0.01 1.35 < 0.001* 
Residuals 227 0.37 – – 
Total 240 – – – 
 d.f. R2 F P-value 
Log(centroid size) 0.16 103.59 < 0.001* 
Species 0.47 49.63 < 0.001* 
Log(centroid size):species 0.01 1.35 < 0.001* 
Residuals 227 0.37 – – 
Total 240 – – – 

The randomized residual permutation procedure used 10 000 permutations.

*Significant P-value.

The ontogenetic allometric trajectories overlapped in most species presenting parallel slopes except for P. tricuspis (Supporting Information, Table S6). The trajectory for P. tetradactyla overlapped with all the others except those of M. javanica and P. tricuspis. When comparison of the intercepts was performed considering shape predictions for x = 0, the ontogenetic allometric trajectories overlapped in all species (Supporting Information, Table S7). The heterochrony test showed that all species with overlapping trajectories presented heterochronic shifts with respect to each other (Supporting Information, Table S8).

The intercepts of the regressions of predicted values of shape against size (Fig. 4) were relatively distinct within Phataginus and Manis, whereas both Smutsia species presented overlapping trajectories, as revealed by the Tukey comparisons (Supporting Information, Table S9). Considering the predicted shapes for minimum and maximum size resulting from the multivariate regression, the main size-related intraspecific morphological change was the increase in length of the rostrum (Figs 4, 5). Landmarks in the anterior part of the nasal and maxilla tended to be more anterodorsally positioned, and the nasal projected more posteriorly (e.g. Fig. 4, 5). The braincase was relatively lower in the adult, across all species, with the dorsal landmarks on the midline of the skull being more ventral when compared with their position in the juveniles (Fig. 5). Additionally, the landmarks placed on the zygomatic process of the maxilla and associated structures (Fig. 3; landmarks 10, 57 and 58) showed a tendency to project more posteriorly in larger specimens (Fig. 5A–D). In contrast, S. temminckii showed no allometric growth of the posterior projection of the zygomatic processes (Fig. 5E). In contrast, this species presented the most significant change in the anterior projection of the zygomatic process of the squamosal, with this structure being noticeably less projected in smaller specimens (Fig. 5E).

Figure 4.

Allometric trajectories among seven pangolin species (N = 241). A, the x-axis values are the log-transformed centroid sizes for each specimen; the y-axis values are the principal component 1 of the predicted values of a multivariate regression of shape on size. B–H, deformed meshes for the maximum (top) and minimum (bottom) shapes predicted from a multivariate Procrustes regression for each species are presented.

Figure 4.

Allometric trajectories among seven pangolin species (N = 241). A, the x-axis values are the log-transformed centroid sizes for each specimen; the y-axis values are the principal component 1 of the predicted values of a multivariate regression of shape on size. B–H, deformed meshes for the maximum (top) and minimum (bottom) shapes predicted from a multivariate Procrustes regression for each species are presented.

Figure 5.

Mapping of ontogenetic variation in shape in seven pangolin species. A, Manis crassicaudata. B, Manis pentadactyla. C, Manis javanica. D, Smutsia gigantea. E, Smutsia temminckii. F, Phataginus tricuspis. G, Phataginus tetradactyla. Grey dots represent relative Procrustes coordinates positions in the juvenile (minimum-sized specimen), and red dots represent coordinates for the adult (maximum-sized specimen). Deformed meshes corresponding to juveniles were superimposed with the three-dimensional Procrustes coordinates.

Figure 5.

Mapping of ontogenetic variation in shape in seven pangolin species. A, Manis crassicaudata. B, Manis pentadactyla. C, Manis javanica. D, Smutsia gigantea. E, Smutsia temminckii. F, Phataginus tricuspis. G, Phataginus tetradactyla. Grey dots represent relative Procrustes coordinates positions in the juvenile (minimum-sized specimen), and red dots represent coordinates for the adult (maximum-sized specimen). Deformed meshes corresponding to juveniles were superimposed with the three-dimensional Procrustes coordinates.

Interspecific variation of the skull shape in extant pangolins

A Procrustes ANOVA revealed that both sex (F1,105 = 1.80, P = 0.085; F1,106 = 2.03, P = 0.057) and the interaction between sex and species grouping (F6,105 = 0.89, P = 0.656; F6,106 = 1.07, P = 0.332) did not significantly influence shape, with and without considering size as a covariate in the model, respectively (N = 120; Supporting Information, Table S10). The Procrustes ANOVA performed with the adult dataset revealed a significant effect of species on skull shape (F6,166 = 45.02, P < 0.001; Table 2). Pairwise comparisons showed that all species presented significantly different skull shapes (Supporting Information, Table S11). The variation in skull shape was visualized using a PCA performed on the raw shape variables of the seven pangolin species (Fig. 6). The results of a PCA including juveniles and M. culionensis (Supporting Information, Fig. S1) and a full analysis of allometry-corrected skull shape are presented in the Supporting Information (Appendix S3; Fig. S3; Tables S12, S13). Despite the non-parallel slopes, the variance explained by the interaction between size and species was relatively low (Table 1), allowing us to use the residuals of a multivariate regression of shape on size as allometry-corrected shapes.

Table 2.

ANOVA of shape ~ taxa/geographical groups of adult specimens of the interspecific and Phataginus tricuspis datasets

Datasets N taxa/geo N d.f. R2 F P-value 
Interspecific 173 166 0.62 45.02 < 0.001* 
P. tricuspis (Gaubert et al., 201670 65 0.16 3.12 < 0.001* 
P. tricuspis (geographical groups) 70 67 0.11 4.23 < 0.001* 
Datasets N taxa/geo N d.f. R2 F P-value 
Interspecific 173 166 0.62 45.02 < 0.001* 
P. tricuspis (Gaubert et al., 201670 65 0.16 3.12 < 0.001* 
P. tricuspis (geographical groups) 70 67 0.11 4.23 < 0.001* 

The randomized residual permutation procedure used 10 000 permutations. Significant P-values indicate differences between skull shapes of taxa/geographical groups.

Abbreviation: N taxa/geo, number of taxa or geographical groups used as factors.

*Significant P-value.

Table 2.

ANOVA of shape ~ taxa/geographical groups of adult specimens of the interspecific and Phataginus tricuspis datasets

Datasets N taxa/geo N d.f. R2 F P-value 
Interspecific 173 166 0.62 45.02 < 0.001* 
P. tricuspis (Gaubert et al., 201670 65 0.16 3.12 < 0.001* 
P. tricuspis (geographical groups) 70 67 0.11 4.23 < 0.001* 
Datasets N taxa/geo N d.f. R2 F P-value 
Interspecific 173 166 0.62 45.02 < 0.001* 
P. tricuspis (Gaubert et al., 201670 65 0.16 3.12 < 0.001* 
P. tricuspis (geographical groups) 70 67 0.11 4.23 < 0.001* 

The randomized residual permutation procedure used 10 000 permutations. Significant P-values indicate differences between skull shapes of taxa/geographical groups.

Abbreviation: N taxa/geo, number of taxa or geographical groups used as factors.

*Significant P-value.

Figure 6.

Principal components (A, PC1 vs. PC2; B, PC1 vs. PC3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3) with associated variation in shape for crania of seven pangolin species (N = 173).

Figure 6.

Principal components (A, PC1 vs. PC2; B, PC1 vs. PC3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3) with associated variation in shape for crania of seven pangolin species (N = 173).

Linear regressions performed on the first 30 PCs (90%) showed that size was significantly correlated with PCs 1–4 and PC7 (Supporting Information, Table S14). Size-related morphological changes captured by PC1 appeared to be associated mainly with evolutionary allometry, whereas change in shape along PC2 recovered differences related to ontogenetic allometry. The first two PCs explained 49.59% of the total variance (33.95 and 15.64%, respectively).

Principal component 1 was positively correlated with a high and wide rostrum, a nasofrontal inflation associated with an orbital constriction, posteriorly projected zygomatic processes of the maxillary and a high and wide braincase, with dorsal squamosal-parietal-frontal junctions. African specimens tended to display mostly negative PC1 scores, whereas the Asian clade exhibited positive PC1 scores (Fig. 6A, B). The sole exception was the African S. gigantea, which presented positive PC1 values and grouped with Asian specimens (Fig. 6). Principal component 1 also separated the two African genera (Phataginus presented the most negative scores). Juvenile specimens of Asian species were characterized by less positive PC1 values, plotting closer to African pangolins (see ‘Results: allometry within extant pangolins’; Supporting Information, Fig. S1).

Principal component 2 was positively correlated with a long rostrum, a long posterior projection of the premaxilla on the midline, an anterolaterally projecting zygomatic process of the squamosal and a braincase with a pseudorectangular shape in dorsal view (Fig. 6). Principal component 2 separated the three species of Manis and the two Smutsia. Smutsia temminckii and M. pentadactyla scored the lowest PC2 average values, whereas M. javanica scored the highest. Phataginus spp., S. gigantea and M. crassicaudata presented PC2 scores ranging in between the two groups.

Principal component 3 scores were positively correlated with the anterior projection of the anterior flanges of the frontal, a wide and long palatine, a long infraorbital canal, an anteroposteriorly elongated dorsal edge of the zygomatic processes of the squamosal and posteriorly projected pterygoid hamuli (Fig. 6B). This axis separated species within Phataginus and Smutsia. Phataginus tetradactyla scored extremely positive PC3 values, whereas S. gigantea scored the most negative PC3. The three Manis species and S. temminckii presented values slightly above zero, on average, and P. tricuspis showed mostly negative values, but not as negative as S. gigantea (Fig. 6B).

An LDA was performed to take 90% of the variance into account (first 30 PCs). Linear discriminant analysis group a posteriori probabilities retrieved 100% accuracy for species attribution (Supplemental Information, Appendix S4). Specimens were grouped by species and were well discriminated on linear discriminants 1 and 2 (LD1, 54.0%; LD2, 15.6%; Fig. 6C, D). Asian pangolin skulls presented negative LD1 values, whereas African pangolins scored both negative (Smutsia) and extremely positive values (Phataginus; Fig. 6C, D). The three species of Maninae Gray, 1821 resembled each other the most, being discriminated only by LD2. Within the African clade, LD2 discriminated the two Phataginus species (P. tetradactyla showed the highest LD2 values). Smutsia temminckii and S. gigantea were well discriminated by LD1. Smutsia spp. showed the highest LD3 values.

The analyses on the allometry-corrected shapes revealed some differences that are discussed in detail in the Supporting information (Fig. S3; Tables S12 and S13).

Intraspecific variation of the skull shape in extant pangolins

Intraspecific variation in P. tricuspis

A multivariate regression revealed that log-transformed centroid size (F1,84 = 11.56, P < 0.001) and geographical distribution (cryptic lineages; F4,84 = 3.71, P < 0.001) had a highly significant effect on the cranial shape (Table 3). It also retrieved a significant effect of the interaction between size and geographical distribution (F4,84 = 1.11, P = 0.001). However, the pairwise matrix effect sizes were relatively small, and corresponding P-values were not significant. Therefore, all cryptic lineages presented parallel allometric trajectories (Supporting Information, Fig. S4; Table S15). A Procrustes ANOVA revealed that cryptic lineages presented different skull shapes (F4,65 = 3,12, P < 0.001; Table 2). Pairwise comparisons showed that all tested cryptic lineages presented significantly different skull shapes except for Ghana and Western Africa (Supporting Information, Table S16) Gabon was not tested owing to lack of replicates.

Table 3.

ANOVA of shape (Procrustes coordinates) ~ log(centroid size)*Phataginus tricuspis cryptic lineages (N = 96)

 d.f. R2 F P-value 
Log(centroid size) 0.10 11.56 < 0.001* 
Cryptic lineages 0.15 3.71 < 0.001* 
Log(centroid size):cryptic lineages 0.04 1.11 0.001* 
Residuals 85 0.71 – – 
Total 95 – – – 
 d.f. R2 F P-value 
Log(centroid size) 0.10 11.56 < 0.001* 
Cryptic lineages 0.15 3.71 < 0.001* 
Log(centroid size):cryptic lineages 0.04 1.11 0.001* 
Residuals 85 0.71 – – 
Total 95 – – – 

The randomized residual permutation procedure used 10 000 permutations.

*Significant P-value.

Table 3.

ANOVA of shape (Procrustes coordinates) ~ log(centroid size)*Phataginus tricuspis cryptic lineages (N = 96)

 d.f. R2 F P-value 
Log(centroid size) 0.10 11.56 < 0.001* 
Cryptic lineages 0.15 3.71 < 0.001* 
Log(centroid size):cryptic lineages 0.04 1.11 0.001* 
Residuals 85 0.71 – – 
Total 95 – – – 
 d.f. R2 F P-value 
Log(centroid size) 0.10 11.56 < 0.001* 
Cryptic lineages 0.15 3.71 < 0.001* 
Log(centroid size):cryptic lineages 0.04 1.11 0.001* 
Residuals 85 0.71 – – 
Total 95 – – – 

The randomized residual permutation procedure used 10 000 permutations.

*Significant P-value.

Principal component 1 explained 12.8% of the variance of cranial shape within P. tricuspis and was positively correlated with a larger skull height and width, an anteroposteriorly short orbit, large tympanic bullae and a relatively round occipital region. Specimens from Central Africa mostly scored negative PC1 values (Fig. 7A, B). On average, specimens from WAF and GHA presented the most positive PC1 scores. The Western Central African cluster also presented mostly positive PC1 scores, whereas specimens from the Dahomey Gap presented a wide range of PC1 scores, varying from negative to positive values. The only specimen from Gabon scored negative PC1 values, plotting near the CAF morphospace.

Figure 7.

Principal components (A, PC1 vs. PC2; B, PC1 vs. PC3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3) with associated variation in shape for crania of six cryptic lineages of Phataginus tricuspis (N = 71; Gaubert et al., 2016).

Figure 7.

Principal components (A, PC1 vs. PC2; B, PC1 vs. PC3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3) with associated variation in shape for crania of six cryptic lineages of Phataginus tricuspis (N = 71; Gaubert et al., 2016).

Principal component 2 explained 9.7% of the variance and was positively correlated with a shorter palate with short maxillary projections, anteriorly projecting squamosal roots and shorter tympanic bullae well separated from the postglenoid foramina. Although specimens from CAF had a wide distribution along PC2, WCA, WAF and GHA presented a much narrower range of PC2 scores in the middle of the distribution. On average, DHG presented the most negative PC2 values.

Principal component 3 explained 6.4% of the variance (Fig. 7B) and did not segregate specimens according to geographical origin.

The LDA performed on the first 34 PCs (90% variance) discriminated WAF–GHA, DHG and CAF–WCA cryptic lineages, along LD1 (59.1%) (Fig. 7C). Western Africa–Ghana presented the most positive LD1 values, whereas CAF–WCA specimens presented mostly negative LD1 values (Fig. 7C, D). Specimens from DHG presented intermediate positive LD1 values. LD2 (23.7%) discriminated DHG skulls (negative values) from all other cryptic lineage specimens (Fig. 7C). LD3 (13.0%) discriminated the WCA specimens (most positive values) from the remaining lineages (Fig. 7D). Group a posteriori probabilities retrieved a 75.4% attribution accuracy (see Supporting Information, Appendix S4). The vast majority of incorrect attributions were found in the major divisions WAF–GHA and CAF–WCA. Based on this result, we performed an additional LDA with a priori attributions of WAF–GHA specimens to a western group (WES) and CAF–WCA specimens to an eastern group (CEN), while keeping DHG as a separate group (Fig. 8). Group a posteriori probabilities of the LDA of the three groups shows an attribution accuracy of 95.7% (see Supporting Information, Appendix S4).

Figure 8.

Linear discriminant analysis of the Phataginus tricuspis sample (N = 70) divided into three management units. Consensus shapes (mean shapes) of the proposed management units in lateral view. A, Central African region (CEN). B, Dahomey Gap region (DHG). C, western African region. Black dots are landmark positions.

Figure 8.

Linear discriminant analysis of the Phataginus tricuspis sample (N = 70) divided into three management units. Consensus shapes (mean shapes) of the proposed management units in lateral view. A, Central African region (CEN). B, Dahomey Gap region (DHG). C, western African region. Black dots are landmark positions.

The additional LDA discriminated CEN from WES and DHG groups along LD1 (73.3%). LD2 (26.7%) discriminated CEN and WES from DHG groups. We tested the statistical significance of the intraspecific variation in the three groups identified above. A Procrustes ANOVA revealed a significant effect of the division of P. tricuspis into three morphological groups, WES, DHG and CEN (F2,67 = 4.23, P < 0.001; Table 2). Pairwise comparisons showed that all groups presented significantly different skull shapes (Supporting Information, Table S17). We calculated the mean shapes for each group (Supporting Information, Fig. S5). The WES skulls presented the shortest and widest rostrum, a relatively elongated infraorbital canal and the longest zygomatic process of the maxillary. The DHG skulls presented the widest nasals posteriorly, the most posterior projections of the maxilla into the palatine and the most posterior ventral margin of the foramen magnum. Skulls from the CEN region were characterized by the narrowest rostra, the smallest tympanic bullae and the most anterior petrosal–squamosal–exoccipital intersection.

The analyses on the allometry-corrected shapes revealed similar results (Supporting Information, Appendix S3; Fig. S6; Tables S12, S18, and S19).

Intraspecific variation in M. javanica

The analyses on the M. javanica dataset showed that differences in shape were solely explained by differences in size (Supporting Information, Figs S7, S8; Tables S12 and S20).

DISCUSSION

Size influences skull shape in extant pangolins

Size explained a significant part of the total variation in skull shape among species within the Pholidota (evolutionary allometry; Table 1; Fig. 4; Supporting Information, Fig. S2; Table S2). The HOS and PTA tests showed that the directions of ontogenetic trajectories were conserved among pangolins, with only M. javanica and M. pentadactyla differing significantly from S. gigantea (HOS) and P. tricuspis (PTA; Supporting Information, Tables S4 and S5). Ontogenetic trajectory angles did not differ within clades, which is consistent with observations that intraspecific (ontogenetic and/or static) allometric trajectories tend to differ as species divergence time increases (Voje et al., 2014; Esquerré et al., 2017).

Significant differences between the intercepts of multivariate regressions (overlap test) were found only between P. tricuspis and all other species and between P. tetradactyla and M. javanica. With the exception of the difference in intercept between the two Phataginus species, significant P-values were relatively close to 0.05 (i.e. same order of magnitude). This indicates that allometric trajectories are still relatively preserved across the Pholidota and suggests that cranial morphology is similar in early developmental stages. The heterochrony test showed that more than half of the pairwise comparisons performed (12 out of 19) revealed heterochronic changes. In fact, when the overlap test was performed taking size = 0 as reference, heterochrony was detected for all 19 comparisons (Supporting Information, Tables S6 and S7). This could suggest a major pattern of heterochrony driving the differentiation between the Asian and African clades, the first being putatively peramorphic by presenting longer rostra and longer zygomatic processes of the maxilla (traits associated with the PC1 of the predicted values; Figs 4, 5). However, given the absence of fetuses in our analyses, predictions of shape at minimum size (x = 0) should be considered with caution. Additionally, the heterogeneity of sampling between species might have influenced the significance of the overlap test, given that the largest differences between intercepts were not necessarily significant (Supporting Information, Table S6). The generally low P-values from the overlap test (P < 0.23) might reflect substantial differences in the intercepts, despite the non-significance yielded (Amrhein et al., 2019).

Moreover, African and Asian clades clearly presented non-overlapping ontogenetic allometric trajectories for traits correlated with PC1 of the predicted values of multivariate regression of shape on size (Fig. 4), with African species sharing a higher intercept relative to Asian ones. We interpret this size–shape space as a good representation of the evolutionary allometry. Considering that parameters of modelled growth trajectories can be used efficiently as continuous characters in phylogenies (Bardin et al., 2017), these different intercept values could theoretically constitute a valuable character to distinguish members of the Asian and African clades. In fact, they corroborate previous morphology- and DNA-based results (Gaubert et al., 2018) that described a split of the extant Manidae into two continental clades. The differences in the ontogenetic allometric patterns between the two size–shape spaces (multivariate vs. PC1 of the predicted values) might also suggest strong cranial modularity (i.e. rostrum module; Goswami, 2006). Cranial modules evolve semi-independently, and allometric patterns detected for highly integrated modules might differ from the overall pattern (Gerber & Hopkins, 2011). Further analyses are needed to confirm this hypothesis but would be beyond the scope of the present study.

Heterochronic changes are better illustrated by the ontogenetic allometric trajectories of the two Smutsia species (Fig. 4; Supporting Information, Tables S6–S8), as they overlap in both size–shape and shape spaces (Supporting Information, Fig. S1; Mitteroecker et al., 2005; Esquerré et al., 2017). Both PC1 of the predicted values and multivariate regressions suggest that S. temminckii is paedomorphic, because it tends to resemble juvenile S. gigantea (peramorphic). Nevertheless, heterochrony is not always associated with close phylogenetic affinities. Differing ontogenetic allometric trajectories between closely related species were reported previously in hominids (Mitteroecker et al., 2004). Despite being sister taxa, humans and chimpanzees differ in skull shape from early ontogenetic stages. The allometric trajectories of the two Phataginus species exhibit a similar pattern (Fig. 4; Supporting Information, Table S6). Both species show rather distinct intercepts (distinct shapes from early stages).

Overall, our results suggest that complex allometric changes played an important role in the morphological evolution of the pangolin skull. All pangolins follow a similar ontogenetic trend characterized by the elongation of the rostrum and a posterior projection of the zygomatic process of the maxilla (Figs 4, 5). As a consequence, the braincase is relatively smaller in larger species. These ontogenetic patterns are in line with described patterns of evolutionary allometry in which large-sized mammals evolve longer rostra (Cardini & Polly, 2013; Cardini et al., 2015; Cardini, 2019). Our allometric and phenotypic trajectories (Fig. 4; Supporting Information, Fig. S2), associated with the thin plate spline deformations, enable us to suggest that the ontogenetic drift of S. gigantea towards the Maninae (Asian pangolins) morphospace is attributable, in part, to the elongation of the skull associated with size (Fig. 5). Larger species with more elongated rostra (Smutsia gigantea and Manis) additionally present deep nasal notches (Supporting Information, Fig. S9).

Evolutionary patterns of ontogenetic allometry should therefore be taken into account in morphology-based studies. For instance, when looking at cranial character states in the morphology-based phylogeny from Gaudin et al. (2009), the depth of the nasal notch (character 306; Gaudin et al., 2009) appears to be associated with size. The small-sized Phataginus species are the only ones presenting shallow nasal notches (Supporting Information, Fig. S9). The relative length of the parietal–squamosal suture (character 385; Gaudin et al., 2009) also appears to be influenced by allometry (Figs 4, 5). The parietal–squamosal suture is relatively longer in skulls with shorter snouts, which is the case in Phataginus species (the smallest pangolins) that show a relatively long parietal–squamosal suture. This is to be expected, because skull length is influenced mostly by elongation of the rostrum in larger species (Figs 4, 5). According to Gaudin et al. (2009), S. temminckii is the only species presenting a multistate for this character [< 25% greatest skull length (0) or > 25% (1)], which is congruent with its short snout and intermediate average size between that of Phataginus spp. and that of Manis spp. and S. gigantea (Fig. 4). These results therefore call for a revision of some characters included in morphological matrices used to reconstruct extant and extinct pangolin phylogenetic relationships.

Patterns of skull shape variation support the classification of extant pangolins in three distinct genera

Based on morphological features, pangolins have been classified from a single genus (Manis) to a maximum of six different genera. In 1882, Jentink published a monograph on the comparative anatomy of extant pangolins, except for M. culionensis, in which he briefly referred to the remarkable differences between the skulls of the seven recognized species (Jentink, 1882). Nevertheless, he ascribed all pangolin species to the same genus, postponing a thorough investigation of skull morphology. A division of pangolins into six different genera was later proposed by Pocock (1924), with the African Smutsia (S. temminckii and S. gigantea), P. tricuspis, Uromanis longicaudata (Linnaeus, 1766) (P. tetradactyla) and the Asian M. pentadactyla, Paramanis javanica and Phatages crassicaudata. This classification was based on external/soft tissue characters only, but neglected cranial osteology.

Although the African–Asian split has been widely accepted, and the Asian genera have generally been merged into the single genus Manis, genus-level classifications have varied within the African clade. The four species have been ascribed either to a single genus, Phataginus (Patterson, 1978; Corbet & Hill, 1991), or kept separate in three distinct genera, as proposed by Pocock (1924; McKenna & Bell, 1997). Recent phylogenies based on morphological traits (Gaudin et al., 2009) and molecular data (Gaubert et al., 2018) have supported the view of an Asian–African split, with the Asian genus, Manis (Maninae), as the sister clade to the African pangolins composed of the two genera Phataginus (Phatagininae Gaubert, 2018) and Smutsia (Smutsiinae Gray, 1873).

Our results support the distinction of three extant pangolin genera, as recognized in previous works (Gaudin & Wible, 1999; Koenigswald, 1999; Gaudin et al., 2009; Gaubert et al., 2018). The PCAs and LDAs (Fig. 6A, B) largely reflect the division of extant Pholidota in Maninae (Asian; mostly positive PC1 values) and Smutsiinae + Phatagininae. The African clade was weakly supported in the most recent phylogeny based on anatomical characters (Gaudin et al., 2009). In all analyses excluding fossil taxa, Smutsia spp. are recovered as a sister group to Maninae (Gaudin et al., 2009). This might be explained, in part, by the detected allometric effect (see ‘Discussion: Size influences skull shape in extant pangolins’; Fig. 4; Supporting Information, Fig. S2). Nevertheless, both the PCA (Fig. 6B) and LDA (Fig. 6C, D) reveal the grouping of all Manis species to the exclusion of Smutsia and Phataginus. These analyses also show a clear separation between Smutsia and Phataginus, with the small pangolins showing lower values of PC1 and higher values of LD1.

We also show some degree of variation in skull shape at the intrageneric level that confirms species-level delineation of extant pholidotans (Fig. 6). Despite the split between S. gigantea and S. temminckii being slightly more recent than in other genera (5.6–13.2 Mya for Smutsia, 10.3–15.6 Mya for Manis and 9.3–16.5 Mya for Phataginus; Gaubert et al., 2018), their skulls appear comparatively more distinct. Some of the differences between S. temminckii and S. gigantea are related to size (see above; Fig. 4). The most extensive morphological phylogenetic work performed to date found that Smutsia was the least-supported modern genus, with only three unique unambiguous synapomorphies, none of which involve the cranium (Gaudin et al., 2009). Our results on variation in skull shape are congruent with the low support for the Smutsia node. As discussed above, elongation of the rostrum largely influences this intrageneric difference in shape. This allometric pattern is present at both ontogenetic and evolutionary levels and explains some of the differences in skull shape between S. gigantea and S. temminckii. Although the shape of the skull of S. gigantea is more similar to that of Asian pangolins, S. temminckii is closer to Phataginus (Fig. 6). In addition to these substantial differences in skull shape within Smutsia, previous molecular analyses reported a relatively large mitogenomic distance within the genus (11.9%; Gaubert et al., 2018), although lower than those reported within Phataginus (see below).

The two Phataginus species present the largest intrageneric mitogenomic distance (14.3%; Gaubert et al., 2017). This distance is patent in PC3 and LD2 scores, which clearly separate P. tetradactyla from P. tricuspis (Fig. 6B, C). However, the cranial shape is more similar between the two Phataginus species than between the two Smutsia species (Supporting Information, Table S11). Phataginus tetradactyla was previously ascribed to the genus Uromanis (Pocock, 1924), but recent cladistic analyses based on osteological characters yielded a strong support for placement in the genus Phataginus, the best supported among all genera (Gaudin et al., 2009). Of the seven unambiguous synapomorphies corresponding to cranial traits, the orientation and size of the zygomatic process of the squamosal (character 355; Gaudin et al., 2009) is coded as ventrally directed and short dorsoventrally for both species of tree pangolin. We confirmed this character state (Fig. 6A; PC1), but additionally found that the shape of this process in the horizontal direction constitutes one of the main differences between the two species, with P. tetradactyla presenting the longest among all pholidotans (Fig. 6B; PC3).

Although three genera have been proposed in Maninae, recent studies have suggested that the three species should be grouped into the single genus, Manis (Gaudin et al., 2009; Gaubert et al., 2018). On average, Manis is the genus with the lowest (but still high) mitogenomic distance among species (mean = 9.3%), with M. javanicaM. culionensis showing the lowest value (3.1%) and M. pentadactylaM. javanica/M. culionensis/ M. crassicaudata showing the highest (12.2%; Gaubert et al., 2017). The three species of Maninae show some overlap in morphospace, but are well segregated by PC2 (Fig. 6A). In contrast, M. culionensis overlaps with M. javanica in morphospace (Supporting Information, Fig. S1), which is congruent with the low mitogenomic distance and the recent divergence time estimated between these two species (0.4–2.5 Mya). However, our data are clearly insufficient to assess the morphological discrimination between M. javanica and M. culionensis (Supporting Information, Fig. S1).

Despite the strength of the Maninae node, infrageneric relationships have been greatly debated. Molecular-based analyses show a well-supported node (Bayesian posterior probability = 1), including M. javanica/M. culionensis and M. crassicaudata, with M. pentadactyla as the sister taxon (Gaubert et al., 2018). The M. crassicaudataM. javanica/M. culionensis clade remains moderately supported in the Bayesian inference excluding mitogenomes (nuclear DNA only). In contrast, a moderately supported sister-group relationship between M. pentadactyla and M. crassicaudata is retrieved by morphology-based phylogenetic analyses (bootstrap value = 76; Gaudin et al., 2009). The list of synapomorphies from Gaudin et al. (2009) for the node gathering M. crassicaudata and M. pentadactyla featured only three cranial traits. Among these, only the position of the foramen ovale at the level of the anterior edge of the ectotympanic is an unambiguous synapomorphy [character 379(1)]. Landmarks 18/43 and 20/45 describe this synapomorphy and contribute to segregate M. crassicaudata and M. pentadactyla from M. javanica along with other traits correlated with PC2 (Fig. 6A). Further anatomical investigation (i.e. internal characters) remains necessary to explore the morphological support for both hypotheses more thoroughly.

Skull shape variation corroborates cryptic phylogeographical lineages in P. tricuspis

Cranial osteological data are extremely useful to unveil patterns of cryptic speciation (Sukumaran & Knowles, 2017). Our results confirm, at least in part, the existence of several geographical groupings within P. tricuspis (Gaubert et al., 2016, 2018; Figs 2, 7; Supporting Information, Fig. S5). The PCA showed some degree of overlap between the different lineages in morphospace (Fig. 7A, B), although we were able to find significant differences among the skulls of four lineages or lineage groups (Western Africa + Ghana, Dahomey Gap, Western Central Africa and Central Africa). Homogeneity in skull shape might be explained by recent divergence times between these lineages (0.8–4.6 Mya for the most recent common ancestor of Western and Central African lineages; Gaubert et al., 2016, 2018). Nevertheless, we showed that two lineages (WAF and GHA) from a molecularly identified western lineage (DHG–GHA–WAF; Gaubert et al., 2016) present a distinct cranial morphology (Figs 7C, D, 8). Furthermore, the Dahomey Gap lineage presents a distinct skull shape. Although the identification of group membership was not particularly high among cryptic lineages (75.4%), we showed that WAF–GHA, DHG and CAF–WCA form three morphologically distinct groups separated along a longitudinal gradient, with a high rate of a posteriori attributions (95.7%; Fig. 8).

Our results recover only partial evidence of the three nuclear groups found by Gaubert et al. (2016) that delineated Western Africa (WAF, GHA and DHG), Western Central Africa (WCA) and Central Africa (CA). Although the grouping of WAF–GHA is congruent with the molecular phylogeny, the Dahomey Gap group is morphologically divergent in this clade. In addition, the similarity between Central Africa and Western Central Africa groups could suggest that they both retained a plesiomorphic cranial shape. However, testing the congruence between molecular and morphological data would require further scrutiny.

In contrast, we posit that different environmental conditions might also explain some parts of the cranial variation among pangolin lineages. The Dahomey Gap corresponds to a savannah-like corridor that divides western and central lowland rain forests (Dupont & Weinelt, 1996), with a longitudinal gradient ranging from dry (Dahomey Gap) to more humid rainforest conditions (western and central rainforest). The adaptation to a drier climate resulting from tropical forest fragmentation (Salzmann & Hoelzmann, 2005) could explain the differentiation of a Dahomey Gap lineage. Concurrently, genetic drift resulting from an isolation-induced reduction of gene flow (Renaud & Millien, 2001) might also have played a role in the Dahomey Gap skull shape differentiation, following a vanishing refuge model of diversification (Vanzolini & Williams, 1981; Damasceno et al., 2014; Gaubert et al., 2016). At least three other endemic mammal species/subspecies have recently been described or confirmed on a genetic basis in the Dahomey Gap (Colyn et al., 2010; Nicolas et al., 2010; Houngbédji et al., 2012). Further analyses are needed to assess the potential interaction between variation in skull shape and environmental conditions in pangolins.

Intraspecific variation in M.javanica

In contrast to P. tricuspis, we did not find solid evidence of skull shape discrimination between molecularly identified lineages within M. javanica, because differences in shape appear to be associated solely with size (Supporting Information, Appendix S3; Figs S7, S8). The lack of differences in skull shape might be explained by introgressions between lineages or by more recent divergence times (Nash et al., 2018) than among P. tricuspis lineages. Additionally, the cryptic lineages within M. javanica might currently lack a well-defined geographical delimitation, attributable, in part, to the lack of precise geographical information for the wild specimens sampled (Nash et al., 2018). Human-induced specimen translocation by the introduction of pangolin seizures of unknown origin might also have resulted in the mixing of different lineages (Pantel & Chin, 2009).

Conclusion

Our results are congruent with the currently accepted genus- and species-level classification of extant pangolins. We found that heterochronic changes explain, in part, the morphological differentiation of the skull at an intrageneric level. However, some species appear to present different allometric trajectories resulting from changes in skull shape during early developmental stages. Asian and African clades can be discriminated on the basis of the allometric trajectories of traits related to PC1. Namely, we identified elongation of the rostrum to be related to ontogenetic allometry, and we hypothesized that this might be also present at the evolutionary level. This might explain the detected differences in rostral proportions between species of different sizes and the apparent morphological convergence between S. gigantea and Manis species. Our results underline the importance of accounting for allometry when performing phylogenetic analyses based on morphological characters.

Our results also show that skull shape differs between cryptic lineages within P. tricuspis, and that these can be circumscribed into three geographical groups from western Africa (WAF–GHA), the Dahomey Gap (DHG) and Central Africa (CAF–WCA). We show that skull shape is potentially useful to determine pangolin species identity and, at a finer scale, the geographical origin of specimens of white-bellied pangolins seized in illegal trade hubs or markets. Such information could help to determine differential poaching pressures, delimitate management units, and thus refine threat status at a regional level.

SUPPORTING INFORMATION

Additional supporting information may be found in the online version of this article at the publisher’s web-site:

Appendix S1. Specimen list and landmark coordinates.

Appendix S2. Discrete traits.

Appendix S3. Additional analyses.

Appendix S4. LDA a posteriori attribution tables.

Table S1. Definitions of the 75 cranial landmarks used.

Table S2. Phylogenetic ANOVA of shape (Procrustes coordinates) ~ log(centroid size). *Significant P-value. The randomized residual permutation procedure used 10 000 permutations.

Table S3. Pairwise comparisons of the allometric trajectory angles for skull shape. P-values (10.000 iterations) are in the upper triangle and angles between slopes (in degrees) in the lower triangle. Significant results are written in bold.

Table S4. Pairwise comparisons of the phenotypic trajectory angles. P-values (10 000 iterations) are in the upper triangle and angles between path distances in the lower triangle. Significant results are written in bold.

Table S5. Pairwise comparisons of the phenotypic trajectory lengths. P-values (10 000 iterations) are in the upper triangle and absolute differences between path distances in the lower triangle. Significant results are written in bold.

Table S6. Pairwise comparison of the intercept (x = 5.125) of the multivariate Procrustes allometric regressions (overlap test). P-values (Hochberg corrected) of the difference between the intercepts computed with 10 000 iterations are in the upper triangle and observed differences in the lower triangle. Taxa with significantly non-parallel trajectories were not included in the test. Significant results are written in bold.

Table S7. Pairwise comparison of the intercept (x = 0) of the multivariate Procrustes allometric regressions (overlap test). P-values of the difference between the intercepts computed with 10 000 iterations are in the upper triangle and observed differences in the lower triangle. Taxa with significantly non-parallel trajectories were not included in the test. Significant results are written in bold.

Table S8. Pairwise comparisons of the predicted head shape (predicted Procrustes residuals) differences at maximum centroid size (heterochrony test). P-values (Hochberg corrected; in black for x = 5.125 and in blue for x = 0) of the difference between them computed with 10 000 iterations are in the upper triangle and observed differences in the lower triangle. Taxa with significantly non-parallel trajectories were not included in the test. Significant results are written in bold.

Table S9. Pairwise comparison of the principal component 1 (PC1) allometric trajectory intercepts for skull shape. P-values of a multiple comparison of means (Tukey’s test) are in the upper triangle and the t-statistics values in the lower triangle. The same intercept in all species is the null hypothesis. The P-values are indicated for species that show the same intercept; *all other pairwise comparisons retrieved a P-value < 0.001. Only species with parallel slopes are included. Significant results are written in bold.

Table S10. ANOVA of shape (Procrustes coordinates) ~ sex*species and shape (Procrustes coordinates) ~ size + sex*species. *Significant P-value. The randomized residual permutation procedure used 10 000 permutations.

Table S11. Pairwise comparison of the Procrustes distances between least-squares (LS) means for species. P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S12. ANOVA of shape ~ taxa/geographical groups of adult specimens of the interspecific, Phataginus tricuspis and Manis javanica datasets. Significant P-values indicate differences between skull shapes of taxa/geographical groups; n groups is the number of taxa or geographical groups used as factors. *Significant P-value. The randomized residual permutation procedure used 10 000 permutations.

Table S13. Pairwise comparison of the allometry-corrected Procrustes distances between least-squares (LS) means for species. P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S14. Significant t-tests of principal components vs. log-transformed centroid size.

Table S15. Pairwise comparisons of the allometric trajectory angles for skull shape in Phataginus tricuspis cryptic lineages [homogeneity of slopes (HOS)]. P-values (10 000 iterations) are in the upper triangle and angles between slopes (in degrees) in the lower triangle. Significant results are written in bold.

Table S16. Pairwise comparison of the Procrustes distances between least-squares (LS) means for cryptic lineages (Gaubert et al., 2016). P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S17. Pairwise comparison of the Procrustes distances between least-squares (LS) means for three geographical groups. P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S18. Pairwise comparison of the allometry-corrected Procrustes distances between least-squares (LS) means for cryptic lineages (Gaubert et al., 2016). P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S19. Pairwise comparison of the allometry-corrected Procrustes distances between least-squares (LS) means for three geographical groups. P-values of the difference between the LS means computed with 10 000 iterations are in the upper triangle and observed distances in the lower triangle. Significant results are written in bold.

Table S20. ANOVA of shape (Procrustes coordinates) ~ log(centroid size)*cryptic lineages (N = 35) in M. javanica. *Significant P-value.

Figure S1. Principal components analysis on the dataset including juveniles and two specimens of Manis culionensis (N = 243).

Figure S2. Phenotypic trajectory analysis among seven pangolin species. The morphospace delimited by principal component (PC)1 and PC2 explaining the variance between adults and juveniles within each species (N = 125) is shown. For each species, a trajectory representing the change in shape between the shape estimates for juveniles (light grey circles) and adults (dark grey circles) is represented. Deformed meshes represent the maximum and minimum shapes of PC1 and PC2.

Figure S3. Principal component (A, PC1 vs. PC2; B, PC1 vs. PC3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3), with associated allometry-corrected variation in shape, for crania of seven pangolin species (N = 173). Shapes are the residuals of a pooled within-group multivariate regression of shape on log-transformed centroid size.

Figure S4. Allometric trajectories of the cryptic lineages (Gaubert et al., 2016) of Phataginus tricuspis (N = 95). The x-axis values are the log-transformed centroid sizes for each specimen; the y-axis values are the principal component 1 of the predicted values of multivariate regression of shape ratios on size. The size of the dots indicates the size of the specimens. Abbreviation: CAF, Central Africa; DHG, Dahomey Gap; GHA, Ghana; WAF, Western Africa; WCA, Western Central Africa. GAB (Gabon) was excluded because our dataset included only one skull.

Figure S5. Mean shapes of proposed management units for the Phataginus tricuspis sample (N = 70) in lateral (left), ventral (middle) and dorsal (right) views. A, Central African region. B, Dahomey Gap region. C, Western African region. Black dots are landmark positions.

Figure S6. Principal components (A, PCres1 vs. PCres2; B, PCresC1 vs. PCres3) and linear discriminant analyses (C, LD1 vs. LD2; D, LD1 vs. LD3) with associated allometry-corrected variation in shape for crania of six cryptic lineages of Phataginus tricuspis (N = 71; Gaubert et al., 2016). Shapes are the residuals of a pooled within-group multivariate regression of shape on log-transformed centroid size.

Figure S7. Allometric trajectories of the lineages (Nash et al., 2018) of Manis javanica (N = 35). The x-axis values are the log-transformed centroid sizes for each specimen; y-axis values are the principal component 1 of the predicted values of multivariate regression of shape ratios on size. The size of the dots indicates the size of the specimens. Abbreviations: BOR, Borneo; JAV, Java; SUM/SIN, Sumatra/Singapore.

Figure S8. Principal components (A, PC1 vs. PC2; B, PC1 vs. PC3) with associated variation in shape for crania of Manis javanica lineages (N = 25; Nash et al., 2018).

Figure S9. Minimum (left) and maximum (right) shape prediction from a multivariate regression on log-transformed centroid size for two species of small (A, B) and two species of large (C, D) pangolins in dorsal view. A, Phataginus tricuspis. B, Phataginus tetradactyla. C, Smutsia gigantea. D, Manis javanica. Grey and red dots mark landmark positions at minimum and maximum sizes, respectively. Black dots are landmark positions for maximum size predictions.

ACKNOWLEDGEMENTS

We thank Steffen Bock, Frieder Mayer, Detlef Willborn (MfN), Roberto Portela Miguez, Louise Tomsett (BMNH), Géraldine Veron (MNHN), Chris Conroy (MVZ), Darrin Lunde, Paula Bohaska, John Ososky (USNM) and Eleanor Hoeger (AMNH) for access to collections. We acknowledge Farah Ahmed, Amin Garbout (BMNH) and Renaud Lebrun (ISEM) for assistance with micro-computed tomography scanning. We also thank Pierre-Henri Fabre, Quentin Martinez, and Julien Claude (ISEM) for fruitful discussions. Finally, we thank two anonymous reviewers for helpful comments on the manuscript. S.F.-C., L.H. and F.D. were supported by a European Research Council (ERC) consolidator grant (ConvergeAnt #683257). L.H. and F.D. were supported by Centre National de la Recherche Scientifique (CNRS). This research received support from the Synthesys Project, https://www.synthesys.info/, financed by the European Community Research Infrastructure Action under the FP7 (GB-TAF-5606 and BE-TAF-5661). This work was supported by PANGO-GO (ANR-17-CE02-0001) and by ‘Investissements d’Avenir’ grants managed by Agence Nationale de la Recherche Labex CEMEB (ANR-10-LABX-0004), Labex NUMEV (ANR-10-LABX-0020). This is contribution ISEM 2019-150 of the Institut des Sciences de l’Évolution de Montpellier.

REFERENCES

Adams
DC
,
Collyer
ML
.
2009
.
A general framework for the analysis of phenotypic trajectories in evolutionary studies
.
Evolution
63
:
1143
1154
.
Adams
D
,
Collyer
M
,
Kaliontzopoulou
A
,
Sherratt
E
.
2017
.
Geomorph: software for geometric morphometric analyses
. R package, Version 3.0.5. Available at: https://cran.r-project.org/package=geomorph.
Adams
DC
,
Nistri
A
.
2010
.
Ontogenetic convergence and evolution of foot morphology in European cave salamanders (Family: Plethodontidae)
.
BMC Evolutionary Biology
10
:
216
.
Amrhein
V
,
Greenland
S
,
McShane
B
.
2019
.
Scientists rise up against statistical significance
.
Nature
567
:
305
307
.
Bardin
J
,
Rouget
I
,
Cecca
F
.
2017
.
Ontogenetic data analyzed as such in phylogenies
.
Systematic Biology
66
:
23
37
.
Benítez
HA
,
Bravi
R
,
Parra
LE
,
Sanzana
MJ
,
Sepúlveda-Zúñiga
E
.
2013
.
Allometric and non-allometric patterns in sexual dimorphism discrimination of wing shape in Ophion intricatus: might two male morphotypes coexist?
Journal of Insect Science (Online)
13
:
143
.
Bickford
D
,
Lohman
DJ
,
Sodhi
NS
,
Ng
PK
,
Meier
R
,
Winker
K
,
Ingram
KK
,
Das
I
.
2007
.
Cryptic species as a window on diversity and conservation
.
Trends in Ecology & Evolution
22
:
148
155
.
Cardini
A
.
2019
.
Craniofacial allometry is a rule in evolutionary radiations of placentals
.
Evolutionary Biology
46
:
239
248
.
Cardini
A
,
O’Higgins
P
.
2004
.
Patterns of morphological evolution in Marmota (Rodentia, Sciuridae): geometric morphometrics of the cranium in the context of marmot phylogeny, ecology and conservation
.
Biological Journal of the Linnean Society
82
:
385
407
.
Cardini
A
,
Polly
PD
.
2013
.
Larger mammals have longer faces because of size-related constraints on skull form
.
Nature Communications
4
:
2458
.
Cardini
A
,
Polly
D
,
Dawson
R
,
Milne
N
.
2015
.
Why the long face? Kangaroos and wallabies follow the same ‘rule’ of cranial evolutionary allometry (CREA) as placentals
.
Evolutionary Biology
42
:
169
176
.
Challender
DWS
,
Waterman
C
,
Baillie
JEM
.
2014
.
Scaling up pangolin conservation. IUCN SSC Pangolin Specialist Group Conservation Action Plan.
London: Zoological Society of London
.
Cheverud
JM
.
1982
.
Relationships among ontogenetic, static, and evolutionary allometry
.
American Journal of Physical Anthropology
59
:
139
149
.
Collyer
ML
,
Adams
DC
.
2013
.
Phenotypic trajectory analysis: comparison of shape change patterns in evolution and ecology
.
Hystrix
24
:
75
83
.
Colyn
M
,
De Rennes
U
,
Hulselmans
J
,
Oude
P
.
2010
.
Discovery of a new duiker species (Bovidae: Cephalophinae) from the Dahomey Gap, West Africa
.
Zootaxa
2637
:
1
30
.
Corbet
GB
,
Hill
JE
.
1991
.
A world list of mammalian species
.
London
:
Natural History Museum Publishing and Oxford University Press
.
Damasceno
R
,
Strangas
ML
,
Carnaval
AC
,
Rodrigues
MT
,
Moritz
C
.
2014
.
Revisiting the vanishing refuge model of diversification
.
Frontiers in Genetics
5
:
353
.
Dryden
IL
,
Mardia
KV
.
1993
.
Multivariate shape analysis
.
Sankhyā: The Indian Journal of Statistics, Series A
55
:
460
480
.
Dupont
LM
,
Weinelt
M
.
1996
.
Vegetation history of the savanna corridor between the Guinean and the Congolian rain forest during the last 150,000 years
.
Vegetation History and Archaeobotany
5
:
273
292
.
Emry
RJ
.
1970
.
A North American Oligocene pangolin and other additions to the Pholidota
.
Bulletin of the American Museum of Natural History
142
:
455
510
.
Emry
RJ
.
2004
.
The edentulous skull of the North American pangolin
,
Patriomanis americanus. Bulletin of the American Museum of Natural History
285
:
130
138
.
Esquerré
D
,
Sherratt
E
,
Keogh
JS
.
2017
.
Evolution of extreme ontogenetic allometric diversity and heterochrony in pythons, a clade of giant and dwarf snakes
.
Evolution
71
:
2829
2844
.
Evin
A
,
Cucchi
T
,
Cardini
A
,
Strand Vidarsdottir
U
,
Larson
G
,
Dobney
K
.
2013
.
The long and winding road: identifying pig domestication through molar size and shape
.
Journal of Archaeological Science
40
:
735
743
.
Ferreira-Cardoso
S
,
Delsuc
F
,
Hautier
L
.
2019
.
Evolutionary tinkering of the mandibular canal linked to convergent regression of teeth in placental mammals
.
Current Biology
29
:
468
475.e3
.
Foth
C
,
Hedrick
BP
,
Ezcurra
MD
.
2016
.
Cranial ontogenetic variation in early saurischians and the role of heterochrony in the diversification of predatory dinosaurs
.
PeerJ
4
:
e1589
.
Galatius
A
,
Kinze
CC
,
Teilmann
J
.
2012
.
Population structure of harbour porpoises in the Baltic region: evidence of separation based on geometric morphometric comparisons
.
Journal of the Marine Biological Association of the United Kingdom
92
:
1669
1676
.
Gaubert
P
.
2011
.
Family Manidae
. In:
Wilson
D
,
Mittermeier
R
, eds.
Handbook of the mammals of the world. Vol. 2. Hoofed mammals.
Barcelona
:
Lynx Edicions
,
82
103
.
Gaubert
P
,
Antunes
A
,
Meng
H
,
Miao
L
,
Peigné
S
,
Justy
F
,
Njiokou
F
,
Dufour
S
,
Danquah
E
,
Alahakoon
J
,
Verheyen
E
,
Stanley
WT
,
O’Brien
SJ
,
Johnson
WE
,
Luo
SJ
.
2018
.
The complete phylogeny of pangolins: scaling up resources for the molecular tracing of the most trafficked mammals on earth
.
The Journal of Heredity
109
:
347
359
.
Gaubert
P
,
Njiokou
F
,
Ngua
G
,
Afiademanyo
K
,
Dufour
S
,
Malekani
J
,
Bi
SG
,
Tougard
C
,
Olayemi
A
,
Danquah
E
,
Djagoun
CA
,
Kaleme
P
,
Mololo
CN
,
Stanley
W
,
Luo
SJ
,
Antunes
A
.
2016
.
Phylogeography of the heavily poached African common pangolin (Pholidota, Manis tricuspis) reveals six cryptic lineages as traceable signatures of Pleistocene diversification
.
Molecular Ecology
25
:
5975
5993
.
Gaudin
TJ
,
Emry
RJ
,
Wible
JR
.
2009
.
The phylogeny of living and extinct pangolins (Mammalia, Pholidota) and associated taxa: a morphology based analysis
.
Journal of Mammalian Evolution
16
:
235
305
.
Gaudin
TJ
,
Gaubert
P
,
Billet
G
,
Hautier
L
,
Ferreira-Cardoso
S
,
Wible
JR
.
2019
.
Evolution & morphology
. In:
Challender
DWS
,
Nash
H
,
Waterman
C
, eds.
Pangolins: science, society and conservation.
Cambridge
:
Academic Press
.
Gaudin
TJ
,
Wible
JR
.
1999
.
The entotympanic of pangolins and the phylogeny of the Pholidota (Mammalia)
.
Journal of Mammalian Evolution
6
:
39
65
.
Gerber
S
,
Hopkins
MJ
.
2011
.
Mosaic heterochrony and evolutionary modularity: the trilobite genus Zacanthopsis as a case study
.
Evolution
65
:
3241
3252
.
Goswami
A
.
2006
.
Cranial modularity shifts during mammalian evolution
.
The American Naturalist
168
:
270
280
.
Gray
JA
,
Sherratt
E
,
Hutchinson
MN
,
Jones
MEH
.
2019
.
Changes in ontogenetic patterns facilitate diversification in skull shape of Australian agamid lizards
.
BMC Evolutionary Biology
19
:
7
.
Gunz
P
,
Mitteroecker
P
,
Neubauer
S
,
Weber
GW
,
Bookstein
FL
.
2009
.
Principles for the virtual reconstruction of hominin crania
.
Journal of Human Evolution
57
:
48
62
.
Hatt
RT
,
Lang
H
,
Chapin
JP
.
1934
.
The pangolins and aard-varks collected by the American Museum Congo Expedition
.
Bulletin of the American Museum of Natural History
66
:
643
672
.
Hautier
L
,
Billet
G
,
Eastwood
B
,
Lane
J
.
2014
.
Patterns of morphological variation of extant sloth skulls and their implication for future conservation efforts
.
The Anatomical Record
297
:
979
1008
.
Hautier
L
,
Billet
G
,
de Thoisy
B
,
Delsuc
F
.
2017
.
Beyond the carapace: skull shape variation and morphological systematics of long-nosed armadillos (genus Dasypus)
.
PeerJ
5
:
e3650
.
Hebert
PDN
,
Penton
EH
,
Burns
JM
,
Janzen
DH
,
Hallwachs
W
.
2004
.
Species realities and numbers in sexual vertebrates: perspectives from an asexually transmitted genome
.
Proceedings of the National Academy of Sciences of the United States of America
96
:
992
995
.
von Helversen
O
,
Heller
KG
,
Mayer
F
,
Nemeth
A
,
Volleth
M
,
Gombkötö
P
.
2001
.
Cryptic mammalian species: a new species of whiskered bat (Myotis alcathoe n. sp.) in Europe
.
Die Naturwissenschaften
88
:
217
223
.
Houngbédji
MG
,
Djossa
BA
,
Adomou
AC
,
Dakpogan
SC
,
Sinsin
B
,
Mensah
GA
.
2012
.
Conservation status of the red-bellied guenon (Cercopithecus erythrogaster erythrogaster) in the western Dahomey Gap in southwestern Benin and the adjacent Togodo Forest Reserve, south Togo
.
African Primates
7
:
184
192
.
Jentink
FA
.
1882
.
Revision of the Manidae in the Leyden Museum
.
Notes from the Leyden Museum
4
:
193
209
.
Klingenberg
CP
.
2011
.
MorphoJ: an integrated software package for geometric morphometrics
.
Molecular Ecology Resources
11
:
353
357
.
Klingenberg
CP
.
2013
.
Cranial integration and modularity: insights into evolution and development from morphometric data
.
Hystrix
24
:
43
58
.
Klingenberg
CP
.
2016
.
Size, shape, and form: concepts of allometry in geometric morphometrics
.
Development Genes and Evolution
226
:
113
137
.
McKenna
MC
,
Bell
SK
.
1997
.
Classification of mammals above the species level.
New York
:
Columbia University Press
.
Miranda
FR
,
Casali
DM
,
Perini
FA
,
Machado
FA
,
Santos
FR
.
2018
.
Taxonomic review of the genus Cyclopes Gray, 1821 (Xenarthra: Pilosa), with the revalidation and description of new species
.
Zoological Journal of the Linnean Society
183
:
687
721
.
Mitteroecker
P
,
Gunz
P
,
Bernhard
M
,
Schaefer
K
,
Bookstein
FL
.
2004
.
Comparison of cranial ontogenetic trajectories among great apes and humans
.
Journal of Human Evolution
46
:
679
697
.
Mitteroecker
P
,
Gunz
P
,
Bookstein
FL
.
2005
.
Heterochrony and geometric morphometrics: a comparison of cranial growth in Pan paniscus versus Pan troglodytes
.
Evolution & Development
7
:
244
258
.
Moraes-Barros
N
,
Miyaki
CY
,
Morgante
JS
.
2007
.
Identifying management units in non-endangered species: the example of the sloth Bradypus variegatus Schinz, 1825
.
Brazilian Journal of Biology
64
:
829
837
.
Murphy
W
,
Eizirik
E
,
Johnson
W
,
Zhang
Y
.
2001a
.
Molecular phylogenetics and the origins of placental mammals
.
Nature
409
:
614
618
.
Murphy
W
,
Eizirik
E
,
O’Brien
S
,
Madsen
O
,
Scally
M
,
Douady
C
,
Teeling
E
,
Ryder
O
,
Stanhope
M
,
De Jong
W
,
Springer
M
.
2001b
.
Resolution of the early placental mammal radiation using Bayesian phylogenetics
.
Science
294
:
2348
2351
.
Nash
HC
,
Low
GW
,
Choo
SW
,
Chong
JL
,
Semiadi
G
,
Hari
R
,
Sulaiman
MH
,
Turvey
TS
,
Evans
TA
,
Rheindt
FE
.
2018
.
Conservation genomics reveals possible illegal trade routes and admixture across pangolin lineages in Southeast Asia
.
Conservation Genetics
19
:
1083
1095
.
Nicolas
V
,
Olayemi
A
,
Wendelen
W
,
Colyn
M
.
2010
.
Mitochondrial DNA and morphometrical identification of a new species of Hylomyscus (Rodentia: Muridae) from West Africa
.
Zootaxa
2579
:
30
44
.
Pantel
S
,
Chin
SY
.
2009
.
Proceedings of the workshop on trade and conservation of pangolins native to South and Southeast Asia
.
Singapore
:
TRAFFIC Southeast Asia
.
Patterson
B
.
1978
.
Pholidota and Tubulidentata.
In:
Maglio
VJ
,
Cooke
HBS
, eds.
Evolution of African Mammals
.
Cambridge
:
Harvard University Press
,
268
278
.
Piras
P
,
Salvi
D
,
Ferrara
G
,
Maiorino
L
,
Delfino
M
,
Pedde
L
,
Kotsakis
T
.
2011
.
The role of post-natal ontogeny in the evolution of phenotypic diversity in Podarcis lizards
.
Journal of Evolutionary Biology
24
:
2705
2720
.
Pocock
RI
.
1924
.
The external characters of the pangolins (Manidae)
.
Proceedings of the Zoological Society of London
94
:
707
723
.
R Development Core Team
.
2013
.
R: a language and environment for statistical computing
.
Renaud
S
,
Millien
V
.
2001
.
Intra- and interspecific morphological variation in the field mouse species Apodemus argenteus and A. speciosus in the Japanese archipelago: the role of insular isolation and biogeographic gradients
.
Biological Journal of the Linnean Society
74
:
557
569
.
Rohlf
F
,
Slice
D
.
1999
.
Extensions of the Procrustes method for the optimal superimposition of landmarks
.
Systematic Biology
39
:
40
59
.
Salzmann
U
,
Hoelzmann
P
.
2005
.
The Dahomey Gap: an abrupt climatically induced rain forest fragmentation in West Africa during the late Holocene
.
The Holocene
15
:
190
199
.
Schlager
S
.
2017
.
Morpho and Rvcg – shape analysis in R
. In:
Zheng
G
,
Li
S
,
Székely
G
, eds.
Statistical shape and deformation analysis
.
London
:
Academic Press
,
217
256
.
Sidlauskas
B
,
Mol
J
,
Vari
R
.
2011
.
Dealing with allometry in linear and geometric morphometrics: a taxonomic case study in the Leporinus cylindriformis group (Characiformes: Anostomidae) with description of a new species from Suriname
.
Zoological Journal of the Linnean Society
162
:
103
130
.
Sukumaran
J
,
Knowles
LL
.
2017
.
Multispecies coalescent delimits structure, not species
.
Proceedings of the National Academy of Sciences of the United States of America
114
:
1607
1612
.
Sveegaard
S
,
Galatius
A
,
Dietz
R
,
Kyhn
L
,
Koblitz
JC
,
Amundin
M
,
Nabe-Nielsen
J
,
Sinding
MHS
,
Andersen
LW
,
Teilmann
J
.
2015
.
Defining management units for cetaceans by combining genetics, morphology, acoustics and satellite tracking
.
Global Ecology and Conservation
3
:
839
850
.
Vanzolini
P
,
Williams
E
.
1981
.
The vanishing refuge: a mechanism for ecogeographic speciation
.
Papéis Avulsos de Zoologia
34
:
251
255
.
Venables
WN
,
Ripley
BD
.
2002
.
Modern applied statistics with S
.
New York
:
Springer
.
Villemant
C
,
Simbolotti
G
,
Kenis
M
.
2007
.
Discrimination of Eubazus (Hymenoptera, Braconidae) sibling species using geometric morphometrics analysis of wing venation
.
Systematic Entomology
32
:
625
634
.
Voje
KL
,
Hansen
TF
,
Egset
CK
,
Bolstad
GH
,
Pélabon
C
.
2014
.
Allometric constraints and the evolution of allometry
.
Evolution
68
:
866
885
.
Von Koenigswald
W
.
1999
.
Order Pholidota
. In:
Rössner
G
,
Heissig
K
, eds.
The Miocene land mammals of Europe
.
Munich
:
Dr Friedrich Pfeil
,
75
80
.
Zhang
H
,
Miller
M
,
Yang
F
,
Chan
HK
,
Gaubert
P
,
Ades
G
,
Fischer
GA
.
2015
.
Molecular tracing of confiscated pangolin scales for conservation and illegal trade monitoring in Southeast Asia
.
Global Ecology and Conservation
4
:
414
422
.
Zhou
ZM
,
Zhou
Y
,
Newman
C
,
Macdonald
DW
.
2014
.
Scaling up pangolin protection in China
.
Frontiers in Ecology and the Environment
12
:
97
98
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com