Geographical Variation in the Sand Cat, Felis margarita (Carnivora: Felidae)

Abstract

Sand cats, Felis margarita, range from northern Africa and the Arabian Peninsula to Central Asia. Their apparently discontinuous distribution is recognized as comprising four subspecies. Recent genetic research found little differentiation between subspecies except for the North African form. In this study, 90 skins and 88 skulls were analyzed from the four subspecies. A discriminant function analysis of the scores, ranging from 1 to 4, of four pelage characteristics revealed differentiation between putative subspecies, except between Turkmenian and Pakistani sand cats. Northern African and Arabian sand cats tend to be spotted and striped, while Turkmenian and Pakistani sand cats are less spotted and have a dorsal crest of fur. Nonmetric multidimensional scaling (NMDS) models generated from 21 skull measurements revealed an overlap in morphospace between all subspecies, except for larger Turkmenian sand cats; northern African sand cats were smallest. Therefore, both pelage characteristics and skull morphometrics support up to three subspecies. However, considering recent genetic research, it is likely that two subspecies should be recognized, F. m. margarita from northern Africa and F. m. thinobia from the Arabian Peninsula, and Southwest and Central Asia. Widening of the dataset and nuclear DNA evidence are required to increase our understanding of geographical variation in this little studied species.

1. Introduction

The sand cat, Felis margarita, has a wide distribution in deserts and semi-deserts, ranging from Morocco and Algeria in the west through the Sahara Desert and the Arabian Peninsula to Southwest and Central Asia, including Israel, Syria, Jordan, Kuwait, Pakistan, Iraq, Iran, Afghanistan, Uzbekistan, Kazakhstan, and Turkmenistan (Figure 1) [1,2]. However, this distribution does not appear to be currently fully contiguous, such that different apparently fragmented populations are currently recognized as up to four subspecies, including F. m. margarita in northern Africa, F. m. harrisoni in the Arabian Peninsula and Middle East, F. m. scheffeli in the Nushki Desert of Pakistan, and F. m. thinobia in Southwest and Central Asia [3,4,5]. However, sand cats are recorded from Iran, Iraq, and Afghanistan [6,7,8,9,10], and it is unclear under this taxonomic arrangement to which subspecies the sand cats from these countries should be assigned.

1.1. Taxonomic History of the Sand Cat

Felis margarita was described by Loche [12] based on a mounted skin of an animal from Algeria, which was in a private collection in Paris. Sadly, this specimen no longer seems to survive, but it was clearly illustrated in Loche’s original description. Ognev [13] described a new genus and species of desert cat, Eremaelurus thinobius, from Repetek in the southeastern Karakum Desert, Turkmenistan, based on several characteristics, including a shorter and rounder skull than in Felis; they also noted that the vertical diameter of the auditory bulla exceeds its transverse diameter, the lower anterior edge of the orbit is thin, and the pelage lacks spots and stripes except for tail bands. Heptner and Dementiev [14] regarded Eremaelurus as a junior synonym of Felis and thinobius as being the same species as and probably a subspecies of F. margarita. Later, Pocock [15] followed Heptner and Dementiev [14] and recognized F. thinobius (sic.) as a subspecies of margarita.

Pocock [16] described a new subspecies, F. m. meinertzhageni, from El Golea, latitude 30° N, in the Algerian Sahara based on pelage differences (a larger black area on the back of the ear and the middle part of ear is olivaceous grey), but regarded Felis thinobia as a distinct species, owing to the large geographical gap between the two known taxa. Pocock [17] also described the subspecies F. m. airensis based on pelage differences of a specimen from In-Abbangarit, west of Air, French Sudan (Niger).

Following its initial discovery by J. Anderson (pers. comm.) and Lay et al. [18], Hemmer [3] described a new subspecies, Felis margarita scheffeli, based on specimens from the Nushki Desert, Pakistan, which had been brought into captivity by Walter Scheffel. Finally, Hemmer et al. [4] described Felis margarita harrisoni from the Arabian Peninsula based on a specimen from the northern edge of Umm as Samin, Oman (21°55′ N, 55°50′ E).

In a review, Schauenberg [5] recognized three subspecies (margarita, thinobia, and scheffeli), synonymizing meinertzhageni and airensis with the nominate margarita, while Hemmer et al. [4], after describing the subspecies harrisoni, recognized the four currently recognized subspecies based on the apparent separation between populations and differences in body size, pelage coloration, and skull proportions. Wozencraft [19] partly followed Hemmer et al. [4], but also listed airensis as a fifth subspecies. More recently, Sliwa [20] also followed Hemmer et al. [4] in recognizing the same four subspecies. Kitchener et al. [21] reviewed the subspecies taxonomy of Felis margarita and proposed two subspecies, F. m. margarita from northern Africa and F. m. thinobia from Southwest and Central Asia, including the Arabian Peninsula, based on pelage characteristics and the preliminary results of a genetic study [22]. Sunquist and Sunquist [23] also recognized these two subspecies.

1.2. Subspecies Recognition

Subspecies recognition has been based on size, pelage coloration and markings, and skull characteristics and measurements. Based on diagnoses from Hemmer et al. [4] and Schauenberg [5], the four currently recognized subspecies of Felis margarita are distinguished as follows:

F. m. margarita: small size; relatively narrow skull with relatively small auditory bulla tail; small carnassials, low narrow occiput; coloration relatively bright, well-marked with stripes and spots, with buffy-white paws and a buffy collar on throat; 2–6 tail bands.

F. m. thinobia: large size; relatively broad skull with low narrow occiput and relatively small auditory bullae; large carnassials; coloration darker and greyer than margarita with reduced spots and stripes; 2–3 tail bands.

F. m. scheffeli: males large, females small; skull broad with very large auditory bullae, but carnassials smaller than those in thinobia; occiput not expanded; coloration usually similar to thinobia, but more strongly striped and spotted individuals more frequent; more than eight tail bands (at least in kittens).

F. m. harrisoni: small size; broad skull with large auditory bullae and high broad occiput and large carnassials; similar bright coloration to margarita, but even more clearly spotted and striped; less extensive and diffuse darkening dorsally; clean white paws; 5–7 tail bands; ear patch smaller and less dark.

Overall, there appears to be a stepped cline in size from small in the west (Sahara) to large in the east (Turkmenistan).

Howard-McCombe et al. [22] have recently carried out a study on genetic variation between the four putative subspecies based on mitochondrial DNA, comprising the control region, NADH subunit 5, and cytochrome b genes (totaling 643 bp), from 47 animals of known geographical origin from across the sand cat’s range collected over the past 100 years. This study found little genetic evidence for the separation of the subspecies harrisoni, scheffeli, and thinobia from each other, but northern African sand cats, subspecies margarita, were more genetically distinct. However, there has been no comprehensive morphological assessment of geographical variation in sand cats throughout their range since the 1970s. In this study, we have re-examined the morphological basis for the recognition of subspecies in F. margarita, including analyses of pelage characteristics and skull measurements, and we consider this intraspecific variation in relation to the results of the recent molecular study.

2. Material and Methods

Sand cats are rare in museum collections, but we were able to measure skulls and assess pelages of all putative subspecies (

Table 1. The numbers of skins and skulls of putative subspecies of F. margarita used in this study.

Subspecies	Skins	Skulls
margarita	10	9
harrisoni	12	13
scheffeli	61	38
thinobia	16	28
Totals	99	88

and Table S1) from seven museum collections, i.e., the Field Museum of Natural History, Chicago (FMNH), the Muséum National d’Histoire Naturelle, Paris (MNHN), National Museums Scotland, Edinburgh (NMS), the Natural History Museum, London (NHMUK), the Senckenberg Forschungsinstitut und Naturmuseum, Frankfurt-am-Main (SMF), the Zoological Institute of the Russian Academy of Sciences, Saint Petersburg (ZIN), and the Zoological Museum of Moscow State University (ZMMU). Studied material included holotype specimens of airensis, harrisoni, meinertzhageni, scheffeli, and thinobia. The holotype of margarita is no longer extant.

2.1. Pelage Characters and Coloration

The study was based on 99 skins of adult F. margarita from localities across the species’ distribution range (Figure 1, Table S1). The skins originated from northwestern Africa and the Sahara (margarita, n = 10), Southwest and Central Asia (thinobia, n = 7), northeastern Africa and Arabia (harrisoni, n = 12), and Pakistan (scheffeli, n = 61).

Pelage characters and coloration were assessed from direct examination of museum skins and photographs taken during their examination. The following characters were recorded and scored on a range of 1 to 4 (Table S2):

• Presence of dark spots and stripes on flanks;
• Number of tail bands;
• Presence of sagittal crest of fur along dorsum;
• Size of black ear patch.

Scores for character state for each of these characters are recorded in

Table 2. Scores for character states for each pelage character (see Supplementary Materials).

Character	Scores for Character States
	1	2	3	4
Presence of dark spots and stripes on flanks	No spots and stripes	Spots and stripes on legs	Faint spots and stripes on body	Clear spots and stripes on body
Number of tail bands including tail tip	1–2 bands	3–4 bands	5–6 bands	-
Presence of dorsal crest of fur along dorsum	No dorsal crest	Weakly differentiated dorsal crest	Clear dorsal crest	-
Size of black ear patch	Small ear patch	-	Large ear patch	-

. Photographs of these character states are in the Table S2. Originally we scored ground coloration, but we decided to exclude this from the analyses, because we were concerned that we would not be able to give each specimen a reliable score based on images photographed under vastly different lighting conditions. Also, coloration differs between summer and winter pelages, and collecting dates were lacking for many of the specimens, especially the scheffeli specimens. Recently, it has been found that sand cats may rest in trees and, hence, their fur is potentially exposed to high levels of sunlight, which may cause fading of coloration [24]. Because pelage characters were scored as categorical/ordinal variables, Kruskal–Wallis tests were used to test if there were statistically significant differences in each pelage character amongst the four putative subspecies, and a discriminant analysis was performed to see if the four pelage characters distinguish putative subspecies (both tests were performed using SPSS 27, IBM, Armonk, NY, USA).

2.2. Skull Morphometrics

The study was based on 88 skulls of adult F. margarita from localities across the species’ distribution range (Figure 1, Table S1). The skulls originated from northwestern Africa and the Sahara (margarita, n = 9), Southwest and Central Asia (thinobia, n = 28), northeastern Africa and Arabia (harrisoni, n = 13), and Pakistan (scheffeli, n = 38).

Only adults were used in the analysis. Age classes were defined by scoring the morphological features of skull structure, such as the development of crests, the obliteration of sutures, e.g., basi-sphenoid suture, and a fully erupted adult dentition.

Twenty-one cranial measurements were made using digital sliding calipers to the nearest 0.1 mm (Figure S1, Supplementary information). To estimate missing data (due to damaged skulls), we used expectation-maximization [25], regression substitution [26], and multiple imputation algorithms [27]. In all cases, the hypothesis of a random distribution of missing values was accepted. In our case, the regression substitution method showed the minimum deviations of the sample statistics (mean, standard deviation) from their initial values. The estimates for missing values were further used to estimate morphological distances between any pair of individuals in the multivariate analysis. All univariate statistics and the sexual size dimorphism index were calculated using only actual measurements.

We used a discriminant analysis to determine the sex of unknown and problem specimens. There were 62 skulls (~70% of the sample) of known sex (26 females and 36 males). The discriminant function was developed on 52 of these specimens (26 females and 26 males) and 10 male skulls were used as a control group. As a result, 88.9% of males and 91.3% of females from the learning sample and 100% of males from the control sample were correctly identified by the best linear combination of six skull measurements (maxillary toothrow length, condylobasal length, occipital height, auditory bulla length, mandible length, and breadth at canine alveoli). We used this function to determine the sex of individuals whose sex was unknown.

As a measure of sexual size dimorphism in the i-th measurement (SSD_i), we chose the ratio 100 × (M_male − M_female)/(M_male + M_female), where M is the sample mean for the males or females. Average sexual size dimorphism (ASSD) was calculated as $\left(\sum _{i=1}^{21}{SSD}_i\right)/21$.

The basis of our approach to using multivariate analysis techniques [28] to describe morphometric variability is described in [29,30,31]. In order to exclude any influence of the “scale” of the different measurements on the results, all variables (v_i) were standardized as follows: $(v_i-v_{min})/(v_{max}-v_{min})$. As part of the multivariate analysis, we introduced a morphological or morphometric space (morphospace), which was constructed to provide a compact representation of the morphological distances between individuals in our sample. The geometry and the properties of the morphospace are determined both by the metric and by the method of dimension reduction on the original data [32]. We used two metrics for morphological distance estimation, namely the Euclidean distance and Kendall’s tau-b [33]. The first metric describes the dissimilarity in the “size” of the individuals, and the second, the rank correlation metric, describes the similarity in proportions (“shape”) between them [29].

Matrices of Euclidean (E) distance and Kendall (K) coefficients were used in non-metric multidimensional scaling (NMDS), which can be used as a non-parametric analogue of the principal components method [28,34,35,36,37]. This resulted in two morphospaces describing the variability in size (size variability model, SZM) and shape (shape variability model, SHM) of individuals. The coordinates of the SZM and SHM are further denoted as E1, E2, … and K1, K2…, respectively.

The “best minimum” dimensions for SZM and SHM were estimated based on a variation of the criterion known as “stress formula 1” or “Kruskal’s stress” [35], calculated for morphospaces with 1, 2, …, 15 coordinates [29,30,31,38,39,40]. We used an approach somewhat similar to the one used in case of the “scree test” in a PCA [41]. It is assumed that if the distribution of dissimilarities/similarities in the input matrix is close to random, the value of Kruskal’s stress must be maximal, and also monotonically and smoothly decreasing as the number of coordinates increases. The best NMDS solution is associated with a local Kruskal’s stress minimum and maximum deviation from the simulated stress value for a random dataset (we used a 100 × 100 random matrix).

We used components of variance analysis [42] on coordinates of SZM and SHM for estimation of an a priori (“putative”) taxonomic effect on their variances. The statistical significance of the ANOVA results was assessed using the F criterion. As we worked with small samples, we used a non-parametric Kruskal–Wallis test, also known as Kruskal–Wallis ANOVA, to compare the means of size coordinates in different putative subspecies [43]. Also, we used Mood’s median test [44], which is useful for comparing two or more groups of data, to determine if they come from populations with equivalent medians.

Statistical analyses of skull variables were performed using SPSS 25 (IBM, Armonk, NY, USA), NCSS 2012 [45] and PAST 4.02 [46].

3. Results

3.1. Pelage Variation

There are statistically significant differences in the presence of dark spots and stripes on flanks (spots) (Kruskal–Wallis test: df = 3, H = 16.27, p = 0.001), and in the presence of a crest of fur along the dorsum (dorsal crest) (Kruskal–Wallis test: df = 3, H = 30.26, p < 0.001). Putative subspecies harrisoni and margarita had clearer spots and stripes in comparison to scheffeli and thinobia, whilst the latter two had clearer dorsal crests than the other two (

Table 3. Sample sizes of putative subspecies and median scores for two pelage characters recorded from skins of sand cats.

Pelage Character	Subspecies	Sample Size	Median Score	Range
Spots	margarita	10	3	1–4
	harrisoni	12	3	1–4
	scheffeli	57	2	1–4
	thinobia	7	2	1–2
Dorsal crest	margarita	10	1	1–3
	harrisoni	11	2	1–2
	scheffeli	57	3	1–3
	thinobia	7	2	1–3

).

Three canonical discriminant functions were extracted by the discriminant analysis, where discriminant functions 1, 2, and 3 explain 86.4%, 8.8%, and 4.8% of the total variance, respectively.

Table 4. Classification results (number of specimens (% of specimens)) by the discriminant analysis of pelage characters; 77.8% of the specimens were classified correctly.

Actual Group	Predicted Scheffeli	Predicted Margarita	Predicted Harrisoni	Predicted Thinobia	Total
scheffeli	28 (90.3)	1 (3.2)	2 (6.5)	0 (0)	31 (100)
margarita	1 (14.3)	5 (71.4)	1 (14.3)	0 (0)	7 (100)
harrisoni	1 (10.0)	1 (10.0)	8 (80.0)	0 (0)	10 (100)
thinobia	5 (83.3)	0 (0)	0 (0)	1 (16.7)	6 (100)

and Figure 2 summarize the results. The discriminant analysis distinguished the four subspecies well in general (Wilks’ Lambda = 0.412, χ² = 43.46, df = 12, p < 0.001). It appears that putative subspecies scheffeli, margarita, and harrisoni are well distinguished from each other with 85.4% of specimens correctly classified, whilst thinobia cannot be distinguished from scheffeli (Table 4).

3.2. Skull Morphometrics

3.2.1. Sexual Size Dimorphism

The ASSD Index (

Table 5. Sexual size dimorphism index (SSD) in F. margarita cranial measurements. n—number of specimens, M ± m—means and error of means, M–W—Mann–Whitney U test Z (absolute values), p—statistical significance, R—average relative difference between males and females in %. n.s.—not significant.

Measure	Males		Females		M–W	p	R	SSDi
	n	M/±m	n	M/±m
Greatest length skull	50	91.5/0.72	34	86.0/0.76	4.73	<0.001	6.0	3.1
Condylobasal length	49	85.9/0.66	35	80.5/0.66	4.99	<0.001	6.3	3.2
Zygomatic width	51	70.2/0.68	35	65.1/0.71	4.49	<0.001	7.2	3.8
Mastoid width	51	43.9/0.26	35	42.0/0.32	4.01	<0.001	4.3	2.2
Neurocranium width	52	44.4/0.19	36	43.3/0.25	3.67	<0.001	2.6	1.3
Occipital height	49	24.3/0.15	33	22.9/0.15	5.55	<0.001	6.1	3.1
Auditory bulla length	50	25.6/0.13	35	24.3/0.16	5.14	<0.001	4.8	2.5
Auditory bulla width	50	16.70.12	35	16.1/0.14	3.10	0.002	3.7	1.9
Auditory bulla height	50	17.9/0.14	34	17.1/0.17	3.49	<0.001	4.4	2.3
Postorbital breadth	52	33.8/0.17	35	34.1/0.18	1.13	n.s.	-	-
Interorbital breadth	52	19.0/0.23	36	17.6/0.23	4.05	<0.001	7.7	4.0
Breadth at canine alveoli	52	22.6/0.20	35	20.80.25	5.06	<0.001	8.2	4.3
P4 length	49	11.0/0.07	34	10.4/0.10	4.14	<0.001	5.3	2.7
P4 width	49	4.70/0.06	34	4.4/0.06	3.59	<0.001	6.4	3.3
Maxillary toothrow length	49	28.5/0.16	34	26.5/0.19	6.10	<0.001	6.9	3.6
Nasal length	52	27.0/0.21	34	25.5/0.31	3.77	<0.001	5.5	2.8
Mandible length	52	61.2/0.49	36	56.9/0.50	5.16	<0.001	7.0	3.6
Mandible height	52	27.8/0.38	36	25.2/0.38	4.22	<0.001	9.3	4.9
M1 length	46	8.1/0.06	34	7.5/0.07	4.61	<0.001	6.7	3.4
Mandibular toothrow length	46	30.7/0.23	34	28.7/0.24	5.26	<0.001	6.5	3.4
Mandible depth at M1	52	10.2/0.13	34	9.2/0.12	5.00	<0.001	10.1	5.3
ASSD							6.2	3.1

and Tables S3–S6) in sand cat skulls is 3.1. On average, males are larger than females in all cranial variables by 3.7–10.1% except postorbital width, which shows no sexual dimorphism. The SSD index varied from 1.3 (neurocranium width) to 5.3 (depth of the mandible at M₁). Males and females differed most significantly in the lengths of the maxillary and mandibular toothrows, occipital height, auditory bulla length, and breadth at canine alveoli.

3.2.2. Multivariate Analysis

The morphological variability in the skulls of females and males was examined separately, in order to eliminate the influence of the effects of sexual dimorphism (Table 5). The number of morphospace coordinates was three and four in SZM for males and females, respectively, and three in SHM for both sexes. The linear combination of SZM and SHM coordinates accounted for 53–96% of the variance of the original cranial measurements (r²,

Table 6. Spearman’s rank-order correlations between coordinates of the SZM and SHM morphological spaces (E1–E4 and K1–K3) and measurements of F. margarita skull, r²—squared multiple correlation coefficients from multiple linear regression model. n.s.—not significant.

Measure	SZM Space				SHM Space			r2
	E1	E2	E3	E4	K1	K2	K3
	males
Greatest length skull	0.95	−0.01	−0.06		0.32	0.08	0.06	0.93
Condylobasal length	0.96	−0.08	−0.03		0.38	0.08	0.09	0.95
Zygomatic width	0.90	−0.28	0.08		0.59	0.18	0.17	0.91
Mastoid width	0.86	0.13	−0.01		0.15	0.21	0.08	0.85
Neurocranium width	0.73	0.46	−0.28		−0.13	0.34	−0.06	0.88
Occipital height	0.71	−0.22	−0.22		0.46	0.13	−0.09	0.76
Auditory bulla length	0.72	−0.04	0.09		0.18	−0.20	0.12	0.70
Auditory bulla width	0.57	0.24	0.31		−0.07	0.09	0.23	0.59
Auditory bulla height	0.75	0.16	0.37		0.15	0.39	0.43	0.82
Postorbital breadth	0.22	0.70	−0.51		−0.49	0.43	−0.31	0.80
Interorbital breadth	0.91	−0.03	0.14		0.41	0.27	0.27	0.89
Breadth at canine alveoli	0.83	−0.14	0.26		0.40	−0.01	0.34	0.79
P4 length	0.65	0.41	0.22		−0.16	0.21	0.28	0.69
P4 width	0.28	0.42	0.38		−0.28	−0.09	0.52	0.58
Maxillary toothrow length	0.86	−0.02	0.01		0.31	0.11	0.16	0.73
Nasal length	0.50	−0.51	−0.09		0.57	−0.08	−0.08	0.62
Mandible length	0.92	−0.29	−0.01		0.57	0.00	0.07	0.95
Mandible height	0.85	−0.25	0.00		0.57	0.25	0.09	0.88
M1 length	0.72	0.41	0.05		−0.14	−0.04	0.13	0.81
Mandibular toothrow length	0.77	−0.07	0.20		0.29	−0.25	0.34	0.79
Mandible depth at M1	0.83	−0.09	0.16		0.44	0.16	0.30	0.84
	Relative variance component, %
“Geographical factor”	48.0	34.7	34.7		15.4	35.1	30.8
	females
Greatest length skull	0.95	−0.01	−0.02	0.10	0.17	0.20	−0.31	0.92
Condylobasal length	0.94	0.02	0.06	0.10	0.09	0.21	−0.28	0.96
Zygomatic width	0.86	0.20	−0.13	−0.22	0.21	0.56	−0.11	0.90
Mastoid width	0.84	−0.16	−0.32	−0.07	0.50	0.32	−0.25	0.91
Neurocranium width	0.85	−0.28	−0.21	−0.12	0.41	0.06	−0.04	0.83
Occipital height	0.76	−0.07	−0.19	0.24	0.31	0.17	−0.58	0.81
Auditory bulla length	0.75	0.06	0.25	0.24	−0.11	0.12	−0.25	0.74
Auditory bulla width	0.75	−0.28	−0.01	−0.44	0.45	0.24	0.17	0.88
Auditory bulla height	0.72	−0.44	0.09	−0.31	0.46	0.09	0.18	0.88
Postorbital breadth	0.51	−0.50	−0.33	0.36	0.32	−0.46	−0.20	0.86
Interorbital breadth	0.86	0.13	0.04	−0.32	0.20	0.46	0.16	0.85
Breadth at canine alveoli	0.85	0.28	0.15	−0.22	0.04	0.54	−0.03	0.92
P4 length	0.76	−0.36	0.46	−0.01	0.29	0.03	0.05	0.88
P4 width	0.50	0.15	0.15	0.07	−0.19	−0.06	0.09	0.53
Maxillary toothrow length	0.87	0.06	0.24	0.10	0.00	0.18	−0.17	0.87
Nasal length	0.44	0.62	−0.06	0.19	−0.46	0.37	−0.14	0.75
Mandible length	0.84	0.32	−0.07	−0.15	0.00	0.53	−0.15	0.93
Mandible height	0.81	−0.01	−0.21	−0.04	0.28	0.21	−0.09	0.77
M1 length	0.74	−0.34	0.33	0.24	0.25	−0.06	−0.17	0.88
Mandibular toothrow length	0.75	0.34	0.13	−0.02	−0.15	0.43	−0.17	0.83
Mandible depth at M1	0.85	−0.14	−0.07	0.18	0.18	0.08	−0.37	0.84
	Relative variance component, %
“Geographical factor”	54.4	40.1	n.s.	31.0	14.8	29.4	12.2

).

In the SZMs, the first coordinates (E1) correlate with the “general size” of the skull. For males, the correlation coefficients between E1 and condylobasal length and zygomatic width were 0.96 and 0.90, respectively, and in females, they were 0.94 and 0.86, respectively. Only postorbital width and P⁴ width of males varied independently of total skull size (Table 6). Several measurements showed partial independence of variability from general cranial variations (Table 6), namely postorbital breadth, length of nasal bones, and lengths of P⁴ and M₁ in males, and postorbital breadth, the length of nasal bones, and the length of P⁴ in females.

In the SHMs, the first coordinate (K1), as well as the others, are not strongly correlated with the skull measurements (Table 6). For some measurements, a significant allometric effect (the relationship between size and shape) can be assumed, including zygomatic width, nasal length, mandible length/height, occipital height, postorbital breadth, and P⁴ width in males and the zygomatic width, nasal length, mastoid width, postorbital breadth, interorbital breadth, mandible length, and breadth at canine alveoli in females. Thus, allometric patterns differ between males and females, which should be further investigated in a future study.

The effects of putative intraspecific taxonomy (the random effect in the ANOVA design) on the variance of the coordinates varied across the width range (Table 6). According to these results, the skulls of the putative subspecies do not differ very significantly. The most significant test was obtained for the E1 coordinate in males and females. This means that the greatest differences between geographical forms of the sand cat are likely to be expressed in terms of their “general” skull size (

Table 7. Pearson correlations between coordinates of the SZM and SHM morphological spaces (E1–E4 and K1–K3).

	Males			Females
	K1	K2	K3	K1	K2	K3
E1	0.32	0.21	0.22	0.16	0.26	−0.13
E2	−0.87	0.18	0.14	0.61	−0.69	0.09
E3	0.01	−0.25	0.76	−0.34	0.01	0.31
E4				−0.33	−0.47	−0.56

). Along coordinate E1, the centroid of the thinobia sample differs from the centroids of the margarita–scheffeli samples, and the centroid of the harrisoni sample differs from the thinobia centroid along coordinate E2. In males, the centroids of margarita and scheffeli are distant from the centroid of harrisoni. A similar pattern was observed in females, but with less statistical significance. In general, there is a strong overlap between the ranges of points labelled scheffeli, harrisoni, and margarita (Figure 3A,B). Nevertheless, the skulls of thinobia are on average larger than those of other putative subspecies. The dendrograms in Figure 3C,D illustrate the separation of the sample centroids, taking into account the coordinates of the SZM and SHM, for which statistically significant ANOVA tests were obtained (Table 6). The thinobia sample centroid is separate from the others. The centroid of the harrisoni sample occupies an intermediate position between margarita, scheffeli and thinobia, but it is not far from the scheffeli centroid.

3.2.3. Univariate Analysis

At the level of individual measurements, the most significant results, corresponding to the most pronounced geographical patterns, were obtained for auditory bulla height, P⁴ length, and interorbital width in both sexes (

Table 8. The Kruskal–Wallis test (χ² (df = 3), p) and means of skull measurements by the geographical samples of F. margarita males and females.

Measure	χ2	p	Mean ± Standard Error
			Margarita	Scheffeli	Harrisoni	Thinobia
	males
Greatest length skull	20.9	0.0001	86.9 ± 2.29	89.8 ± 0.85	88.8 ± 1.32	95.8 ± 0.93
Condylobasal length	16.4	0.0009	81.9 ± 2.21	84.1 ± 0.82	84.5 ± 1.27	89.4 ± 0.90
Zygomatic width	13.2	0.004	66.1 ± 2.49	68.1 ± 0.92	71.3 ± 1.44	73.1 ± 1.02
Mastoid width	22.8	<0.0001	41.3 ± 0.83	43.5 ± 0.31	42.9 ± 0.48	45.4 ± 0.34
Neurocranium width	28.3	<0.0001	43.1 ± 0.57	44.2 ± 0.21	43.1 ± 0.33	45.7 ± 0.23
Occipital height	6.4	0.09	24.2 ± 0.54	23.9 ± 0.20	24.1 ± 0.31	25.0 ± 0.22
Auditory bulla length	8.7	0.03	24.4 ± 0.48	25.5 ± 0.18	25.3 ± 0.28	26.0 ± 0.20
Auditory bulla width	15.4	0.001	14.9 ± 0.36	16.6 ± 0.13	16.6 ± 0.21	17.3 ± 0.15
Auditory bulla height	33.5	<0.0001	17.0 ± 0.35	17.2 ± 0.13	18.1 ± 0.20	18.9 ± 0.14
Postorbital breadth	16.4	0.0009	33.5 ± 0.54	33.8 ± 0.20	32.7 ± 0.31	34.4 ± 0.22
Interorbital breadth	24.8	<0.0001	17.4 ± 0.54	18.1 ± 0.20	18.9 ± 0.31	20.5 ± 0.22
Breadth at canine alveoli	12.6	0.005	20.8 ± 0.54	22.1 ± 0.20	23.0 ± 0.31	23.4 ± 0.22
P4 length	23.2	<0.0001	10.6 ± 0.25	10.8 ± 0.09	11.0 ± 0.14	11.5 ± 0.10
P4 width	12.3	0.006	4.3 ± 0.22	4.6 ± 0.08	4.6 ± 0.12	5.0 ± 0.09
Maxillary toothrow length	15.2	0.002	27.2 ± 0.58	28.2 ± 0.22	28.2 ± 0.34	29.3 ± 0.24
Nasal length	0.83	0.82	26.9 ± 0.90	26.8 ± 0.33	27.0 ± 0.52	27.1 ± 0.37
Mandible length	9.0	0.03	57.2 ± 1.87	60.1 ± 0.69	61.3 ± 1.08	63.1 ± 0.76
Mandible height	15.1	0.001	26.5 ± 1.36	26.5 ± 0.50	27.6 ± 0.79	29.9 ± 0.56
M1 length	20.1	0.0002	7.5 ± 0.19	8.0 ± 0.07	7.8 ± 0.11	8.4 ± 0.08
Mandibular toothrow length	6.0	0.109	29.1 ± 0.85	30.6 ± 0.31	30.5 ± 0.49	31.4 ± 0.35
Mandible depth at M1	17.1	0.0007	9.6 ± 0.45	9.8 ± 0.17	9.9 ± 0.26	11.0 ± 0.18
	females
Greatest length skull	14.5	0.002	81.7 ± 1.44	84.7 ± 0.91	85.3 ± 1.58	90.1 ± 1.11
Condylobasal length	12.1	0.007	76.9 ± 1.31	79.6 ± 0.83	80.2 ± 1.44	83.8 ± 1.02
Zygomatic width	11.8	0.007	62.1 ± 1.43	63.6 ± 0.90	67.0 ± 1.56	68.3 ± 1.11
Mastoid width	12.8	0.005	40.4 ± 0.62	41.7 ± 0.39	41.5 ± 0.68	43.7 ± 0.48
Neurocranium width	13.9	0.003	42.4 ± 0.47	42.8 ± 0.30	42.7 ± 0.51	44.9 ± 0.36
Occipital height	11.1	0.011	22.2 ± 0.30	22.7 ± 0.19	22.30.33	23.60.24
Auditory bulla length	7.9	0.05	23.4 ± 0.34	24.5 ± 0.22	24.1 ± 0.38	24.8 ± 0.27
Auditory bulla width	15.3	0.001	15.5 ± 0.27	15.8 ± 0.17	16.2 ± 0.29	16.9 ± 0.21
Auditory bulla height	19.7	0.0002	16.6 ± 0.25	16.6 ± 0.16	16.8 ± 0.28	18.3 ± 0.20
Postorbital breadth	9.0	0.03	33.6 ± 0.39	34.1 ± 0.25	32.9 ± 0.43	34.7 ± 0.30
Interorbital breadth	21.1	0.0001	15.9 ± 0.38	17.2 ± 0.24	18.1 ± 0.41	18.9 ± 0.29
Breadth at canine alveoli	18.0	0.0004	19.0 ± 0.44	20.4 ± 0.28	21.9 ± 0.48	21.8 ± 0.34
P4 length	19.8	0.0002	9.9 ± 0.17	10.3 ± 0.10	10.3 ± 0.18	11.0 ± 0.13
P4 width	4.04	0.25	4.2 ± 0.14	4.4 ± 0.09	4.6 ± 0.16	4.6 ± 0.11
Maxillary toothrow length	12.7	0.005	25.4 ± 0.37	26.4 ± 0.23	26.6 ± 0.41	27.5 ± 0.29
Nasal length	2.6	0.45	24.8 ± 0.76	25.5 ± 0.48	26.7 ± 0.84	25.7 ± 0.59
Mandible length	10.6	0.014	54.0 ± 1.04	56.2 ± 0.66	58.6 ± 1.14	58.8 ± 0.81
Mandible height	12.4	0.006	23.1 ± 0.76	24.9 ± 0.48	24.7 ± 0.83	27.3 ± 0.59
M1 length	18.6	0.0003	7.1 ± 0.12	7.4 ± 0.08	7.2 ± 0.13	8.0 ± 0.09
Mandibular toothrow length	5.7	0.12	27.6 ± 0.44	28.7 ± 0.22	28.9 ± 0.36	29.2 ± 0.63
Mandible depth at M1	15.2	0.002	8.6 ± 0.24	9.1 ± 0.11	8.7 ± 0.14	9.8 ± 0.23

). Nasal length, maxillary toothrow length, and auditory bulla length showed no significant geographical pattern in either males or females. On average, the sand cats from the northeastern part of the Asian range (thinobia) had the largest skull size, including length and width, and the smallest skull size was recorded in cats living in the Sahara, in the southwest of the species’ range (margarita) (Figure 1).

4. Discussion

Previous studies have demonstrated a lack of morphological and genetic differentiation amongst putative subspecies of larger mammals in North Africa, the Arabian Peninsula, and Southwest Asia, including lions, Panthera leo [47,48], dorcas gazelles, Gazella dorcas [49]. and dama gazelles, Nanger dama [50]. However, these previous studies have focused on large, mobile species, so greater differentiation might be expected in smaller, less mobile species, such as the sand cat, as has been suggested for the desert hedgehog (Paraechinus aethiopicus) [51]. The home ranges of sand cats vary greatly from 7.78 km² to 1363.08 km² for a 100% minimum convex polygon [52], but there are no data on dispersal distances, and it was unclear whether some of the very large home ranges reflected dispersing individuals or not. It is clear that sand cats may move over long distances, but whether these are sufficient to promote widespread gene flow seen in larger species is unclear. Our results suggest that putative subspecies of sand cat are only weakly differentiated or overlap in some characteristics.

Currently four subspecies of sand cat are recognized based on four apparently separated populations, i.e., F. m. margarita in northern Africa, F. m. harrisoni in the Arabian Peninsula, F. m. scheffeli in Pakistan, and F. m. thinobia in Southwest and Central Asia [1,3,4]. More recently, based on a review of current evidence that took into account range extensions in the Middle East and Southwest Asia, and a phylogeographical study by Howard-McCombe et al. [22], this was tentatively reduced to two subspecies, i.e., F. m. margarita in northern Africa and F. m. thinobia in Southwest and Central Asia, including the Arabian Peninsula [21].

Schauenberg [5] was unable to find consistent pelage differences to distinguish between subspecies but did provisionally propose three subspecies based on differences in the size of skulls, i.e., margarita, scheffeli, and thinobia. Hemmer et al. [4] identified various morphological differences between four putative subspecies based on pelages and skulls. These subspecies designations were based on a much smaller overall sample size than the one used in this study and did not involve direct examination of any specimens of thinobia. In our study of almost 100 skins and skulls of sand cats, these suggested diagnostic characters were more variable in each putative subspecies and, hence, we were unable to separate them clearly. Overall, we observed that margarita has smaller skulls than those of the other putative subspecies and that the pelages of margarita and harrisoni are more likely to have distinct markings and lack a distinct dark crest along the mid-line of the dorsum, whereas thinobia has the largest skulls of all putative subspecies and the pelages of both thinobia and scheffeli are less well marked, and both tend to have a distinct crest of darker fur along the mid-line of the dorsum. The number of tail bands ranges from 2 to 7 in scheffeli with a mean of 3.7, which encompassed the variation in all other putative subspecies (

Table 9. Means and ranges of tail bands (TB) in putative subspecies of sand cat.

Subspecies	n	Mean TB	TB	Max TB
margarita	6	3.3	3	4
harrisoni	12	3.7	3	6
scheffeli	56	3.7	2	7
thinobia	6	4	3	5

), so that tail band number is not diagnostic for any putative subspecies (see [4]).

Hemmer et al. [4] distinguished between putative subspecies based on a number of relative skull measurements, including the greatest length of skull, width of skull (zygomatic width), lengths of upper and lower carnassials, occiput height, and auditory bulla height, width, and length. We found that skull breadth did not vary between putative subspecies except for harrisoni, which has a relatively broader skull (zygomatic breadth) (Figure 4a). However, we should note that all the harrisoni specimens we measured were from zoos and had either been bred in captivity or lived most of their lives in captivity. Although most of the scheffeli specimens were also from captivity, these were almost all from the wild and had lived only a short time in captivity. Big cats’ skull morphology appears to be influenced by captivity [53,54,55,56] and even within one generation we must not overlook possible effects of captivity, which may have influenced our results. Typically, captive-bred cats have wider skull breadths (zygomatic breadth) than wild-living cats of the same species, owing to differences in diet and killing/feeding behaviour [53,57].

However, there is no gradual (clinal) geographical pattern of skull size variability running from northeast (thinobia) to southwest (margarita) within the species range, as the skull sizes of scheffeli and harrisoni are very similar in most cases, and the range of scheffeli is significantly further east than that of harrisoni. Our data also do not support the hypothesis of a latitudinal gradient in skull size variability, as the ranges of margarita, scheffeli, and harrisoni are approximately within the same latitudinal range.

Looking more closely at the relative skull measurements previously used to separate subspecies [4], there appear to be no clear patterns in the relative lengths of the upper or lower carnassials, auditory bullae and occiput height between putative subspecies (

Table 10. Means and ranges of lengths of (a) upper (P⁴) and lower (M₁) carnassials and (b) occiput height in putative subspecies of sand cat.

(a)
Subspecies	n	Mean P4 Length	Minimum P4 Length	Maximum P4 Length	n	Mean M1 Length	Minimum M1 Length	Maximum M1 Length
margarita	9	10.14	9.27	10.85	9	7.27	6.67	8
harrisoni	9	10.47	9.27	11.06	9	7.4	6.85	7.65
scheffeli	27	10.61	9.66	11.33	27	7.78	6.99	8.37
thinobia	28	11.32	10.67	12.15	26	8.25	7.54	9.01
(b)
Subspecies	n	Mean Occipital Height	Minimum Occipital Height	Maximum Occipital Height
margarita	7	23.18	21.9	25.53
harrisoni	10	23.23	21.83	24.88
scheffeli	27	23.42	21.52	25.67
thinobia	28	24.48	22.43	26.63

a, Figure 4b–d and Figure 5). However, it should be noted that individuals in the margarita sample have, on average, the shortest carnassials relative to greatest length of the skull. The height of the auditory bulla may be greater in thinobia than in other putative subspecies (Figure 5), but this measurement is difficult to determine accurately and consistently by calipers, even by the same recorder [58], so we must treat these results with caution.

Therefore, in summary, almost none of the diagnostic characters and measurement ratios proposed by Hemmer et al. [4] was able to distinguish between the four putative subspecies in the larger sample available to us, except for zygomatic width, which may have been due to effects of captivity on skull development, and auditory bulla height, which may have been due to inter-recorder error. Indeed, our findings contradicted many of the diagnostic characters proposed by Hemmer et al. [4], including carnassial size, skull breadth, occiput height, and auditory bulla size. The only characters that do provide some distinction are the small overall skull size of margarita and the large overall skull size of thinobia compared with other putative subspecies, and more frequent and bolder pelage markings in margarita and harrisoni compared with fewer, less distinct markings in thinobia and scheffeli, which tend to have a distinct mid-dorsal crest.

Our analyses of geographical variation in the sand cat focused on analyses of four pelage characters and 21 linear skull measurements. The discriminant analysis of four pelage character scores appeared to separate the four putative subspecies well, but it did not reveal any clear geographical patterns of variation consistent with putative subspecies, showing that there was no difference in the four pelage characters between scheffeli and thinobia, although there appears to be discrimination between these and margarita and harrisoni. There is a tendency for northern African and Arabian specimens to be marked with distinct spots and stripes and to lack a dorsal crest, which distinguishes them from specimens from Pakistan and Turkmenistan. The latter have a dorsal crest and are weakly marked with spots and stripes.

In the SZM, only margarita and thinobia occupy clearly separate parts of the morphospace. In addition, thinobia is separated substantially from harrisoni (Figure 3) and the centroid of the margarita sample is significantly separated from the scheffeli centroid along coordinates E3 (males) and E4 (females) (Table 7). However, when interpreting these results, it should be considered that we only examined three males in the margarita sample. We can cautiously suggest the possibility of three subspecies, based on morphometric data only (F. m. margarita, F. m. thinobia, and F. m. scheffeli including harrisoni), but the arguments regarding sampling biases could affect interpretation of the results.

Breton et al. [59] have recently published data on external measurements and weights of Saharan sand cats and found that overall body size was not markedly smaller than in other subspecies, but it is unclear how body size and weight is correlated with skull size.

Howard-McCombe et al. [22] found no significant genetic differentiation between Arabian and Asian (Pakistani and Turkmenian) putative subspecies, especially when captive harrisoni were included, whereas northern African margarita was much more clearly distinct genetically, although its divergence is quite recent, i.e., possibly within the Holocene. Our sample did not include any specimens from intervening areas between thinobia and harrisoni, such as Israel, Jordan, Syria, Iran, and Iraq, and we had only one skin from Egypt; these have mostly not been assigned to putative subspecies and could fill any apparent morphological gaps separating these putative subspecies. Our sample also did not include specimens from Afghanistan and Uzbekistan. The close genetic relationship between thinobia, scheffeli, and harrisoni suggests that the ancestry of these putative subspecies would prove to be much more mixed if samples from intervening geographical areas were to be included. The very low genetic diversity in northern African margarita, compared with sand cats from Arabia and Asia, suggests a very recent range expansion there, reflecting the species’ biogeographical history.

Phylogenetic studies have shown that F. margarita diverged from the common ancestors of F. silvestris/lybica/bieti about 2.37–2.64 Mya (95% confidence intervals 1.79–3.44 Mya) [2], which is approximately at the Pliocene–Pleistocene boundary [60], when the first major glaciations caused expansion of deserts [61,62]. The Pleistocene was dominated by a series of glacial cycles [63]. In general, glaciations increased desert habitats in Asia, thus providing more habitats for the sand cat, whereas interglacials led to a reduction in deserts [63]. During the period from the end of the Pleistocene to early/mid-Holocene (c. 15–5 kya), northern Africa and the Arabian Peninsula experienced monsoon conditions, such that the predominant habitats were savanna and steppe grasslands with large lakes and rivers [64,65,66,67,68,69,70,71]. Savanna mammals were widespread in what is now the Sahara Desert, as evidenced by cave art and fossils, including African elephants, Loxodonta africana, hippopotamus, Hippopotamus amphibius, rhinoceroses, Diceros bicornis and Ceratotherium simum, giraffe, Giraffa camelopardalis, hartebeest, Alcelaphus buselaphus, bohor reedbuck, Redunca redunca, and warthog, Phacochoeurus africanus [72,73,74,75]. It is unlikely that there was much, if any, desert habitat available in the Saharan region and the Arabian Peninsula until climatic deterioration occurred in the mid-Holocene, coincident with the rise of the Ancient Egyptian civilization [67,70]. The refugium for the current sand cat population during the early to mid-Holocene seems to have been Southwest and Central Asia, from where the species recolonized the Arabian Peninsula and subsequently northern Africa, perhaps via Sinai, resulting in a severe population bottleneck and the current weak differentiation of the northern African population. Alternatively, a small, isolated population may have survived somewhere in northern Africa during the Green Sahara period and recolonized the region within the last c. 5000 years. Both these hypotheses are supported by mtDNA, which shows some evidence of weak isolation by distance [22] that may reflect post-Green Sahara and post-Green Arabia recolonization westwards in the early-to-mid Holocene as deserts re-expanded. The sand cat’s possible colonization of northern Africa via Sinai, resulting in a severe bottleneck, could explain the slightly stronger genetic differentiation of nominate margarita and its smaller size, but we cannot rule out small sample sizes for these apparent differences. However, given that our sample included specimens from the Arabian Peninsula and eastern North Africa, there are no “missing” specimens that could fill this morphological gap, so that the size difference is more likely to be real compared with the apparent difference in size between thinobia, harrisoni, and scheffeli. We cannot exclude the possibility that our sample of Arabian sand cat specimens, which was derived from captives, are larger than wild specimens, owing to better nutrition. No wild-caught Arabian sand cat specimens were available to us to test this hypothesis, but we did not find significant differences in sizes between wild and captive specimens of other felids (e.g., [53,76]). Although the Arabian Peninsula is a centre of endemism, the recent arrival of the sand cat and its likely current contiguous distribution with other populations to the north and west suggest that isolation and subsequent subspeciation of Arabian populations are very unlikely, which is supported by genetic data [22].

The onset of Pleistocene glaciations was also coincident with the diversification of small desert rodents, which are the main prey of sand cats, including spiny mice, Acomys spp., and jirds, Meriones spp. and Gerbillus spp. Radiations of spiny mice in dry open habitats of Afro-Arabia may be similar to the speciation of F. margarita; in particular, the divergence and colonization of the Arabian Peninsula by the Acomys russatus complex and its ancestors also occurred at the Pliocene–Pleistocene boundary [77]. The four main sublineages of Gerbillus spp. also diversified at this time [78]. The Libyan jird, Meriones libycus, also has a similar distribution to that of F. margarita; mtDNA analysis and geometric morphometric data [79] identify three allopatric lineages within M. libycus; a western lineage in North Africa (M. l. libycus), a central lineage in Saudi Arabia, Jordan, and Syria (M. l. syrius), and an eastern lineage in Iran, Afghanistan, and China (M. l. erythrourus), which broadly mirror the three main skull morphotypes of sand cat identified in this study. However, based on divergence time estimates, all divergence events within M. libycus probably occurred during the Pleistocene, after 1.597 Ma, which is probably much earlier than any apparent differentiation within F. margarita.

Taking into account both the morphological evidence in this study and genetic evidence from Howard-McCombe et al. [22], it seems likely that only two subspecies of sand cat can be recognized, i.e., Felis m. margarita from northern Africa and F. m. thinobia from the Arabian Peninsula and in Southwest and Central Asia. F. m. margarita can be distinguished from F. m. thinobia by its smaller overall size, and its pelage is more frequently marked with distinct spots and stripes. However, there is some overlap with specimens from the Arabian Peninsula in these characters. Larger samples from the whole geographical distribution of the species may show that no subspecies can be identified, such that Felis margarita is a monotypic species with some ecotypical variation. Nuclear DNA evidence, including whole-genome sequencing, would also help to build on the currently available research by Howard-McCombe et al. [22]. We urge further opportunistic collecting of natural casualties of sand cats from the wild and zoos to expand sample sizes for future research, and welcome information about existing specimens not included in this study, especially from countries unrepresented in our sample.