Cranial variation in allactagine jerboas (Allactaginae, Dipodidae, Rodentia): a geometric morphometric study.

Allactaginae is a subfamily of dipodids consisting of four- and five-toed jerboas ( Allactaga, Allactodipus, Orientallactaga, Pygeretmus, Scarturus) found in open habitats of Asia and North Africa. Recent molecular phylogenies have upended our understanding of this group's systematics across taxonomic scales. Here, I used cranial geometric morphometrics to examine variation across 219 specimens of 14 allactagine species ( Allactaga major, A. severtzovi, Orientallactaga balikunica, O. bullata, O. sibirica, Pygeretmus platyurus, P. pumilio, P. shitkovi, Scarturus aralychensis, S. euphraticus, S. hotsoni, S. indicus, S. tetradactylus, S. williamsi) in light of their revised taxonomy. Results showed no significant sexual size or shape dimorphism. Species significantly differed in cranial size and shape both overall and as species pairs. Species identity had a strong effect on both cranial size and shape. Only a small part of cranial shape variation was allometric, with no evidence of unique species allometries, and most specimens fit closely to the common allometric regression vector. Allactaga was the largest, followed by Orientallactaga, Scarturus, and finally Pygeretmus. Principal component 1 (PC1) separated O. bullata+ O. balikunica+ S. hotsoni (with inflated bullae along with reduced zygomatic arches and rostra) from A. major+ A. severtzovi+ O. sibirica (with converse patterns), while PC2 differentiated Orientallactaga (with enlarged cranial bases and rostra along with reduced zygomatic arches and foramina magna) from Scarturus+ Pygeretmus (with the opposite patterns). Clustering based on the unweighted pair group method with arithmetic mean (UPGMA) contained the four genera, but S. hotsoni clustered with O. bullata+ O. balikunica and O. sibirica clustered with A. major+ A. severtzovi, likely due to convergence and allometry, respectively.

INTRODUCTION

Allactaginae Vinogradov, 1925 is a subfamily of four- and five-toed jerboas and is currently divided into five genera (represented by at least 16 species): Allactaga F. Cuvier, 1836, Allactodipus Kolesnikov, 1937, Orientallactaga Shenbrot, 1984, Pygeretmus Gloger, 1841, and Scarturus Gloger, 1841 (Michaux & Shenbrot, 2017). Allactagines are arid adapted. They use their chisel-shaped incisors, sharp-clawed forelegs, and well-developed nasal disk to dig burrows, which they inhabit during the day but leave at night to forage for seeds, vegetation, roots, and insects with the aid of their large eyes and ears and habitat-matching pelage color (likely acting as camouflage) (Lebedev et al., 2013; Miljutin, 2008; Michaux & Shenbrot, 2017; Nowak & Paradiso, 1983; Shenbrot et al., 1999, 2008).

Allactagines and other jerboas (i.e., three-toed, pygmy, and long-eared jerboas) are members of the family Dipodidae Fischer de Waldheim, 1817, and are specialized for bipedal hopping (Michaux & Shenbrot, 2017). Although this gait originated in humid, forested environments, it is now a hallmark adaptation to arid, open environments as all modern jerboas and most other bipedal rodents occur in such regions (Alhajeri, 2016; Berman, 1985; Lebedev et al., 2013; Mares, 1975; Nowak & Paradiso, 1983; Wu et al., 2014). Jerboas achieve bipedalism through hindlimb elongation and associated modification in muscular morphology, central foot bone fusion (commonly) forming a strong “cannon bone”, lateral digit reduction (or loss), tail lengthening (for balance), and pelvis expansion (to support enlarged hindlimbs) (Lebedev et al., 2013; Michaux & Shenbrot, 2017; Miljutin, 2008; Moore et al., 2015; Nowak & Paradiso, 1983; Shenbrot et al., 1999, 2008; Wu et al., 2014).

Allactagines are among the most specialized dipodids for bipedality, with associated skeletomuscular modifications such as unfused cervical vertebrae and highly modified hindlimbs and pelvises (Allen, 1940; Holden & Musser, 2005; Lebedev et al., 2013). The degree of hindlimb specialization is one of the main criteria for dipodid taxonomy and morphological phylogenetics, in addition to the morphology of teeth (e.g., molar coronal structure), skulls (particularly auditory bullae), skeletons (e.g., cervical vertebrae), and male reproductive system (glans penis and accessory glands) (see Nowak & Paradiso, 1983 and references therein). Recent molecular phylogenies support the monophyly of Dipodidae as a family, as well as the monophyly of its constituent subfamilies (Lebedev et al., 2013; Moore et al., 2015; Pisano et al., 2015; Upham et al., 2019; Wu et al., 2012; Zhang et al., 2013). Molecular phylogenetic studies have also greatly improved our understanding of the relationships among allactagine species (Figure 1 and references therein) and contributed to the recent large-scale revision of the group’s taxonomy (see Michaux & Shenbrot, 2017).

Figure 1

Allactagine relationships based on several molecular phylogenies published in the last decade (Dianat et al., 2013; Hamidi et al., 2016; Lebedev et al., 2013; Moore et al., 2015; Pisano et al., 2015; Upham et al., 2019; Wu et al., 2012; Zhang et al., 2013)

Relationships are shown as cladograms (with arbitrary branch lengths), with conspecific clades of individuals (based on authors of each study) collapsed into a single tip (see notes). For most studies, the allactagine clade was extracted from a larger phylogeny, and whenever available, the tree based on concatenation of all genes was used. Taxonomy of species names was updated following Michaux & Shenbrot (2017). Cladograms were constructed using R library ape (Paradis & Schliep, 2019).

Traditionally, certain allactagines are easy to recognize, e.g., Pygeretmus, which usually lacks premolars and has simple-patterned enamel, shorter ears, and fatter tails; others are differentiated based on less obvious differences such as fur color, body size, skull measurements, and glans penis spines (Allen, 1940; Kryštufek et al., 2013; Miljutin, 2008; Nowak & Paradiso, 1983; Shenbrot, 2009). This interspecific variation, particularly in morphological measurements, is a much-studied topic, particularly in the Iranian allactagines (e.g., Darvish et al., 2008; Dianat et al., 2010; Lebedev et al., 2013; Miljutin, 2008; Shenbrot, 2009; Tarahomi et al., 2010). In the present study, I sampled allactagine species throughout their extensive range from Libya to China (Figure 2) to assess interspecific cranial variation. Geometric morphometrics were used to compute and visualize size-independent shape variation (see Bookstein, 1986; Goodall, 1991), with the usefulness of this technique for allactagine species identification also investigated.

Figure 2

Map of localities sampled in this study

Genera are designated with different colors, and species are designated with different symbols; darker tones are more heavily sampled localities. The x- and y-axes are longitude and latitude, respectively (using geographic WGS84 CRS). Supplementary Data S1 lists geographic coordinates of localities. Locality of a single specimen (A. major, AMNH 227) was unknown and does not appear in this map. Map was generated using ggplot2 (Wickham, 2016), rnaturalearth (South, 2017), and ggspatial (Dunnington, 2020).

MATERIALS AND METHODS

Taxonomic sampling

The crania of vouchered specimens used in this study are housed in the American Museum of Natural History (AMNH; New York City), Field Museum of Natural History (FMNH; Chicago), Museum of Vertebrate Zoology (MVZ; University of California, Berkeley), and United States National Museum of Natural History (USNM; Washington, DC) (Appendix I; Supplementary Data S1). The samples excluded juveniles with incompletely erupted third mandibular molars (see Auffray et al., 2011; Shenbrot et al., 2008). Damaged specimens were also omitted to ensure shape estimate reliability (see below). All allactagine species available in these museums (n=14) were sampled: Allactaga major Kerr, 1792 (n=6), A. severtzovi Vinogradov, 1925 (n=3), Orientallactaga balikunica Hsia and Fang, 1964 (n=3), O. bullata Allen, 1925 (n=26), O. sibirica Forster, 1778 (n=36), Pygeretmus platyurus Lichtenstein, 1823 (n=1), P. pumilio Kerr, 1792 (n=9), P. shitkovi Kuznetsov, 1930 (n=2), Scarturus aralychensis Satunin, 1901 (n=5), S. euphraticus Thomas, 1881 (n=6), S. hotsoni Thomas, 1920 (n=25), S. indicus Gray, 1842 (n=32), S. tetradactylus Lichtenstein, 1823 (n=40), and S. williamsi Thomas, 1897 (n=25), totaling 219 specimens (Appendix I; Supplementary Data S1). Subspecific epithets were not considered in this study.

The taxonomy of the specimens used by the museum databases was primarily based on Holden & Musser (2005). In the current study, I updated this taxonomy (particularly genera) using Michaux & Shenbrot (2017). This reference was also used to re-identify several sampled specimens, particularly Scarturus elater Lichtenstein, 1825, as follows: S. elater sampled from Afghanistan and Pakistan was re-identified as S. indicus, while those sampled from central Iran (all other S. elater specimens) were re-identified as S. aralychensis. Recently, Bannikova et al. (2019) found evidence for several cryptic species within the S. elater species complex, and thus the S. indicus and S. aralychensis specimens identified in this study based on Michaux & Shenbrot (2017) may include other cryptic species. All sampled O. sibirica specimens were from localities that conformed to the “Sibirica” phylogenetic group of Cheng et al. (2020).

Other changes included re-identifying the O. bullata specimen (AMNH 46403) from Ulaanbaatar in Mongolia, which lacks the characteristic hyperinflated tympanic bulla of this species, as O. sibirica. The O. s. annulata specimen (AMNH 58832) with a hyperinflated tympanic bulla was re-identified as O. bullata. Finally, the A. major specimen from Turkmenistan (USNM 1445) was re-identified as A. severtzovi, as the former does not inhabit this country according to Michaux & Shenbrot (2017). Two S. firouzi specimens from FMNH (112349 and 112351) were included as S. hotsoni (Dianat et al., 2010, 2013; Shenbrot, 2009). All other synonymous names were standardized following Michaux & Shenbrot (2017). All revisions to species/genus names (including the original identifications provided by the museums) are listed in the Supplementary Data S1.

Specimens were assigned geographic coordinate values based on museum database locality information: AMNH (sci-web-001.amnh.org/db/emuwebamnh), FMNH (collections-zoology.fieldmuseum.org), MVZ (arctos.database.museum/mvz_mamm), and USNM (collections.nmnh.si.edu/search/mammals) (Supplementary Data S1). Several specimens were close to or over a century old and lacked precise locality information and/or geographic coordinates and were thus assigned latitudes and longitudes based on the most precise localities listed in the museum databases using the distance-measuring tool in Google Maps (Google, 2020), as described in Alhajeri (2019). Geographic coordinates were obtained in decimal degrees to the fourth decimal place (precision of ~11 m at the equator). The locality of a single specimen (A. major, AMNH 227) was unknown.

Geometric morphometric analysis

All analyses and visualizations used R libraries (R Core Team, 2020). For all analyses, the significance level (P-value) was set at α=0.05 and evaluated with permutations (with a random starting seed). Table legends note the specimens missing from each analysis.

For each specimen, the ventral cranial view (see Discussion) was photographed using a Nikon AF-S DX Micro NIKKOR 40 mm f/2.8G lens attached to a D3200 DSLR Nikon F-mount camera (Nikon, Tokyo, Japan) mounted on a copy stand. Specimens were illuminated using a ring light and external light sources, and photographs were captured at 24 megapixels with a resolution of 6,016×4,000 in JPEG format. The standardized photography protocol is described elsewhere (Alhajeri & Steppan, 2018a; Alhajeri, 2019) and ensured that the photographic plane was parallel to the frontal plane of the cranium. Photographs included 1 mm graph paper to convert pixels into millimeters and estimate centroid size (see below).

Thirty-two landmarks were digitized on the left side of each cranium, except when it was damaged; in these instances, the photograph was reflected, prior to digitizing the right side. The chosen landmarks are common in mammals and found in most rodents (Alhajeri & Steppan, 2018b), but also included allactagine-specific landmarks and evenly covered most of the ventral cranial view (Figure 3; Supplementary Table S1). Landmarks were digitized in ImageJ (Schneider et al., 2012) resulting in raw coordinate files, which were converted into tps format (Rohlf, 2015) using base R functions. These files were combined using the readmulti.tps() function in the geomorph library (Adams et al., 2020). ImageJ coordinates appear upside down in geomorph plots, and thus were flipped along the y-axis using the geomorph function rotate.coords (type=“flipY”).

Figure 3

Digitized landmarks on ventral cranial view of a small five-toed jerboa (S. indicus; FMNH 103883) from Herāt, Afghanistan (Landmarks are described in Supplementary Table S1)

Twenty-seven out of the 219 specimens were missing at least one landmark (12.3%), and a total of 52 landmarks were missing from the entire dataset (0.7%). The positions of these landmarks were imputed using thin-plate spline (TPS; Gunz et al., 2009) with the estimate.missing() function in geomorph. The references for missing landmarks from incomplete specimens were interpolated from complete specimens and carried out separately for each species. Specimens were then aligned using generalized Procrustes analysis (GPA) (Rohlf & Slice, 1990), a process of rescaling, translation, and rotation. The Procrustes-aligned specimens were projected onto tangent space using the geomorph function gpagen(Proj=T). This resulted in a three-dimensional array of Procrustes shape coordinates and a vector of centroid sizes—the square root of summed variances of each landmark around the centroid (Rohlf & Slice, 1990)—which were combined after converting the former into a two-dimensional matrix using the geomorph function two.d.array(). The resulting matrix was used in all analyses, as shown in Supplementary Data S1.

Morphometric data quality (e.g., absence of digitization errors) was ensured by the geomorph function plotOutliers() to detect and scrutinize outliers (with Procrustes distances above the upper quartile of species mean shape), which were then re-digitized or dropped (final 219 specimen dataset contained no outliers). In most instances, the dropped outlier specimens contained several missing landmarks that were imputed using TPS (see above). Centroid sizes were (natural) log-transformed and used as measures of cranial size, while Procrustes shape coordinates were the proxies of cranial shape.

Statistical analyses

The effect of sex on cranial size and shape and its interaction with species identity was tested by evaluating the fit of the models using the randomized residual permutation procedure (RRPP) via the lm.rrpp() function in the RRPP library (Collyer & Adams, 2018, 2020). These preliminary models found no significant sexual dimorphism in size or shape (see Results)—males, females, and those of unknown sex were combined in the final models. In these final models, RRPP was used to test the effects of: (1) species identity on cranial size (while not considering sex), and (2) cranial size, species identity, and their interaction on cranial shape (while not considering sex). All RRPP models were fit using ordinary least squares (OLS), and significance testing was based on 999 iterations of residual randomization, with type-II (hierarchical) sum of squares and cross-product computations used for multifactorial models. Analysis of variance (ANOVA) tables were computed for each model using anova.lm.rrpp() based on random statistical distributions for each model, with the F distribution as the basis for effect size calculation.

For cranial size, pairwise comparisons were conducted for the only significant factor, species identity (see Results), using the RRPP function pairwise(). For cranial shape, pairwise comparisons were conducted for species means, including and excluding cranial size as a covariate (if included it compared the least squares means of species, i.e., performed a correction for common allometry). The RRPP function summary.pairwise() was then used to compute the associated statistics, including distances among species pairs, effect sizes, P-values, and one-tailed 95% confidence limits, based on distances between means (Collyer & Adams, 2020). Differences in cranial size between species were visualized using boxplots.

The geomorph function gm.prcomp() was used for principal component analysis (PCA) of Procrustes shape variables to determine the main axes of shape variation. The first two PCs were visualized as scatterplots (with points color-coded by genera and species identity designated with different symbols) using the plot.gm.prcomp() function of the same library. Cranial shapes at PC axis extremes were visualized in geomorph by plotting their differences from the mean shape of the entire sample—estimated using the function mshape()—using TPS deformation grids generated via the function plotRefToTarget(); to aid visualization, the define.links() function was used to define links between landmarks, and deformations were magnified two-fold to exaggerate differences.

Comparable species-mean shapes were obtained in geomorph by subsetting the dataset by species, obtaining means using mshape(), and realigning resulting means via GPA using gpagen(). Species-mean shapes were depicted as TPS deformations compared to the mean shape of all 14 species means (thus giving equal weight to each species, regardless of sample size); grids were magnified two-fold to make shape differences more apparent. GPA also output a Procrustes distance matrix, which was summarized in a dendrogram based on UPGMA (unweighted pair group method with arithmetic mean) clustering implemented in the base R function hclust(method=“average”).

Because RRPP analysis found shape allometry, but no evidence that it differed among species (i.e., insignificant interaction between size and species identity, see Results), shape-size covariation was visualized using a simple (common) allometry scatterplot. This plot was based on the common allometric component (CAC; Mitteroecker et al., 2004) implemented using the function plotAllometry (method=“CAC”) in geomorph, which uses two-block partial least squares (PLS; Rohlf & Corti, 2000) to estimate the major axis of covariation between size and shape (i.e., CAC). Allometric shape changes were visualized using TPS deformation grids of minimum and maximum CAC scores (relative to the mean shape of the entire sample), using the same process outlined for PCA. Deformations were magnified four-fold to make shape differences more visible. The configurations of the shape coordinates at the CAC and PCA extremes were determined using the geomorph function shape.predictor().

RESULTS

Of the 219 specimens, 112 were females, 100 were males, and seven were of unknown sex—similar frequencies of males and females were sampled for most species (Supplementary Data S1). Specimens came from 11 countries, including Mongolia (n=43), Afghanistan (n=37), Egypt (n=37), Iran (n=33), China (n=19), Turkey (n=17), Kazakhstan (n=10), Pakistan (n=10), Russia (n=6), Libya (n=3), Turkmenistan (n=3), and a single specimen of unknown origin (Supplementary Data S1). A total of 103 localities were sampled (Figure 2; Supplementary Data S1).

Preliminary models found no significant sexual dimorphism in size or shape (effects of sex and its interaction with species identity were non-significant) (all P>0.05; Supplementary Table S2). The final cranial size model found significant and large differences among species (P=0.001; Z=10.658), with species identity explaining 87.5% of the variation in cranial size (R2=0.875; Table 1A). The boxplot showed some clustering among genera relative to their size, with both members of the genus Allactaga being the largest, followed (mostly) by Orientallactaga, then Scarturus, and finally Pygeretmus (Figure 4). Post-hoc pairwise comparisons (which excluded P. platyurus, see Table 2 legend notes) found significant differences in cranial size between almost all species pairs, with most insignificant pairs belonging to the same genus (Table 2A).

Table 1

ANOVA table for linear model of effect of species identity (‘Species’) on logged cranial size (A) and logged cranial size (‘Size’), species identity, and their interaction on cranial shape (Procrustes shape coordinates) (B)

	df	SS	MS	R2	F	Z	P
df: Degrees of freedom; SS: Sum of squares—hierarchical (type-II) was used for (multifactorial) Procrustes shape model; MS: Mean squares; R2: R-squared value; F: F-value; Z: Effect size (standard deviation of F sampling distribution); P: P-values based on 999 permutations. Fit of each linear model was evaluated using randomized residual permutation procedure (RRPP). Model coefficients were estimated using ordinary least squares (OLS). Significant model terms are in bold. All 219 specimens are included in these models.
(A) Centroid size (log)
Species	13	3.8117	0.2932	0.8754	110.88	10.658	0.001
Residuals	205	0.5421	0.0026	0.1245
Total	218	4.3537
(B) Procrustes shape
Size	1	0.0071	0.0071	0.0128	6.81	5.419	0.001
Species	13	0.2552	0.0196	0.4594	18.73	16.966	0.001
Size×species	12	0.0151	0.0012	0.0273	1.20	1.602	0.061
Residuals	192	0.2012	0.0010	0.3622
Total	218	0.5555

Figure 4

Boxplot of recorded cranial sizes

Table 2

Pairwise distances between species means and associated statistics based on logged centroid size (A), Procrustes shape coordinates (B), and size-corrected Procrustes shape coordinates (least squares means with logged centroid size as a covariate and no interaction) (C)

	(A) Size				(B) Shape (no correction)				(C) Shape (size correction)
	d	UCL	Z	P	d	UCL	Z	P	d	UCL	Z	P
d: Distance between (least squares) means; UCL: 95% upper confidence limit for d (one-tailed); Z: Effect size; P: P-values based on 999 permutations. Significantly different pairs are in bold. P. platyurus consists of a single specimen, and thus could not be included in these pairwise comparisons.
A. major — A. severtzovi	0.16	0.19	1.4	0.105	0.03	0.05	–0.3	0.593	0.03	0.05	–0.3	0.573
A. major — O. balikunica	0.46	0.20	6.2	0.001	0.10	0.05	7.1	0.001	0.08	0.06	3.3	0.008
A. major — O. bullata	0.35	0.13	7.3	0.001	0.09	0.03	11.1	0.001	0.08	0.04	5.9	0.001
A. major — O. sibirica	0.22	0.12	4.5	0.002	0.05	0.03	5.8	0.001	0.05	0.03	4.5	0.003
A. major — P. pumilio	0.63	0.15	11.4	0.001	0.05	0.04	4.0	0.002	0.06	0.07	1.3	0.105
A. major — P. shitkovi	0.54	0.23	6.3	0.001	0.05	0.06	1.1	0.133	0.08	0.07	2.4	0.029
A. major — S. aralychensis	0.56	0.17	8.7	0.001	0.08	0.04	6.4	0.001	0.07	0.06	2.2	0.031
A. major — S. euphraticus	0.40	0.16	6.6	0.001	0.07	0.04	4.7	0.001	0.06	0.05	2.6	0.015
A. major — S. hotsoni	0.48	0.13	10.3	0.001	0.10	0.03	12.0	0.001	0.08	0.05	4.6	0.002
A. major — S. indicus	0.51	0.13	11.2	0.001	0.07	0.03	7.6	0.001	0.06	0.06	2.2	0.030
A. major — S. tetradactylus	0.46	0.12	10.3	0.001	0.07	0.03	8.3	0.001	0.06	0.05	2.5	0.017
A. major — S. williamsi	0.35	0.12	7.3	0.001	0.06	0.03	5.8	0.001	0.05	0.04	2.5	0.020
A. severtzovi — O. balikunica	0.30	0.23	2.8	0.014	0.09	0.06	4.8	0.001	0.07	0.06	3.1	0.007
A. severtzovi — O. bullata	0.19	0.16	2.4	0.021	0.08	0.04	6.8	0.001	0.07	0.04	6.1	0.001
A. severtzovi — O. sibirica	0.06	0.16	–0.3	0.521	0.04	0.04	1.9	0.049	0.04	0.04	2.2	0.032
A. severtzovi — P. pumilio	0.46	0.18	6.5	0.001	0.05	0.05	2.7	0.013	0.06	0.06	1.9	0.047
A. severtzovi — P. shitkovi	0.38	0.24	3.7	0.002	0.06	0.06	1.0	0.145	0.07	0.07	2.4	0.024
A. severtzovi — S. aralychensis	0.39	0.20	5.0	0.001	0.07	0.05	3.6	0.003	0.06	0.06	2.1	0.039
A. severtzovi — S. euphraticus	0.24	0.19	2.6	0.019	0.06	0.05	2.4	0.021	0.05	0.05	2.0	0.038
A. severtzovi — S. hotsoni	0.32	0.16	4.9	0.001	0.08	0.04	7.1	0.001	0.07	0.05	4.2	0.002
A. severtzovi — S. indicus	0.35	0.16	5.5	0.001	0.06	0.04	4.1	0.001	0.05	0.05	2.4	0.032
A. severtzovi — S. tetradactylus	0.29	0.17	4.4	0.001	0.06	0.04	4.5	0.002	0.05	0.05	2.7	0.017
A. severtzovi — S. williamsi	0.19	0.17	2.3	0.029	0.05	0.04	2.2	0.027	0.04	0.04	1.5	0.077
O. balikunica — O. bullata	0.11	0.17	0.7	0.222	0.04	0.04	1.5	0.077	0.04	0.04	1.6	0.070
O. balikunica — O. sibirica	0.24	0.17	3.3	0.009	0.08	0.04	6.4	0.001	0.07	0.04	4.6	0.002
O. balikunica — P. pumilio	0.17	0.19	1.6	0.081	0.08	0.05	6.1	0.001	0.09	0.05	7.5	0.001
O. balikunica — P. shitkovi	0.08	0.27	–0.3	0.512	0.10	0.07	4.6	0.001	0.10	0.06	6.0	0.001
O. balikunica — S. aralychensis	0.09	0.22	0.2	0.358	0.06	0.05	3.1	0.006	0.07	0.05	4.0	0.004
O. balikunica — S. euphraticus	0.06	0.20	–0.3	0.536	0.06	0.05	3.2	0.006	0.06	0.05	3.5	0.004
O. balikunica — S. hotsoni	0.02	0.17	–0.9	0.822	0.04	0.04	0.8	0.189	0.04	0.04	1.3	0.098
O. balikunica — S. indicus	0.05	0.17	–0.4	0.558	0.06	0.04	4.0	0.002	0.06	0.04	4.9	0.002
O. balikunica — S. tetradactylus	0.00	0.17	–1.2	0.975	0.07	0.04	5.0	0.001	0.07	0.04	5.8	0.001
O. balikunica — S. williamsi	0.11	0.18	0.8	0.205	0.06	0.04	4.4	0.001	0.06	0.04	4.4	0.003
O. bullata — O. sibirica	0.13	0.07	4.7	0.001	0.08	0.02	16.6	0.001	0.07	0.02	14.8	0.001
O. bullata — P. pumilio	0.27	0.11	6.7	0.001	0.08	0.03	12.7	0.001	0.10	0.04	10.9	0.001
O. bullata — P. shitkovi	0.19	0.19	1.7	0.064	0.09	0.05	6.3	0.001	0.10	0.05	8.0	0.001
O. bullata — S. aralychensis	0.20	0.14	3.4	0.003	0.07	0.03	8.4	0.001	0.08	0.04	8.7	0.001
O. bullata — S. euphraticus	0.05	0.13	–0.1	0.480	0.07	0.03	8.0	0.001	0.07	0.03	9.5	0.001
O. bullata — S. hotsoni	0.13	0.08	3.9	0.001	0.05	0.02	10.6	0.001	0.05	0.02	9.8	0.001
O. bullata — S. indicus	0.16	0.07	5.3	0.001	0.07	0.02	14.1	0.001	0.07	0.02	13.8	0.001
O. bullata — S. tetradactylus	0.10	0.07	3.5	0.004	0.08	0.02	16.1	0.001	0.08	0.02	17.0	0.001
O. bullata — S. williamsi	0.00	0.08	–1.2	0.922	0.07	0.02	13.2	0.001	0.07	0.02	14.8	0.001
O. sibirica — P. pumilio	0.41	0.10	11.1	0.001	0.07	0.03	11.2	0.001	0.08	0.05	6.0	0.001
O. sibirica — P. shitkovi	0.32	0.19	4.1	0.001	0.08	0.05	4.7	0.001	0.09	0.05	5.6	0.001
O. sibirica — S. aralychensis	0.34	0.14	6.8	0.001	0.07	0.03	7.8	0.001	0.07	0.04	4.9	0.001
O. sibirica — S. euphraticus	0.18	0.12	3.5	0.004	0.05	0.03	5.3	0.001	0.05	0.03	4.8	0.002
O. sibirica — S. hotsoni	0.26	0.07	10.6	0.001	0.07	0.02	16.1	0.001	0.06	0.03	8.5	0.001
O. sibirica — S. indicus	0.29	0.07	11.9	0.001	0.06	0.02	13.6	0.001	0.06	0.03	6.4	0.001
O. sibirica — S. tetradactylus	0.24	0.06	10.7	0.001	0.07	0.02	16.3	0.001	0.06	0.03	9.7	0.001
O. sibirica — S. williamsi	0.13	0.07	4.6	0.001	0.05	0.02	10.2	0.001	0.04	0.02	8.4	0.001
P. pumilio — P. shitkovi	0.09	0.21	0.1	0.395	0.05	0.06	0.8	0.192	0.04	0.05	0.5	0.261
P. pumilio — S. aralychensis	0.07	0.15	0.2	0.363	0.05	0.04	3.6	0.002	0.05	0.04	4.6	0.002
P. pumilio — S. euphraticus	0.23	0.15	3.7	0.004	0.06	0.04	5.3	0.001	0.07	0.04	5.6	0.001
P. pumilio — S. hotsoni	0.15	0.11	3.1	0.007	0.08	0.03	11.4	0.001	0.08	0.03	12.4	0.001
P. pumilio — S. indicus	0.12	0.10	2.4	0.024	0.04	0.03	5.5	0.001	0.05	0.03	6.4	0.001
P. pumilio — S. tetradactylus	0.17	0.10	4.1	0.003	0.04	0.03	4.8	0.001	0.05	0.03	5.8	0.001
P. pumilio — S. williamsi	0.28	0.11	6.8	0.001	0.06	0.03	7.3	0.001	0.06	0.04	6.3	0.001
P. shitkovi — S. aralychensis	0.01	0.22	–1.1	0.903	0.07	0.06	2.9	0.012	0.07	0.05	3.4	0.004
P. shitkovi — S. euphraticus	0.14	0.21	0.8	0.200	0.06	0.06	2.4	0.025	0.07	0.06	3.6	0.002
P. shitkovi — S. hotsoni	0.06	0.19	–0.3	0.557	0.09	0.05	6.3	0.001	0.09	0.05	7.6	0.001
P. shitkovi — S. indicus	0.03	0.19	–0.8	0.760	0.06	0.05	3.4	0.006	0.07	0.05	4.1	0.002
P. shitkovi — S. tetradactylus	0.08	0.19	0.1	0.407	0.06	0.05	3.2	0.007	0.07	0.05	4.4	0.001
P. shitkovi — S. williamsi	0.19	0.20	1.8	0.058	0.06	0.05	2.9	0.010	0.07	0.05	4.3	0.002
S. aralychensis — S. euphraticus	0.15	0.17	1.6	0.080	0.04	0.04	1.2	0.124	0.04	0.04	1.3	0.107
S. aralychensis — S. hotsoni	0.07	0.13	0.5	0.289	0.04	0.04	3.5	0.001	0.05	0.03	4.5	0.002
S. aralychensis — S. indicus	0.04	0.13	–0.2	0.489	0.02	0.03	–0.2	0.518	0.02	0.03	0.0	0.432
S. aralychensis — S. tetradactylus	0.10	0.13	1.1	0.143	0.02	0.03	0.2	0.368	0.03	0.03	0.8	0.187
S. aralychensis — S. williamsi	0.20	0.13	3.6	0.003	0.04	0.04	2.7	0.010	0.04	0.04	2.4	0.020
S. euphraticus — S. hotsoni	0.08	0.12	0.7	0.224	0.05	0.03	4.9	0.001	0.05	0.03	4.9	0.002
S. euphraticus — S. indicus	0.11	0.13	1.5	0.092	0.03	0.03	2.2	0.040	0.03	0.03	2.4	0.024
S. euphraticus — S. tetradactylus	0.06	0.12	0.1	0.401	0.04	0.03	3.4	0.006	0.04	0.03	4.0	0.003
S. euphraticus — S. williamsi	0.05	0.12	0.0	0.441	0.02	0.03	–0.5	0.640	0.02	0.03	0.0	0.461
S. hotsoni — S. indicus	0.03	0.08	–0.1	0.471	0.04	0.02	9.0	0.001	0.05	0.02	10.4	0.001
S. hotsoni — S. tetradactylus	0.02	0.07	–0.2	0.513	0.05	0.02	11.5	0.001	0.05	0.02	12.6	0.001
S. hotsoni — S. williamsi	0.13	0.08	4.1	0.001	0.05	0.02	10.4	0.001	0.05	0.02	8.6	0.001
S. indicus — S. tetradactylus	0.05	0.06	1.4	0.115	0.03	0.02	4.9	0.001	0.03	0.02	5.8	0.001
S. indicus — S. williamsi	0.16	0.07	5.6	0.001	0.03	0.02	5.9	0.001	0.03	0.02	5.0	0.001
S. tetradactylus — S. williamsi	0.11	0.07	3.7	0.003	0.04	0.02	9.0	0.001	0.04	0.02	7.9	0.001

For each box, inner line is the median, hinges are the first and third quartiles, whiskers extend from each hinge to the maximum/minimum value (not exceeding 1.5 times the interquartile range), and outliers are plotted as points. Plot was generated using ggplot2.

The final cranial shape model also found significant and large differences between the species (P=0.001; Z=16.966), with species identity explaining 45.9% of the variation in cranial shape (R2=0.459; Table 1B). Cranial shape was also significantly predicted by its size (P=0.001; Z=5.419), with size explaining only 1.2% of the variation in cranial shape (R2=0.012)—size did not interact with species identity (P>0.05;Table 1B), indicating a common allometric pattern across species. Post-hoc pairwise comparisons found significant differences in cranial shape between most species pairs (Table 2B), even after size correction (Table 2C).

Genera and species clearly differed in the first two PC axes, which accounted for 50.1% of cranial shape variation (Figure 5). PC1 distinguished most O. bullata+O. balikunica+S. hotsoni specimens with low scores from most Allactaga (A. major+A. severtzovi)+O. sibirica specimens with high scores—most other species had intermediate PC1 scores (Figure 5). Most Pygeretmus specimens had positive, but intermediate PC1 scores (Figure 5). According to the TPS grids, the former group (PC1 min) was characterized by wider cranial bases, and especially anteriorly inflated bullae, and a reduced facial region, particularly a less flared (and anteriorly placed) zygomatic arch, a smaller and shorter rostrum, and a somewhat larger molar row. Converse patterns were observed in the latter group (PC1 max) (Figure 5). PC2 differentiated most Orientallactaga specimens (PC2 max), with caudolaterally expanded cranial bases (and in the case of O. bullata extremely inflated bullae), less flared zygomatic arches, reduced foramina magna, and enlarged rostra, from most Scarturus+Pygeretmus specimens (PC2 min), which showed the opposite patterns (Figure 5). Comparisons of species-mean shapes (Figure 6) agreed with the PCA (Figure 5). The UPGMA clusters were mostly consistent with the four genera; except for convergence of S. hotsoni to the O. bullata+O. balikunica cluster, and convergence of O. sibirica to the Allactaga (A. major+A. severtzovi) cluster (Figure 6).

Figure 5

Scatterplot of Procrustes shape coordinates based on PCA of 219 specimens

Figure 6

UPGMA dendrogram showing similarity in species-mean shapes (based on Procrustes shape coordinates) along with associated TPS deformation grid for each species (compared to average of all 14 species-mean shapes)

Axes show principal components (PC) 1–2 and percentage of total shape variation explained by each. Genera are designated with different colors, and species are designated with different symbols (following legend in Figure 2). Quadrants are demarcated by dotted lines passing through PC origins. TPS deformation grids show shapes at PC extremes (compared to sample mean). To aid visualization, landmarks (in same order and orientation as Figure 3) are linked, and deformations are magnified two-fold. Plot was generated using geomorph.

To aid visualization, landmarks (in same order and orientation as Figure 3) are linked, and deformations are magnified two-fold. Node ‘height’ is indicated below dendrogram. Figure was generated using base R functions and geomorph.

Except for a few specimens of Pygeretmus+Allactaga, most specimens fit closely to the common allometric regression vector in the CAC plot (Figure 7). Small crania (CAC min) had anteriorly inflated bullae, larger foramina magna, shorter snouts, wider zygomatic arches+molars+incisive foramina, and rostrocaudally compressed crania, with large crania (CAC max) primarily showing the reverse patterns (Figure 7).

Figure 7

Scatterplot of shape-size covariation of 219 specimens

Logged centroid sizes and common allometric component (CAC) are on x- and y-axes, respectively. Genera are designated with different colors, and species are designated with different symbols (following Figure 2). TPS deformation grids show shapes at CAC extremes (compared to sample mean). To aid visualization, landmarks (in same order and orientation as Figure 3) are linked, and deformations are magnified four-fold. Plot was generated using geomorph.

DISCUSSION

I found that allactagine species differed in cranial size (Figure 4) and shape (Figures 5, 6; Tables 1, 2). This is in agreement with recent geometric morphometric analysis of Iranian allactagines by Tarahomi et al. (2010), despite that earlier study examining the dorsal cranial view. The differences in the ventral cranial view in the present study are particularly noteworthy, considering the presence of ecologically pertinent components such as teeth and tympanic bullae, also traditionally used for rodent species delineation. Allactagines are commonly delineated based on hindlimbs, crania, teeth, skeletons, glans penis spines, fur color, body size, and/or similar traits (Allen, 1940; Kryštufek et al., 2013; Nowak & Paradiso, 1983; Shenbrot, 2009). This study showed that allactagine species can also be differentiated using ventral cranial geometric morphometrics, bolstering evidence that this technique contains sufficient resolution for species identification and for elucidating interspecific differences, particularly in rodents (Alhajeri, 2018; Beolchini & Corti, 2004; Boroni et al., 2017; Colangelo et al., 2010; Fadda & Corti, 2001; Tabatabaei Yazdi & Adriaens, 2013; Tarahomi et al., 2010; Zelditch et al., 2012).

Despite finding significant differences in cranial size among species overall (Table 1A) and between most species pairs (Table 2A), Figure 4 shows a highly continuous pattern of size variation among species, with some clustering among genera. This indicates that cranial size may have sufficient resolution to distinguish allactagine genera, but not species. Some size variation may be associated with latitude, and thus could be explained by Bergmann’s rule (1847). For example, the largest sampled species A. major (Figure 4) also had the northernmost range (Figure 2); similarly, the relatively large-sized Orientallactaga species also inhabited northern latitudes.

On the other hand, cranial shape seems to better distinguish both genera and species (Figures 5, 6; Tables 1A, 2B, 2C). While some shape variation does seem to be driven by size (Figure 7), the low effect size of this factor relative to species identity (comparing R2 and Z values in Table 1B) seems to indicate that allometry plays a relatively minor role in interspecific shape differences at the taxonomic scale of the Allactaginae subfamily. While both PCA and UPGMA analyses found large-scale clustering among congeners, interspecific differentiation (in PCA) was less obvious for some species (see below) (Figures 5, 6).

Based on UPGMA of the mean shapes of species, all Pygeretmus species were clustered together (as also reflected in the PCA plot), which is expected considering the distinctive crania of this genus, e.g., relatively extreme zygomatic flaring (Figures 5, 6). The uniqueness of Pygeretmus among allactagines extends to other skull characters (e.g., lacking premolars) and non-skull characters (e.g., fat tail), which is why it was considered as a distinct allactagine by early taxonomists. The most morphologically close genus to Pygeretmus was Scarturus (five of the six congeners clustered together). These two genera seem to cluster together due to shared cranial features, such as reduced cranial bases and rostra along with enlarged zygomatic arches and foramen magna (Figures 5, 6).

Within the Scarturus cluster (excluding S. hotsoni), there was substantial overlap among the species in morphospace (Figures 5, 6). This is in accordance with recent evidence of the presence of cryptic species in this genus, particularly within the S. elater species complex (see Bannikova et al., 2019). Interestingly, S. indicus and S. aralychensis, which were originally considered as subspecies of S. elater (see Michaux & Shenbrot, 2017), were clustered together (Figure 6) within this group and showed considerable overlap in morphospace (Figure 5). The close association between S. williamsi and S. euphraticus within the Scarturus cluster (Figure 6) and their overlap in morphospace (Figure 5) may reflect their close phylogenetic relationship (Figure 1; also see Michaux & Shenbrot, 2017).

The remaining species were divided into two clusters: (1) A. major+A. severtzovi+O. sibirica and (2) O. bullata+O. balikunica+S. hotsoni (Figure 6). For PCA, these two clusters were most differentiated along PC1, with the former cluster having more deflated bullae and more enlarged zygomatic arches and rostra (Figure 5). In the first group, the clustering of O. sibirica with A. major+A. severtzovi (rather than with its congeners) (Figures 5, 6) could be driven partly by allometry, as these three species were the largest among the sampled allactagine species (Figures 4), and increased size according to CAC analysis was mostly associated with reduced bullae and foramina magna as well as enlarged rostra (Figure 7). Within this cluster, A. major greatly differed from the rest (and was also among the most distinct of all sampled allactagines), which could partly be explained by its much larger size compared with other species (Figure 4) and because it diverged somewhat from the common allometric regression vector (Figure 7).

The O. bullata+O. balikunica+S. hotsoni cluster was the most distinct among all sampled allactagines. Furthermore, it clustered away from all other species (Figure 6) and occupied a largely different region in PCA morphospace (Figure 5). This cluster was mainly characterized by extremely inflated bullae (especially for the case of O. bullata+O. balikunica); however, there was a gradual increase in the relative size of both the bullae and rostra in the following order: S. hotsoni – O. balikunica – O. bullata (Figure 6), which may partly be driven by allometry as it was associated with species size (Figures 4, 7). The clustering of S. hotsoni with O. bullata+O. balikunica (rather than with its congeners) (Figures 5, 6) may partly be driven by convergent adaptations to similar habitats, which for these three species include rocky-gravel deserts (see Michaux & Shenbrot, 2017). The bulla inflation in this cluster of species is a hallmark of desert adaptation in mammals, particularly rodents (Alhajeri et al., 2015; Alhajeri & Steppan, 2018c; Lay, 1972; Webster & Webster, 1975). Interestingly, this cluster of rodents also shared other adaptations, including densely furred soles (characteristic of Orientallactaga, but also shared by other allactagines) (see Michaux & Shenbrot, 2017), which is presumably a convergent adaptation to life in flat sandy-gravel deserts. Among all sampled species, O. bullata had the most distinctive cranium (i.e., greatly separated from other species in PC morphospace), which was primarily driven by its greatly inflated tympanic bullae (Figure 5).

Overall, the cranial variation identified here (Figures 5, 6) accords with earlier studies (e.g., Allen, 1940; Michaux & Shenbrot, 2017; Miljutin, 2008). Rodent crania are highly adaptable and readily evolve in response to selective pressure, such as diet (Martin et al., 2016). This is especially the case for allactagines, with species showing elaborate cranial specializations for traits such as digging (e.g., developed nasal disk and proodont incisors) and mastication (associated with zygomatic breadth and molar row length) (Miljutin, 2008). The cranial differences across allactagines likely reflect phylogenetic history (i.e., phylogenetic conservatism), adaptation (e.g., to local climate, diet, habitat), and random genetic drift. As close relatives usually also share adaptations (i.e., phylogenetic niche conservatism), which explains the strong congeneric clustering detected in this study (Figures 5, 6), cranial morphology likely manifests as an interplay of several of these factors.

CONCLUSIONS

The taxonomy of the subfamily Allactagine has been revised extensively in recent years with the aid of molecular phylogenetics. This study sampled 14 allactagine species and found significant differences in their cranial size and shape, indicating that geometric morphometrics may have sufficient resolution for allactagine species identification. Most congeners clustered together, reflecting phylogenetic conservatism in cranial shape. The two exceptions (S. hotsoni clustering with O. bullata+O. balikunica and O. sibirica clustering with A. major+A. severtzovi) were likely influenced by convergent cranial adaptations and allometry, respectively. Cranial variation may also be explained by other factors, such as random genetic drift.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The author declares that he has no competing interests.

AUTHOR’S CONTRIBUTIONS

B.H.A. conceived and designed this study. B.H.A. analyzed the data and wrote the manuscript. The author read and approved the final version of the manuscript.

ACKNOWLEDGEMENTS

I thank the museum curators, collection managers, and staff for providing access to their collections and assistance on multiple visits: AMNH (Ms. Marisa Surovy, Ms. Eleanor Hoeger, Ms. Eileen Westwig); FMNH (Dr. Bruce Patterson, Dr. Lawrence Heaney, Dr. Adam Ferguson, Mr. John Phelps, Mr. William Stanley, Ms. Lauren Smith); MVZ (Dr. Christopher Conroy, Dr. James Patton, Dr. Eileen Lacey); and USNM (Mr. Darrin Lunde and Dr. Michael Carleton). I also thank Dr. Miriam Zelditch for discussions on R code and analyses. The manuscript benefitted from comments by Dr. Julien Claude and five anonymous reviewers. Some travel expenses were covered by Kuwait University (2019–2020 scientific leave).