Phylogenetic relationships of the zokor genus Eospalax (Mammalia, Rodentia, Spalacidae) inferred from whole-genome analyses, with description of a new species endemic to Hengduan Mountains.

Zokors in the genus Eospalax, which are endemic to northern and western China, are subterranean rodents that inhabit various niches, including grasslands, high-altitude meadows, forests, and farmlands. Six species in Eospalax were described a century ago but their taxonomy and phylogeny remain controversial. In this study, we performed high-depth whole-genome sequencing of 47 zokor samples, comprising all six previously described species. Genomic analyses revealed a reliable and robust phylogeny of Eospalax and supported the validity of the six named species. According to the inferred phylogenetic relationships, Eospalax first divergent into two clades in the early Pliocene (ca. 4.68 million years ago (Ma)), one inhabiting the high-altitude Qinghai-Xizang (Tibet) Plateau (QTP) and adjacent regions, and the another inhabiting the low-altitude Loess Plateau and Qinling-Daba Mountains. The most recent divergences occurred between E. baileyi and E. smithii and between E. rufescens and E. rothschildi in the late Pliocene (ca. 2.09 and 2.19 Ma, respectively). We also collected specimens of zokors in the southern Hengduan Mountains (Muli County, Sichuan Province), far from the known distributions of all other zokors. Morphological and molecular analyses strongly suggested that the specimens represent a new species, formally described here as Eospalax muliensis sp. nov. The new species belongs to the high-altitude clade and diverged from closely related species (ca. 4.22 Ma) shortly after the first divergence in Eospalax. Interestingly, Eospalax muliensis sp. nov. possesses more supposedly plesiomorphic characters, suggesting a possible origin of the genus in the Hengduan Mountains.

INTRODUCTION

Zokors (subfamily Myospalacinae) are a group of strictly subterranean rodents endemic to East Asia, which belong to the family Spalacidae and are closely related to the Rhizomyinae and Spalacinae subfamilies (Bouckaert et al., 2019; Guo et al., 2021; Lin et al., 2014; Norris et al., 2004). Zokors are cylindrical in shape, with stocky and strong forelimbs and hardened nose pads to excavate their complex burrow systems. Based on morphological and molecular evidence, extant zokors are classified into two genera: i.e., Myospalax, which includes species with flat occiputs, and Eospalax, which includes species with convex occiputs (Su et al., 2014; Zheng, 1994; Zhou & Zhou, 2008). Myospalax species are mainly distributed in the grasslands of northeast China, Mongolia, and Russia, while Eospalax species inhabit various niches in northern and western China (Figure 1A) (Fan & Shi, 1982), including the high-altitude Qinghai-Xizang (Tibet) Plateau (QTP). The phenotypic and genetic adaptations of the plateau zokor (E. baileyi) to the hypoxic conditions on the QTP have been well researched (Wei et al., 2006; Xu et al., 2021; Zhang et al., 2021).

Figure 1

Sampling of species in Eospalax

A: Geographic distribution of sampling sites. B: Image of typical Eospalax muliensis sp. nov. (upper panel) and aboveground mound (lower panel). Photo by Hao Zhou.

However, despite studies on the physiology, ecology, and genetics of certain species, the taxonomy and phylogeny of Eospalax remain controversial. Six species of Eospalax were described in the late 19th and early 20th century, but the taxonomic status of Eospalax species has changed with time, especially for E. baileyi, E. cansus, and E. rufescens (Table 1). Eospalax fontanierii, E. rothschildi, and E. smithii were recognized as distinct species by Allen (1940), with this classification adopted by many monographs until recently (Smith & Xie, 2008; Wilson & Reeder, 2005). Based on comprehensive morphological comparisons, Fan & Shi ( 1982) found that E. cansus is stably smaller than E. fontanierii in body size and skull and E. baileyi is distinct from E. fontanierii in nose pad shape and presence of a postero-external lobe on M3. Thus, they elevated E. baileyi and E. cansus to species level. Subsequent studies using mitochondrial DNA markers have generally validated the status of all six Eospalax species (He et al., 2012; Su et al., 2014; Zhou & Zhou, 2008). However, the classification of E. baileyi has remained somewhat controversial due to insufficient molecular markers and incomplete sampling (Jiang et al., 2017; Wei et al., 2021). In addition, the potential existence of new undescribed species within the Eospalax genus may complicate established taxonomy. Here, based on a literature review of species distribution, we found reports of zokors in the southern Hengduan Mountains (Muli County, Sichuan Province) (Liu et al., 2007), far from the known distribution of all other zokors (Figure 1A). Thus, in the current study, we collected 14 Eospalax specimens from Muli County between October 2020 and October 2021. Externally, these zokors differed considerably from all other known species due to their small size and long tail (Figure 1B), indicating that they may represent a distinct taxon, which we name herein as Eospalax muliensis sp. nov.

Table 1

Previous and proposed classifications of Eospalax genus

Allen (1940)	Fan & Shi ( 1982)	Song (1986)	Li & Chen ( 1989)	Jiang et al (2017)	Wei et al (2021)	Present study
Recognized species (subspecies) are presented, but not synonyms. *: These two subspecies were not verified in present study but followed Li & Chen ( 1989).
fontanierii	fontanierii	fontanierii	fontanierii	fontanierii	fontanierii	fontanierii
(fontanierii)	cansus	cansus	(fontanierii)	rothschildi	cansus	cansus
(cansus)	(cansus)	rothschildi	(baileyi)	rufescens	baileyi	baileyi
(baileyi)	(rufescens)	rufescens	(cansus)	smithii	rothschildi	smithii
rothschildi	rothschildi	(rufescens)	(rufescens)	cansus	rufescens	rothschildi
smithii	smithii	(baileyi)	rothschildi		smithii	(rothschildi*)
	baileyi		(rothschildi)			(hubeiensis*)
			(hubeiensis)			rufescens
			smithii			muliensis

In addition to taxonomic controversies, the phylogeny of Eospalax remains incomplete and inconsistent. Conflicting phylogenies among Eospalax species have been reported in previous studies based on both mitochondrial DNA markers (Cai et al., 2020; Su et al., 2014; Zhou & Zhou, 2008; Zou et al., 2020) and morphological traits (Li & Chen, 1986; Li & Wang, 1996). Notably, phylogenetic relationships inferred from mitochondrial markers are unstable when different markers or tree-construction methods are applied (Su et al., 2014; Zhou & Zhou, 2008). Therefore, the ambiguous and conflicting phylogenetic relationships reported in Eospalax hamper studies of this remarkable lineage and underscore the need for a reliable and robust phylogeny.

In the current study, we aimed to clarify the taxonomy and phylogeny of Eospalax. All six known species of Eospalax were sampled and whole genomes were sequenced to high depth. Whole-genome analysis revealed reliable and robust phylogenetic relationships and validated the status of all six previously described species. Furthermore, by integrating morphological and genomic evidence, the zokor specimens collected from Muli County were found to represent a new species. Thus, the phylogeny and diversification of the new species were also studied.

MATERIALS AND METHODS

Sample collection

Fieldwork was carried out from April 2017 to October 2021. Skeletal muscle was collected for DNA extraction. All collections followed the animal use protocols approved by the Animal Care and Ethics Committee of the Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences (CAS) (Approval No. SMKX-SQ-2021-067).

Sampling sites of each species were located in their generally acknowledged distribution areas (Figure 1A; Supplementary Table S1). The seven sampling sites of E. baileyi were reported in our previous work (Zhang et al., 2021). Because the distribution of E. cansus is wide and adjacent to several other species in Eospalax, three sampling sites of E. cansus were chosen across Shaanxi, Ningxia, and Gansu provinces to ensure representativeness. For other species, one or two sampling sites were selected. Two zokors collected near Hulun Lake (Inner Mongolia Autonomous Region, China) were recognized as Myospalax psilurus according to their skull and hairless tails. Thus, these two zokors were used as an outgroup.

Whole-genome sequencing and de novo assembly-based variant calling

In this study, 47 zokor samples were collected and whole genomes were sequenced to a high depth, with an average of 137.2 Gb of data for each sample (Figure 1A; Supplementary Table S1). Combined with previously reported whole-genome sequences of 38 E. baileyi samples (Zhang et al., 2021), our dataset covered all six previously described Eospalax species, as well as Eospalax muliensis sp. nov. and two samples of closely related M. psilurus.

DNA extraction, resequencing library construction, and quality control of raw sequencing reads were performed following Zhang et al (2021). Briefly, the resequencing library was constructed following standard protocols and sequenced using Illumina HiSeq or BGISEQ paired-end mode. Raw sequencing reads were trimmed using Trimmomatic v0.36 (Bolger et al., 2014) to remove low-quality reads. Although a E. baileyi genome assembly was recently published (Zhang et al., 2021), all other zokor genomes are currently unavailable. Thus, we sequenced all samples to high depth and applied de novo assembly-based variant calling following the FermiKit pipeline (Li, 2015b). For comparability, variant calling of all E. baileyi samples was also applied using this pipeline. First, clean reads of each sample were error corrected using BFC v.r181 (Li, 2015a). Second, Fermi2 v.r178 (Li, 2012) was used to de novo assemble the genome of each sample. Third, the assembly of each sample was mapped to the plateau zokor genome using minimap2 v.2.17 (Li, 2018) with parameters “-ax asm5/asm10/asm20 -k 23 -w 11 -f 1000 -n 2 -s 100 --end-bonus 5”. According to the genetic distance between samples and plateau zokors, different parameters of -ax were applied to balance the sensitivity and specificity of mapping. Specifically, asm5, asm10, and asm20 were applied for samples of E. baileyi, E. smithii, and all other species, respectively. After mapping, all sample variants were called jointly using the “pileup” functions in HTSbox v.r310 with parameters “-d -V 0.15 -S 300 -q 20 -Q 3 -s 5”. Raw genotypes were first filtered according to allele balance. For homozygous genotypes, at least 90% of reads should support these genotypes. For heterozygous genotypes, the two types of alleles should have 25%–75% of read-support simultaneously. If the above criteria were not satisfied, the genotypes were set as missing. After allelic balance filtering, we only retained those biallelic single nucleotide variants (SNVs) with a missing rate of no more than 20% using VCFtools v.0.1.17 (Danecek et al., 2011).

Phylogenomic analyses

Phylogenetic relationships among all zokor samples were inferred using whole-genome SNVs and maximum-likelihood (ML) algorithms. Considering the computational resources required, we pruned high-quality SNVs by randomly choosing one SNV in each 100 bp window across the genome. ML analyses were conducted in RAxML v.8.2.12 (Stamatakis, 2014) with the nucleotide substitution model of GTRGAMMA. To test the reliability of inferred phylogenetic relationships, 100 bootstraps were applied.

In addition to the construction of phylogenetic trees based on whole-genome pruned SNVs, we also applied a summary of sliding-window ML trees. First, we split the whole genomes into over 50 000 non-overlapping windows (each 50 kb long). For each window, ML trees were constructed using RAxML based on all high-quality SNVs. Weighting and summary of all sliding-window ML trees across whole genomes were explored using topology weighting by iterative sampling of subtrees (Twisst) (Martin & van Belleghem, 2017).

Divergence time among all species in the genus Eospalax was determined using SNAPP v2.6.2 (Bryant et al., 2012) in BEAST2 (Bouckaert et al., 2019). Secondary calibration was applied following He et al (2020), which estimated that the divergence between Myospalax and Eospalax occurred 8.8 million years ago (Ma), with a 95% highest posterior density (HPD) interval from 7.43 to 10.23 Ma. The estimate from He et al (2020) is congruent with fossil records of the genus Myotalpavus, which occurred in the Late Miocene and may be the ancestor of both Eospalax and Myospalax (Zheng, 1994). Thus, we constrained the divergence time between Myospalax and Eospalax with a lognormal distribution centered at 8.8 Ma and a standard deviation of 0.085 in SNAPP analysis. For each species, the individual with the lowest proportion of missing data was used to represent that species. The SNVs were first pruned as described previously, with those identified as monomorphic or containing missing data then filtered. Finally, 100 000 randomly selected SNVs were used, and 10 million Markov Chain Monte Carlo (MCMC) iterations were set for SNAPP analysis. Convergence was assessed using Tracer v1.6 (Rambaut et al., 2018) with 10% burn-in. DensiTree v.2.6.2 (https://www.cs.auckland.ac.nz/~remco/DensiTree/) was used to visualize all posterior trees sampled by SNAPP. TreeAnnotator v2.6.2 (https://beast.community/treeannotator) was used to quantify the posterior probabilities of clades as node support in a maximum clade credibility tree with 10% burn-in. Divergence time and phylogenetic relationships were visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). SNAPP analysis was implemented three times independently and obtained nearly identical results.

Absolute genetic distance (DXY) among species was estimated in 50 kb windows with a 50 kb step size along the genome using popgenWindows.py in genomics_general (https://github.com/simonhmartin/genomics_general). We excluded windows located in scaffolds <1 Mb due to possible poor assembly of short scaffolds.

Mitochondrial genome (mitogenome) analysis

For all whole-genome sequenced samples, we assembled their mitogenomes. We randomly selected 5 Gb of resequencing reads for each sample. Assembly was performed using MitoZ v.2.3 (Meng et al., 2019) or GetOrganelle v.1.7.2b (Jin et al., 2020) with default parameters. The mitogenome of a published M. psilurus sample (JX014234) (Li et al., 2016) served as a reference in GetOrganelle analysis. Annotation of the assembled mitogenomes was performed with MitoZ. The quality of each mitogenome was carefully checked based on length (>16 kb), gene count (37 genes), and circularity.

In addition to the mitogenomes assembled in this study, we downloaded all other mitogenomes in Eospalax and Myospalax reported in previous studies. Thirteen protein-coding genes and two ribosomal RNA (rRNA) genes were included in the following analyses. Gene alignments were performed by CLUSTAL v.2.1 (Larkin et al., 2007) and then concatenated. The best partition schemes for the dataset were identified using PartitionFinder v.2.1.1 (Lanfear et al., 2017) with recommended corrected Akaike Information Criterion (AICc). This analysis partitioned 15 mitochondrial genes into eight subsets. Mitochondrial phylogenomic analyses were inferred using RAxML with 1 000 bootstrap replicates.

Morphological data and analyses

A total of 62 specimens of Eospalax were collected. All specimens and tissue samples were deposited in KIZ, CAS. Specimens were identified following Fan & Shi ( 1982) and Luo et al (2000). Forty-eight specimens were assigned to the six recognized species, including E. baileyi (n=34), E. cansus (n=5), E. fontanierii (n=1), E. rothschildi (n=2), E. rufescens (n=2), and E. smithii (n=4). Fourteen specimens collected from Muli County, Sichuan Province, were assigned to the putative new species (Eospalax muliensis sp. nov.).

The body weight (BW) of each specimen was weighed up to 0.01 g using an electronic scale. External measurements, including head and body length (HB), tail length (TL), and hindfoot length (HF), of some individuals were taken in the field using a ruler to the nearest 0.1 mm. Twenty-one cranial and dental measurements were taken using a digital caliper graduated to 0.01 mm following Pan et al (2007). These measurements included greatest length of skull (GLS), condylobasal length (CBL), basal length (BL), palatal length (PL), zygomatic width (ZMW), interorbital breadth (IOB), foramen infraorbital breadth (FIB), rostrum width (RSW), mastoid width (MTW), maximum width across upper second molars (M2–M2), length of upper tooth row (LUTR), length of upper molars (LUM), nasal length (NSL), braincase height (BCH), greatest breadth of foramen magnum (GBFM), length of auditory bulla (LAB), distance between auditory bulla (DAB), mandibular length (MDL), length of below toothrow (LBTR), length of lower molars (LLM), and length of lower incisor (LLI). We performed principal component analysis (PCA) to analyze morphometric variation of these specimens in SPSS v19.0 (SPSS Inc., USA). PCA was only based on the 21 craniodental measurements of adult specimens. All craniodental variables were log10-transformed before PCA. We considered any component with an eigenvalue exceeding 1.0 to be meaningfully interpretable. We further examined specimen morphology as per Fan & Shi ( 1982), Song (1986), and Luo et al (2000), and followed their terminologies for morphological descriptions.

RESULTS

Phylogenetic relationship and molecular dating

A total of 85 whole-genome sequenced samples were included in this study (Figure 1A; Supplementary Table S1), covering all six previously described Eospalax species as well as Eospalax muliensis sp. nov. and two samples of closely related Myospalax psilurus.

After stringent quality control of genome alignment and variant calling, ~227 million high-quality SNVs were identified. Based on the genome-wide SNVs, ML analysis yielded a highly resolved phylogeny. All nodes that separated species within the genus Eospalax exhibited 100% bootstrap support (Figure 2A). Accordingly, species in Eospalax could be divided into two major clades (i.e., high- and low-altitude clades). The high-altitude clade, comprised of E. baileyi, E. smithii, and Eospalax muliensis sp. nov., generally inhabited the high-altitude QTP and adjacent regions (2 700–4 300 m a.s.l.), while the low-altitude clade, comprised of the other four Eospalax species, generally inhabited the relatively low-altitude Loess Plateau and Qinling-Daba Mountains. Within the high-altitude clade, Eospalax muliensis sp. nov. was found in a basal position and inhabited extremely high altitudes (3 700 m a.s.l.). Within the second clade, E. fontanierii split first from the other species. Notably, three E. cansus populations located remotely from each other formed a tight monophyletic cluster.

Figure 2

Phylogenetic relationships among species of Eospalax

A: ML phylogenetic tree from whole-genome SNVs; 100 bootstraps were applied; scale bar represents level of similarity. B: Pairwise genetic distance among species in genera Myospalax and Eospalax. Genetic distances were calculated in over 50 000 windows (50 kb in size) across whole genome. Data are mean±SD.

The genetic distances between the M. psilurus outgroup and species in Eospalax ranged from 3.5% to 3.7% (Figure 2B). Within the genus Eospalax, the genetic distances among species ranged from 1.0% to 2.2% (Figure 2B). For each species pair separated by the same internal node, their genetic distances were comparable (Figure 2B), demonstrating the reliability of the inferred phylogenetic relationships and that the applied methods enabled phylogenomic analysis.

Phylogenetic relationships inferred from the mitogenomes were also constructed (Supplementary Figure S1). Consistent with nuclear genomic analyses, Eospalax muliensis sp. nov. was also identified as a distinct genetic lineage based on mitogenomic analyses (Supplementary Figure S1). However, significant incongruences among phylogenies inferred from nuclear genomes and mitogenomes were found, especially in the positions of Eospalax muliensis sp. nov. and E. fontanierii (Figure 2A; Supplementary Figure S1) (see Discussion).

To confirm the inferred phylogenetic relationships and investigate the alternative topology, we split the genomes into >50 000 windows (50 kb in size) and constructed window trees based on the corresponding SNVs in each window. Most window trees (75.8%) were consistent with the whole-genome-based phylogenetic tree ( Figure 3A, B). In addition, four alternative topologies were supported by at least 1% of windows (Figure 3A). The top two alternative topologies differed from the majority tree in the position of Eospalax muliensis sp. nov. (Figure 3C,D). The other two less-supported topologies revealed an uncertain relative position of E. fontanierii and E. cansus within the second clade (Supplementary Figure S2).

Figure 3

Alternative topologies and divergence times of genus Eospalax

Average weighting of ranked topologies in over 50 000 windows (50 kb in size) across whole genome (A). Three best topologies and their weighting are shown (B–D). Divergence times estimated using SNAPP based on whole-genome SNVs (E). Node numbers refer to divergence time, blue shadows represent 95% confidence intervals of divergent time. Red circle indicates node for split ofMyospalax and Eospalax as a calibration point.

In divergence time analysis, all posterior trees sampled by SNAPP were identical to the whole-genome-based phylogenetic tree (Supplementary Figure S3). The two major high- and low-altitude clades diverged from each other ~4.68 Ma (95% CI=3.94–5.48 Ma) (Figure 3E). Shortly after, Eospalax muliensis sp. nov. diverged from E. baileyi and E. smithii (~4.22 Ma; 95% CI=3.57–4.97 Ma). Within the low-altitude clade, E. fontanierii first diverged from other species ~3.87 Ma (95% CI=3.26–4.54 Ma). Thus, the speciation events within the genus Eospalax generally occurred in the Pliocene (Figure 3E).

Morphometric comparison and morphological diagnosis

The external and craniodental measurements of each species are given in Table 2. A total of 57 intact skulls of adult specimens were used for PCA. The first three principal components (PC1–3) had eigenvalues greater than 1.0, and together accounted for >79.78% of total variance (Supplementary Table S2). PC1 accounted for 44.08% of the variation (eigenvalue=9.26) and was strongly loaded with most variables, except IOB, GBFM, DAB, and M 2–M2, indicating it may be a marker of size. PC2 accounted for 21.57% of the variation (eigenvalue=4.53) and was strongly dominated by RSW, IOB, FIB, and LLI (loading>0.62), while PC3 accounted for 14.13% of the variation (eigenvalue=2.97) and was dominated by GBFM, DAB, and M2–M2.

Table 2

Body weight and external and cranial measurements (mm) (mean±SD and ranges) of Eospalax specimens (n) examined in this study

Variable	Eospalax muliensis		E. baileyi	E. cansus	E. fontanierii	E. rothschildi	E. rufescens	E. smithii
Character abbreviations (unless specifically, all in mm): BW: Body weight (grams); HB: Head and body length; TL: Tail length; HF: Hindfoot length; GLS: Greatest length of skull; CBL: Condylobasal length; BL: Basal length; PL: Palatal length; ZMW: Zygomatic width; IOB: Interorbital breadth; FIB: Foramen infraorbital breadth; RSW: Rostrum width; MTW: Mastoid width; M2–M2: Maximum width across the upper second molars; LUTR: Length of upper tooth row; LUM: Length of upper molars; NSL: Nasal length; BCH: Braincase height; GBFM: Greatest breadth of the foramen magnum; LAB: Length of auditory bulla; DAB: Distance between auditory bulla; MDL: Mandibular length; LBTR: Length of below toothrow; LLM: Length of lower molars; LLI: Length of lower incisor. *: The smaller value is from a subadult individual. N/A: Not available.
	n=14		n=34	n=5	n=1	n=2	n=2*	n=4
BW	155.64±43.78 118.24–280.22; 12		248.44±73.20 133.63–384.33; 34	N/A	N/A	N/A	N/A	N/A
HB	168.92±14.64 145.00–204.00; 12		191.47±16.03 160.00–220.00; 34	N/A	N/A	N/A	N/A	N/A
TL	52.86±8.82 45.00–70.00; 7		N/A	N/A	N/A	N/A	N/A	N/A
HF	28.00±2.08 25.00–32.00; 7		N/A	N/A	N/A	N/A	N/A	N/A
GLS	40.52±1.83 38.01–45.52; 14		44.43±2.46 40.57–50.18; 32	44.38±3.07 41.45–49.05; 5	50.94	41.83±3.17 39.58–44.07; 2	40.60±4.10 37.70–43.50; 2	48.22±1.88 45.56–49.63; 4
CBL	38.15±1.90 35.55–43.34; 13		42.46±2.53 37.84–48.67; 32	42.17±3.32 39.33–47.09; 5	49.36	39.62±2.26 38.02–41.21; 2	37.82±4.38 34.72–40.91; 2	45.81±1.79 43.43–47.30; 4
BL	35.87±1.83 33.30–40.76; 13		40.14±2.59 35.55–46.44; 32	39.88±3.15 37.26–44.62; 5	46.29	37.58±2.31 35.95–39.21; 2	35.85±4.71 32.52–39.18; 2	43.20±1.86 40.78–44.86; 4
PL	27.02±1.21 25.60–30.39; 14		30.63±1.83 26.62–34.73; 32	30.74±2.32 28.72–34.18; 5	35.93	31.33±2.29 29.71–32.95; 2	26.84±2.49 25.08–28.60; 2	33.07±1.87 30.46–34.59; 4
ZMW	26.25±1.75 23.57–30.60; 14		30.83±2.91 25.75–36.58; 32	29.51±3.21 26.90–33.78; 5	37.01	29.06±3.73 26.42–31.70; 2	27.52±4.51 24.33–30.71; 2	34.76±2.93 30.47–36.88; 4
IOB	7.11±0.32 6.61–7.64; 14		7.80±0.37 6.95–8.37; 32	7.34±0.36 7.01–7.91; 5	8.53	7.81±0.33 7.57–8.04; 2	7.03±0.48 6.69–7.37; 2	8.04±0.27 7.69–8.33; 4
FIB	6.24±0.23 5.80–6.72; 14		8.67±0.76 7.32–10.28; 32	7.73±0.57 7.22–8.35; 5	8.62	7.60±1.03 6.87–8.33; 2	7.10±0.59 6.68–7.52; 2	9.07±0.79 8.04–9.75; 4
RSW	8.99±0.50 8.32–10.31; 14		11.54±0.77 10.22–13.31; 32	10.43±0.95 9.51–11.68; 5	13.84	10.10±1.19 9.26–10.94; 2	9.54±1.75 8.30–10.78; 2	12.26±0.97 11.13–13.23; 4
MTW	23.54±1.31 21.55–26.74; 14		27.70±2.25 23.44–31.77; 32	26.28±1.93 24.09–28.38; 5	35.60	26.11±3.39 23.71–28.50; 2	22.42±2.75 20.47–24.36; 2	29.99±2.31 26.59–31.69; 4
M2–M2	8.82±0.45 8.21–9.74; 14		9.02±0.44 8.37–10.33; 32	8.36±0.34 7.96–8.86; 5	10.21	8.68±0.52 8.31–9.04; 2	8.16±0.73 7.64–8.67; 2	9.24±0.62 8.42–9.88; 4
LUTR	23.94±1.28 22.33–27.45; 14		27.24±1.56 24.50–30.55; 32	27.19±2.12 24.97–29.96; 5	32.47	25.38±1.57 24.27–26.49; 2	24.29±3.00 22.17–26.41; 2	29.29±1.57 27.61–30.86; 4
LUM	8.95±0.46 8.33–10.26; 14		10.28±2.88 9.14–25.93; 32	10.34±0.65 9.70–11.19; 5	11.93	9.80±0.02 9.78–9.81; 2	9.34±1.02 8.62–10.06; 2	10.20±0.66 9.51–10.99; 4
NSL	14.55±1.17 12.41–17.37; 14		17.09±1.24 14.67–19.29; 32	16.40±1.08 15.17–17.98; 5	18.13	15.68±2.14 14.16–17.19; 2	15.60±2.14 14.08–17.11; 2	18.96±1.15 17.38–20.02; 4
BCH	16.20±0.66 15.35–17.80; 14		17.53±1.10 15.86–19.63; 32	16.87±1.51 15.69–18.99; 5	20.94	17.68±2.73 15.75–19.61; 2	14.47±0.78 13.91–15.02; 2	17.63±0.67 16.80–18.18; 4
GBFM	6.51±0.20 6.25–6.99; 13		6.63±0.28 6.04–7.05; 32	6.52±0.33 6.08–6.84; 5	6.96	6.25±0.18 6.12–6.37; 2	6.21±0.18 6.08–6.33; 2	6.58±0.52 5.91–7.19; 4
LAB	8.32±0.39 7.81–9.13; 14		9.27±0.57 8.44–11.08; 32	10.06±0.49 9.26–10.52; 5	11.09	9.19±0.81 8.61–9.76; 2	9.48±0.64 9.03–9.93; 2	10.48±0.25 10.14–10.75; 4
DAB	4.24±0.28 3.75–4.76; 14		4.93±0.89 4.01–9.21; 31	3.73±0.48 3.23–4.29; 5	5.77	3.75±0.29 3.54–3.95; 2	3.15±0.62 2.71–3.58; 2	4.91±0.76 3.88–5.72; 4
MDL	25.17±1.25 23.63–28.48; 14		26.95±1.50 23.71–29.86; 32	28.12±2.49 25.61–31.44; 5	33.40	26.28±1.75 25.04–27.51; 2	25.93±2.91 23.87–27.99; 2	28.59±1.53 27.01–30.08; 4
LBTR	19.65±1.27 18.05–22.68; 14		22.12±1.43 18.55–25.18; 32	21.67±1.43 20.11–23.33; 5	26.23	20.14±2.05 18.69–21.59; 2	19.93±2.86 17.90–21.95; 2	25.15±1.05 24.17–26.25; 4
LLM	9.43±0.46 8.81–10.34; 14		10.02±0.35 9.29–10.57; 32	10.39±0.61 9.68–11.23; 5	12.39	10.05±0.06 10.00–10.09; 2	10.25±1.16 9.43–11.07; 2	10.73±0.59 9.96–11.39; 4
LLI	8.59±1.51 6.40–10.84; 14		11.61±2.00 7.97–15.95; 32	10.11±1.72 7.80–12.05; 5	13.53	8.82±2.65 6.94–10.69; 2	10.22±2.79 8.24–12.19; 2	15.63±0.86 14.83–16.53; 4

The PC1–2 plot (Figure 4) showed that Eospalax muliensis sp. nov. occupied the lower-left corner and was well separated from the other species. Specimens of Eospalax muliensis sp. nov. were plotted on the negative field of both PC1 and PC2, reflecting its generally smaller size, narrower RSW, IOB, and FIB, and shorter LLI. Specimens of E. baileyi occupied a large area, from –1.34 to 1.45 in PC1 and –0.49 to 1.70 in PC2, indicating great variation in size, with a relatively wide RSW, IOB, and FIB, and long LLI. Specimens of E. cansus, E. fontanierii, E. rufescens, and E. smithii were mostly plotted on the positive region of PC1, indicating relatively large skulls.

Figure 4

Plot of principal components 1 and 2 from analysis of 21 craniodental measurements of Eospalax genus

Both molecular and morphological evidence supported Eospalax muliensis sp. nov. as a new species, which is formally described below.

Taxonomic account

Zhang, Chen & Shi, sp. nov.

Suggested common name: Muli zokor; 木里鼢鼠(Muli Fenshu)

Holotype: KIZ 040324 (field No.: SC 2110385). Adult female collected by Tao Zhang, Meng-Long Lei, Zhong-Zheng Chen, and Hao Zhou on 21 April 2021. Dried skin, cleaned skull, and alcohol-preserved carcass are deposited in KIZ. Measurements (mm): BW=133.00 g, HB=172.0, TL=56.5, HF=28.0, GLS=39.41, CB=37.46, BL=35.40, PL=27.27, ZMW=25.58, IOB=7.27, FIB=6.12, RSW=8.84, MTW=22.95, M2–M2=9.30, LUTR=23.93, LUM=9.11, NSL=14.65, BCH=15.84, GBFM=6.25, LAB=8.33, DAB=4.12, MDL=24.75, LBTR=20.36, LLM=9.63, LLI=10.51.

Type locality: Kangwu Ranch (N28.135°, E101.196°), Muli County, Sichuan Province, China, adjacent to Shangri-La Lake at an altitude of ~3 700 m a.s.l.

Paratypes: Seven specimens were collected from the type locality. Four specimens (KIZ 040325–040328) were deposited in KIZ, CAS; three specimens (SC 2110384, 2110386, and 2110451) were deposited in the Biological Museum of Anhui Normal University. All specimens were prepared as dried skins with cleaned skulls and alcohol-preserved carcasses.

Etymology: The species name muliensis is derived from Muli County, the type locality of the new species in Sichuan Province, China, and the Latin adjectival suffix -ensis means “belonging to”.

Diagnosis: Eospalax muliensis sp. nov. can be distinguished from other described species of Eospalax by a combination of the following characters: size small (mean BW=155.64 g; mean GLS=40.52 mm); tail relatively long, densely hairy; nose-pad trifoliate; rostrum weak, nearly rectangular, nasals small, posterior border of nasals transverse, anterior halves of premaxilla nearly aligned with nasals; braincase well-domed, temporal ridges not conspicuous, parallel in front, lambdoid crests only present on side of skull; occipital shield with developed occipital ridges, extending well posteriorly, almost forming plane with occipital condyle; about 1/3 of incisive foramina included within maxillae, remaining 2/3 of incisive foramina included within premaxillae; palate and pterygopalatine fossa large, M2–M2 almost equal to LUM; M3 with two reentrant angles on outer side.

Description: Small-sized zokor (Table 2), with head and body length <175 mm (range: 145.0–175.0 mm) and skull length <42.20 mm (range: 38.01–42.17 mm). Only one elderly male (KIZ 0403327) significantly larger, with head and body length of 204 mm and skull length of 45.52 mm. Dorsal pelage dark grayish brown, with slightly cinnamon-colored tips, ventral pelage slightly paler; most individuals without white blaze on forehead, small when present. Lips and muzzle white, surrounded by short white hairs. Nose-pad trifoliate. Eyes very small and external ears absent ( Figure 1B). Tail relatively long, densely hairy, proximal 2/3 grayish-brown, and distal 1/3 white. Forepaws strong and powerful; third claw longest and stoutest; second and fourth claws almost equal in length, about 2/3 of third claw; fifth claw stout, about half length of fourth finger; thumb very small. Dorsal surfaces of hands and feet covered with short white hair.

Skull relatively weaker than other zokors (Figure 5). Rostrum relatively weak, thin, and nearly rectangular; nasals relatively small and trapezoidal, extending far beyond incisors, posterior border of nasals transverse. Premaxilla undeveloped, anterior half not expanded, nearly aligned with nasals. Zygomatic arches slightly expanded and thin; zygomatic plates narrow and simple. Interorbital region slightly constricted, almost same width as rostrum. Braincase well-domed, rounded ventrally. Temporal ridges not conspicuous (but much more developed in old male KIZ 0403327), parallel in front. Lambdoid crests sharp, but only on side of skull; occipital shield and occipital ridges developed, extending well posteriorly, almost forming plane with occipital condyle; small medial occipital ridge present between occipital ridges. Incisive foramina small, about 2/3 included within premaxilla, posterior 1/3 within maxillaries. Palate and pterygopalatine fossa broad, relatively wider than other Eospalax species. Body of mandible curved. Coronoid process weak, pointed, and curved to posterior. Condyloid process stout and squared. Angular process small, curved upward almost joining with ramus.

Figure 5

Dorsal, ventral, and lateral views of skull and mandible of Eospalax muliensis sp. nov. (KIZ 040324, holotype) (A), E. baileyi (DF 004) (B), E. cansus (LM 001) (C), E. fontanierii (FS 001) (D), E. rothschildi (ZB 001) (E), E. rufescens (FP 001) (F), and E. smithii (XZ 001) (G)

Teeth relatively light. Incisor teeth stout with tips curving slightly backward, outer side orange. Upper molar tooth rows slightly divergent to posterior; M1 with two reentrant angles on outer and inner side; M2 and M3 with two reentrant angles on outer side, and one inside, no postero-external lobe on M3; lower incisor short, slightly curved, outer side orange; M1 with two reentrant angles on outer side and three on inner side; M2 and M3 with two reentrant angles on outer and inner sides.

Comparison: According to the measurements, Eospalax muliensis sp. nov. is similar in size to E. rothschildi, but much smaller than other species in all dimensions (Table 2; Supplementary Table S3). The new species has some features unique in the genus. The skull of Eospalax muliensis sp. nov. is obviously weaker than that of the other Eospalax species, especially in the rostrum region (Figure 5). Therefore, the rostrum-related measurements (i.e., RSW, FIB, and NSL) of Eospalax muliensis sp. nov. are the smallest amongEospalax species (Table 2). The nasals of the new species are very short and narrow, the smallest in the genus. The premaxilla of Eospalax muliensis sp. nov. is depressed and very thin, with the anterior half parts nearly aligned with the nasals, in contrast with other species, in which the premaxillae are well-developed and broader than the nasals. The palate and pterygopalatine fossa are relatively wider in Eospalax muliensis sp. nov than in other species, thus M2–M2 is almost equal to LUM in the new species but is much smaller than LUM in other species. The braincase of Eospalax muliensis sp. nov. are well-domed and temporal ridges being ambiguous, characteristics only seen in some individuals of E. cansus and different from the flattened braincase with developed temporal ridges seen in other species.

Eospalax muliensis sp. nov. can be distinguished from small-sized E. rothschildi by its strong and powerful forepaws (vs. slender); longer tail TL=45–70 mm (vs. short tail TL=23–38 mm, Supplementary Table S3); and proximal 2/3 of tail grayish-brown and distal 1/3 white (vs. gray above white below). Eospalax muliensis sp. nov. can be distinguished from E. fontanierii and E. cansus based on trifoliate nose-pad and densely hairy tail (vs. oval nose-pad and nearly naked tail). Eospalax muliensis sp. nov. can be clearly distinguished from E. baileyi and E. smithii based on incisive foramina about 1/3 expanding premaxillary boundaries (vs. incisive foramina wholly contained in premaxillary). Eospalax muliensis sp. nov. can be distinguished from E. baileyi by the postero-external lobe on M3; the parallel temporal ridges of E. muliensis sp. nov. is distinct from the temporal ridges meeting in the median line of E. smithii. The zygomatic arches of E. rufescens are much stouter and flare wider than Eospalax muliensissp. nov.

Distribution and habitat: Eospalax muliensis sp. nov. is currently known only from the Kangwu Ranch, Muli County, Sichuan Province, China (elevational range ~3 700 m a.s.l.). The habitat is alpine meadow, surrounded by shrubs (Supplementary Figure S4). This species may be adapted to high-elevation habitats, and thus may occur in other alpine meadows in the southern Hengduan Mountains.

DISCUSSION

Taxonomic implications

Both morphological and molecular analyses strongly suggested that the zokors inhabiting Muli County represent a new species (Eospalax muliensis sp. nov.) of Eospalax. The new species can be diagnosed from all others within Eospalax by the combination of small size, long tail, weak rostrum and nasal, and large palate and pterygopalatine fossa. Genomic analyses also revealed that Eospalax muliensis sp. nov. diverged from other recognized zokor species ~4.22 Ma, with high genetic distances in nuclear and mitochondrial genomes.

Our results also strongly supported the status of all six previously described Eospalax species. The most recent divergences in Eospalax occurred between E. baileyi and E. smithii, and between E. rufescens and E. rothschildi. Evidence revealed substantial differences among the two pairs of sister species and supported their species status. First, based on whole-genome analysis, the two species pairs diverged 2.09 Ma and 2.19 Ma (Figure 3E), with relatively high genetic distances of 1.05% and 1.02% (Figure 2B), respectively. Genetic distances of the mitochondrial cytb gene between sister species of these pairs reached 11.9% and 12.5%, respectively (Supplementary Table S4) (Bradley & Baker, 2001). Second, the morphological characteristics that distinguish these sister species have been long recognized. For example, E. smithii can be distinguished from E. baileyi based on temporal ridges meeting at the median line (Figure 5) (Fan & Shi, 1982); and E. rothschildi can be distinguished from E. rufescens based on temporal ridges conspicuous and raised (Figure 5) (Fan & Shi, 1982). Third, the distribution of habitats differs among these sister species. Eospalax baileyi is dominant on the QTP (altitudinal range of 2 700–4 300 m a.s.l.), while E. smithii inhabits a small area around the Min Mountains (eastern edge of QTP) at an altitude of 2 700 m a.s.l. and possibly around the Liupanshan Mountains (within the Loess Plateau) at much lower altitudes (Fan & Shi, 1982). In the second sister species pair, E. rufescens occurs in the Qinling Mountains, while E. rothschildi occurs in the Daba Mountains and probably the Min Mountains. Although the forest habitats of these two species are similar, their local climate and flora differ, and they are likely geographically isolated by the Hanjiang River. In summary, genomic, morphological, and biogeographical evidence well supports the delimitation of these two pairs of recently diverged sister species. Furthermore, the status of other species within Eospalax that were divergent longer before should also be retained.

Phylogenetic implications

Utilizing whole-genome sequencing data, our study revealed reliable and robust phylogenetic relationships among species in the genus Eospalax. Before discussing the implications of the whole-genome-based phylogenetic relationships, we first explored the significant incongruences among phylogenies inferred from nuclear genomes, mitogenomes, and previous studies. Regarding phylogenies inferred from mitogenomes and mitochondrial genes, our results and previous studies show general similarity in the close relationship between E. baileyi and E. smithii and the basal position of E. fontanierii (Supplementary Figure S1) (Cai et al., 2020; He et al., 2012, 2020; Su et al., 2014). However, extensive inconsistencies appear in mitochondrial analyses, which may be due to the limited informative sites of mitochondrial genes and complex speciation and post-speciation processes in the genus Eospalax. For the latter, some clues may be present in our mitogenomic analysis that integrates the mitogenomes in our and previously published data. Our mitogenomic analyses revealed a robust phylogeny but embedded one E. cansus (KC514112) individual tightly within E. baileyi (Supplementary Figure S1). This individual was reportedly sampled from Lintan County, Gansu Province (around Min Mountains) (Su et al., 2013), which contains the type localities of E. cansus, E. rothschildi, and E. smithii. Furthermore, E. baileyi habitats are also adjacent to this region. Besides the possible misidentification of the specimen or habitat expansion of E. baileyi in Lintan, another plausible explanation may be post-speciation gene flow between E. cansus and E. baileyi around the Min Mountains. If this hypothesis is true, the evolutionary history of zokors around the Min Mountains may be more complex, which requires further research.

Several incongruences in the phylogenetic topologies inferred from nuclear and mitochondrial genomes were also found. Notably, nuclear genomic analysis placed E. fontanierii within the low-altitude clade (Figure 2A), whereas mitogenomic analysis and most previous studies placed E. fontanierii in a basal position, while Eospalax muliensis sp. nov. were put aside (Supplementary Figure S1). In addition, nuclear genomic analysis placed Eospalax muliensis sp. nov. within the high-altitude clade (Figure 2A), while mitogenomic analysis placed it in the basal position of the genus (Supplementary Figure S1). Nuclear and mitochondrial discordance appears to be prevalent in phylogenetic studies and the possible roles of rapid speciation and incomplete lineage sorting have been discussed extensively (He et al., 2010; Mckay & Zink, 2010). Rapid speciation and incomplete lineage sorting may have contributed to the inconsistent relationships found among E. cansus, E. rufescens, and E. rothschildi, given the rapid cladogenesis of these three species (Figures 2A, 3E). Other factors may also have contributed to the inconsistent positions of E. fontanierii and Eospalax muliensis sp. nov. For example, the clear pattern of mitochondrial phylogenomic analyses is the relevance between geographical distributions and phylogenetic positions of the species in Eospalax. Eospalax fontanierii and Eospalax muliensis sp. nov. are distributed in the remote northeast and southwest regions (Figure 1A), and their phylogenetic positions are basal. Furthermore, other species are geographically and phylogenetically adjacent. The possible explanation for this pattern may be ancient cross-species hybridization between geographically adjacent species or some other geographically relevant factors (Toews & Brelsford, 2012). However, these hypotheses lack solid evidence and require further comprehensive investigations. Regardless of the actual source of incongruence, the phylogeny inferred from nuclear genomes was more robust and thus serves as our current ‘‘best” estimate (Figure 2A).

From the inferred nuclear whole-genome-based phylogeny, Eospalax was first divided into high- and low-altitude clades. Because uplift of the QTP may predate the separation of these two clades (Spicer et al., 2021) and no other species in the subfamily Myospalacinae or family Spalacidae inhabit the QTP, adaptation to high altitude may have occurred in the common ancestor of Eospalax muliensis sp. nov., E. baileyi, and E. smithii. Nowadays, E. smithii inhabits areas at relatively low altitudes, implying possible re-adaptation to low elevations. Notably, given the possible existence of other unrecognized species in the Hengduan Mountains, adaptation to high altitude could have occurred in the common ancestor of all Eospalax species. Within the low-altitude clade, successful habitation of forest niches appeared in E. rufescens and E. rothschildi. Because these two species are closely related and separated more recently (Figures 2A, 3E), forest adaptation likely occurred in their common ancestor. Thus, although the dispersal capacity of zokors is extremely limited (Wei et al., 1997), species in Eospalax successfully and rapidly adapted to diverse ecological niches in the Pliocene, and further investigations on their adaptive mechanisms are warranted.

Implications of systematics and macroevolution from Muli zokors

The discovery of Eospalax muliensis sp. nov. in the southern Hengduan Mountains is unexpected and of great interest because the habitats of Eospalax muliensis sp. nov. are southernmost of extant zokors and it is the second extant zokor species that inhabit extremely high altitudes (~3 700 m a.s.l.). The Hengduan Mountains are considered a global biodiversity hot spot (Myers et al., 2000). The extremely diverse ecological habitats in three dimensions of Hengduan Mountains (Sky islands) protected species from climate fluctuations and are known inhabited by some relict species (Chen et al., 2017; Chen et al., 2021; He & Jiang, 2014). The discovery of Eospalax muliensis sp. nov. in such a biodiverse area suggests the possible existence of other unknown zokor species in the Hengduan Mountains. On the other hand, the discovery could also provide clues for deciphering the origin and evolutionary history of the genus Eospalax.

Based on the paleontological rules and evolutionary trajectories of zokor fossils, Li & Chen ( 1986) proposed nine plesiomorphic characters of ancient Eospalax species. Strikingly, of the nine characters, six are found in Eospalax muliensis sp. nov., in contrast to the five or fewer in other species in Eospalax. For example, body size trajectory in zokor fossils is from small- to large-sized (Zheng, 1994), whereas Eospalax muliensis sp. nov. is among the smaller ones in Eospalax. Eospalax muliensis sp. nov. also possesses other plesiomorphic characters such as: long tail; weak rostrum and nasal; slightly expanded and thin zygomatic arch; ambiguous temporal ridge; and separated temporal ridge. Absent of Eospalax muliensis sp. nov., previous studies proposed the Loess Plateau or Qinling-Daba Mountains as the possible origin center of Eospalax (Li & Chen, 1986; Li & Wang, 1996). Considering the plesiomorphic characters and relatively basal phylogenetic position of Eospalax muliensis sp. nov., it is possible that Eospalax may have originated from the Hengduan Mountains. This hypothesis seems to contradict fossil evidence, as the earliest known Eospalax fossils are recorded from the Nihewan Basin, Hebei Province (Zheng, 1994) and are dated to ~2 Ma (MN18) (Li et al., 1984). However, the estimated diversification of Eospalax species occurred ~4.68 Ma, and thus the contradiction may be caused by an incomplete Eospalax fossil record, especially in the Hengduan Mountains. Therefore, the origin and evolutionary history of Eospalax remain obscure and require further interdisciplinary evidence.

SCIENTIFIC FIELD SURVEY PERMISSION INFORMATION

All collections followed the animal use protocols approved by the Animal Care and Ethics Committee of the Kunming Institute of Zoology, Chinese Academy of Sciences (Approval No. SMKX-SQ-2021-067). Permission for field sampling in each sampling sites was granted by the local Forestry Bureau.

DATA AVAILABILITY

All sequences reported in this study were deposited in the Genome Sequence Archive database (http://gsa.big.ac.cn/) under Accession ID (CRA005933) and in Science Data Bank (https://www.scidb.cn/) under DOI: 10.11922/sciencedb.o00023.00002. Mitogenomes were deposited in GenBank under Accession Nos. OM329951–OM330035.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

P.S. conceived and supervised the study. T.Z., M.L.L., H.Z., and Z.Z.C. collected samples. Z.Z.C. and T.Z. performed morphometric comparison and morphological diagnosis. T.Z. analyzed whole-genome resequencing data. T.Z., Z.Z.C., and P.S. wrote the paper. All authors read and approved the final version of the manuscript.