A total of 214 specimens of N. confucianus species complex was collected from 35 sampling sites in China using Sherman living cages from March 2009 to August 2018. The sampling sites covered the distribution range of N. confucianus (Fig. 1), and sample information is shown in Suppl. material 1: Table S1. Chromosome preparations were made by the live bone marrow method (Searle 1986). All specimens including the pelt, carcass, and skull were stored at Shandong University (Weihai). All animal sample collection protocols complied with the current laws of China and all animals were handled in a manner consistent with the guidelines approved by the American Society of Mammalogists (Sikes et al. 2016).

Figure 1. 

Distribution of phylogenetic clades of N. confucianus species complex obtained from Cytb. The numbers correspond to the locality code in Suppl. material 1, Table S1.

DNA was extracted from muscle samples using the Easy Pure Genomic DNA Kit (TransGen Biotech Co., Ltd., Beijing, China). The complete mitochondrial cytochrome b gene (Cytb, 1140 bp) was amplified by PCR using the primers described by Irwin et al. (1991). The primer sequences were as follows: Nivicob1 (5'-TGTCATTATTTCTACACAGCACTTA-3') and Nivicob2 (5'-TTTGGGTGTTGATGGTGGG-3'). PCR was performed in a volume of 50 μL containing 30 ng template DNA, 2×EasyTaq PCR SuperMix 25 μL, and 0.5 μM primers. The thermocycling protocol was as follows: initial denaturation of 4 min at 94 °C; 32 cycles of 94 °C for 30 s, annealing temperature (Tm) for 30 s, 72 °C for 70 s; and final extension for 6 min at 72 °C. Detection was carried out by 1% agarose gel electrophoresis, and PCR products were directly sequenced by Sanger sequencing.

We supplemented our new Cytb data with homologous sequences (>1,140 bp) of Niviventer available in GenBank (Balakirev and Rozhnov 2010; Chen et al. 2011; Zhang et al. 2016). Here, six species were selected as outgroup taxa: Bandicota indica, Berylmys bowersi, Leopoldamys nielli, Maxomys surifer, Rattus norvegicus, and Rattus andamanensis (Michaux et al. 2007; Balakirev et al. 2012, 2013; Koma et al. 2013; Conroy et al. 2013). GenBank accession numbers for the original sequences used in this study are MT333860-MT334073 (Suppl. material 1: Table S1).

All sequences were aligned with Clustal X 2.0 (Larkin et al. 2007), manually edited in BioEdit 7.2.5 (Hall 1999), and corrected to eliminate interference from degenerate bases. We used the Akaike Information Criterion (AIC) in jModeltest 1.0 (Posada 2008) to select the best-fit model of sequence evolution for locus alignment. We calculated the population haplotypes in DnaSP (Librado and Rozas 2009) and performed phylogenetic reconstructions based on Cytb using the maximum likelihood (ML) and neighbor-joining (NJ) approaches with the TN93+G+I model and K2P model in MEGA 7 (Kumar et al. 2016), respectively. Bootstraps were obtained using a rapid bootstrapping algorithm with 1000 replicates. We also constructed a phylogenetic tree using Bayesian inference (BI) in MrBayes 3.2 (Ronquist et al. 2012, Balakirev et al. 2013) based on the TN93 model. This step was repeated twice, the replacement rate of the sequence with invgamma was determined, and the process was conducted four Markov chain Monte Carlo runs with four chains for 10 million generations, sampling every 1000 trees and discarding the first 25% as burn-in. We calculated Kimura-2-parameter (K2P) distances of Cytb in MEGA 7 (Kumar et al. 2016) for pairwise comparisons of genetic differentiation within and between different phylogenetic lineages, and standard error was analyzed using 1000 bootstrap tests.

We used ABGD (Puillandre et al. 2012) to recover candidate species. All aligned haplotype sequences (166 Cytb) were uploaded to the web interface (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) and run with the following settings: P (prior limit to intraspecific diversity) range of 0.001–0.1 and relative gap widths (X) of 0.5, 1.0, 1.5, 2.0, and 2.5. Transition/transversion bias (TS/TV) was estimated using MEGA 7. We selected the Kimura 80 model to analyze our data and set the number of both steps and bins to 25.

We captured and analyzed the mitotic phase with improved chromosome dispersion using Cytovision System (Applied Imaging, Newcastle upon Tyne, UK). The diploid number (2n) and chromosome fundamental arm number (FN) were determined in each karyotype. Chromosomes were classified according to Levan et al. (1964) and Motokawa et al. (2001) to analyze the differences in karyotypes among different clades.

To understand the morphological diversity of N. confucianus, N. c. sacer, N. bukit and N. lotipes, which are closely related species in N. confucianus species complex, we explored differences in external and skull morphology among the four species, in which the measurements data of N. bukit was reference to Ge et al. (2018b), the other three were collected and measured by our laboratory. We determined four external indices (head and body length (HBL); tail length (TL); ear length (EL); hind foot length (HFL)) and eight skull indices (greatest length of skull (LS); zygomatic width (ZW); interorbital breadth (IOB); breadth of rostrum (BR); length of incisive foramen (LIF); length of upper tooth row (LUTR); length of auditory bulla (LAB); length of upper diastema (LD)) according to Huang et al. (1995), Yang et al. (2005), Xia et al. (2006), and Ge et al. (2018b). Skull indices were measured with a digital vernier caliper (0.01 mm). Adult specimens were identified by the full eruption of the molars and the frequency histogram of head and body length (Yang 1990, Li et al. 1989, 1990). Male and female specimens were mixed in analyses since previous studies identified no significant differences in the external and cranial measurements between sexes (Stefen and Rudolf 2007; Yang et al. 2011).

To characterize the differences among the four species, standard statistics including the mean and standard error were applied. Pairwise differences between major species groups were tested by an analysis of variance (ANOVA) using least significant difference (LSD) tests, as LSD is more commonly used and sensitive to obtain statistical differences, and multivariate analysis (principal component analysis, discriminant analysis, cluster analysis) was also performed. These analyses were performed using SPSS Statistics 24.0 (SPSS, Chicago, IL, USA).

We photographed the skulls of specimens of N. confucianus, N. c. sacer and N. lotipes for quantitative analysis (as N. bukit was not sampled). The dorsal, ventral, and lateral sides of the skull and lateral view of the mandible were analyzed separately. Landmarks are homologous site of geometric morphology with biological significance on the specimen, which were selected to reflect the shape of the mandible (Bookstein 1991). Landmarks and semi-landmarks were digitized in tpsDig 2.30 (Rohlf 2017), and their location is shown in Suppl. material 2: Figure S1. The raw datasets for each of the above four views were examined to evaluate whether the specimens greatly deviated from the average values. Generalized procrustes analysis (GPA) was used to carry out the superimposition of landmarks, in order to remove shape-irrelevant variables like size, orientation and position from the original landmark configurations, leaving the real shape information. This is a necessary step in geometric morphometric analysis (Gower 1975, Rohlf and Slice 1990). Relative distortion analysis was performed with tpsRelw, and the relative warp score (RW) was determined. Principal component analysis was conducted to visualize shape differences between individuals, and thin-plate spline transformation grids of extreme value were used to show skull shape differences (Bookstein 1997; Slice 2007).

We also analyzed the dorsal hair color, spiny-ness of hairs, yellow patches, and white tail tip of N. c. sacer, N. confucianus and N. lotipes to compare external morphological features by Chi-square test. The dorsal hairs color and spiny-ness of hairs were identified by observing and touching on pelage of specimens on three-category records, as the color of dorsal hairs were all brown, tan and all yellow; the spiny-ness of hairs were hard, medium and soft; the yellow patches on the chest observed by direct observations of pelage specimen with confirmation of specimen photos, which were recorded by dichotomy, yes or no; and the white tail tip was calculated based on the ratio of the measured tail tip length to the total tail length, divided into 4 ranks: 0, 1/4, 1/3, 1/2.

We obtained 1140 bp of mitochondrial Cytb sequences from 312 individuals with 166 haplotypes in this study. Among them, there were 708 conserved loci, 395 parsimoniously informative loci, and 37 single-variant loci; no insertion, deletion, or termination codons were found. The transition/transversion bias was 5.40, and nucleotides in all sequences were accurately translated into amino acids.

The phylogenetic trees constructed based on haplotype data using the NJ, ML, and BI methods showed essentially the same topology with high confidence values (Fig. 2, Suppl. material 2: Figs S2, S3): (1) the N. confucianus individuals clustered into three well-supported clades (i.e., clades C1, C2, C3). Clade C1 is distributed in southwestern Yunnan and southeastern Tibet; Clade C2 is found from central to northern China, and includes two haplotypes of Linyi and Zibo from Shandong; Clade C3 extended from the north of Southeast Asia to central China. (2) N. c. sacer split into two deep subclades in the central Shandong and Yantai, Weihai regions; interestingly, N. bukit from Vietnam and two haplotypes from Xishuangbanna in Yunnan clustered into one clade, forming a sister clade to N. c. sacer.

Figure 2. 

Phylogenetic analyses of Cytb gene from all haplotypes by maximum likelihood.

According to the constructed phylogenetic relationship, the K2P genetic distances within and between each clade were calculated. Genetic distances were found to range from 0.011 to 0.022 within the four clades and 0.053 to 0.084 between the four clades (N. confucianus, N. sacer, N. lotipes, and N. bukit). Among them, the genetic distance between N. c. sacer and N. bukit showed the lowest value (0.053); however, the genetic distance between N. c. sacer and other species was higher than 0.058 (Table 1).

Five ABGD analyses with different relative gap width values (X = 0.5, 1, 1.5, 2, and 2.5) were performed on 166 Cytb sequences, and two gaps (distance = 0.05 and 0.11) were observed (Suppl. material 2: Fig. S4). All analyses consistently supported a 17-group scenario when intraspecific divergence(p) = 0.001-0.0083 (Suppl. material 2: Table S2).

The species tree constructed by ABGD based on genetic distance is shown in Figure 3. The N. confucianus specimens are included in seven groups (groups 11–17), the haplotypes of N. c. sacer appeared as a monophyletic group (group 12), and the haplotypes of N. confucianus from Xishuangbanna in the Yunnan (group 16) and N. bukit from Vietnam (group 4) corresponded to N. bukit clade in the phylogenetic tree.

Figure 3. 

Comparison of ABGD species tree based on genetic distance and phylogenetic tree.

Karyotype analysis showed that the karyotype of N. c. sacer (♀2, ♂4) differed from that of N. confucianus (♂3). The diploid number (2n) is 46 for both, but the karyotype characteristics of N. c. sacer is FN = 55, 8m+4st+32t+X(sm)Y(t), chromosome composition: four pairs with metacentric chromosomes, two pairs with subtelocentric chromosomes, 16 pairs telocentric chromosomes, and two sex chromosomes; in contrast, the karyotype of N. confucianus is FN = 59, 6m+4sm+2st+32t+X(sm)Y(t), chromosome composition: three pairs with metacentric chromosomes, two pairs with submetacentric chromosomes, one pair with subtelocentric chromosomes, 16 pairs telocentric chromosomes, and two sex chromosomes (Fig. 4).

Figure 4. 

Karyotype of N. sacer (S3307) and N. confucianus (S1101).

A total of 98 adult individuals of N. confucianus, N. c. sacer and N. lotipes was screened by age identification, and the complete external indices of 84 individuals and skull indices of 72 individuals were obtained. Most characteristics showed normal distributions (P > 0.05), and thus we performed parametric statistics analysis (Suppl. material 2: Table S3). Information on general variation in body form of the four species is given in Table 2 and Suppl. material 2: Tables S4 and S5. The ANOVA results showed significant morphological differences between the species in external indices and most skull indices (Suppl. material 2: Table S4). Multiple comparison results showed the skull indices (LIF and LAB) of N. c. sacer are significantly larger than those of other three species (Table 2, Suppl. material 2: Table S5).

In the PCA of external indices, two factors had eigenvalues exceeding 1.0, and the first two axes captured 36.6% and 27.6% of the total variation, accordingly (Table 3). The TL and EL provided the greatest contribution to the factor loadings of PC1(0.764 and 0.619, respectively; Table 3). PC2 was mainly influenced by HBL (-0.711; Table 3). The main scatter plots showed that the specimens from the four species had greatly mixed external indices, indicating they cannot be distinguished by the external measurement method alone (Fig. 5a). In the PCA of skull indices, the first two axes captured 54.1% and 16.4% of the total variation, accordingly. LS, LD and LIF are the three measurements that have the highest correlation with PC1(0.936, 0.875 and 0.874, respectively; Table 3), which is obviously associated with size, as all the measurements have the same sign, and most have the same magnitude. PC2 was mainly influenced by LAB and BR (0.873 and -0.449, respectively; Table 3), which seems to be a shape factor given the different magnitudes and signs. The changes of length of auditory bulla and breath of rostrum are associated with PC2. The main scatter plots showed that the four species mix with each other, but part of N. c. sacer can separated from those of others (Fig. 5b, Table 3). In discriminant analysis, 68.3–80.3% of the individuals were correctly classified (Suppl. material 2: Table S6). The scatter plots of the discriminant function showed that individuals between clades were more likely to be confused based on external indices; however, based on skull indices, most individuals of N. c. sacer were accurately identified and classified, whereas a few individuals were easily confused with N. confucianus but completely separated from N. bukit and N. lotipes (Fig. 5c, d).

Figure 5. 

Principal component analysis and discriminant analysis of external and skull morphological indices. Principal component plots of external and skull indices are shown in a and b. Discriminant function plots of external and skull indices are shown in c and d, respectively.

Cluster analysis of external and skull measurement indices showed that the distributions of N. c. sacer and other three species are mixed and mosaic in the dendrogram (Suppl. material 2: Fig. S5), indicating that N. c. sacer and its close relatives cannot be distinguished based on traditional morphological characteristics.

The average configuration of the superimposition on the dorsal, ventral, lateral view of the skull and lateral view of the mandible of the three species is shown in Suppl. material 2: Figure S6. The main scatter plots constructed using RW1 and RW2 revealed no clear line between the three clades, indicating that variation in the samples from the three clades is not obvious, and geometric morphometric cannot distinguish the samples in the three clades (Fig. 6).

Figure 6. 

Principal component analysis of dorsal view (a), ventral view (b), lateral view (c) of skull, and lateral view of the mandible (d) of the three clades.

The thin-plate spline transformation grids of extreme value showed that there are some variations between the three clades (highlighted in red boxes): the maxilla of N. c. sacer is slightly wider than those of N. confucianus and N. lotipes in the dorsal view of the skull (Fig. 7a). The auditory vesicle of N. c. sacer is slightly larger than those of N. confucianus and N. lotipes in the ventral view of skull (Fig. 7b). The height of basion of N. confucianus is slightly higher than those of N. c. sacer and N. lotipes, and the top of the skull tended to be rounder in the lateral view of the skull (Fig. 7c). Samples from N. c. sacer showed a narrower coronal process of the mandibular teeth in the lateral view of the mandible (Fig. 7d). However, each deformation was not obvious enough to distinguish samples between the three species.

Figure 7. 

Thin plate splines of dorsal view (a), ventral view (b), lateral view (c) of skull, and lateral view of the mandible (d) of N. sacer, N. confucianus, and N. lotipes.

In a comparison of the external morphological features of N. c. sacer, N. confucianus, and N. lotipes, we found that samples in the three clades could not be distinguished based on the dorsal hair color and yellow patches, but there were significant differences (P < 0.05) in the spiny hairs and white tail tip: the spiny hairs of N. c. sacer is softer than those in the other clades; the tail color of N. c. sacer is the upper brownish black, the lower white, whereas the proportion of the white tail tip is more than 1/3, which is the same as the holotype specimens first found by Thomas (1908) in Yantai, Shandong. The ventral surface of the tail in N. confucianus is brownish black, the proportion of the white tail tip is less than 1/3, and the tail is more often without white hairs in N. lotipes, with a proportion of white hair of less than 1/4 (Table 4, Suppl. material 2: Fig. S7).

We analyzed specimens from across China and surrounding countries and recover evidence that Niviventer confucianus sacer should be elevated to Niviventer sacer. Molecular phylogenetic analysis indicated that N. sacer formed a sister branch with N. confucianus from Yunnan and N. bukit from Vietnam rather than with N. confucianus from adjacent areas in Shandong (Shanxi, Jiangsu, Hebei, Henan, etc.), which is consistent with the results of Ge et al. (2018a, b). Niviventer sacer is also distinct according to the species delimitation (ABGD method), and the genetic distance (K2P) between N. sacer and N. confucianus was as high as 5.8%, whereas that with N. bukit was 5.3%, both exceeding the empirical threshold of dividing rodent species (5%, Avise and Walker 1999; Baker and Bradley 2006; Zhang et al. 2019), revealing that N. sacer is a distinct species.

In contrast to the conclusions of Ge et al. (2018a, b), who suggested that N. c. sacer may be a population of N. bukit, we consider that they are two independent species based on the following four reasons:

(1) in terms of geographical distribution, N. sacer and N. bukit are distributed in Shandong and Vietnam, respectively, separated by a distance of more than 2,000 kilometers;

(2) Zhang et al. (2016) used nuclear genes to construct a multilocus phylogeny tree, indicating that N. c. sacer was nested among several geographically-adjacent subspecies of N. confucianus to the exclusion of N. bukit;

(3) Ge et al. (2018b) only used Cytb and the number of specimens of N. c. sacer was limited, and for this reason N. c. sacer was clustered with N. bukit rather than N. confucianus, it may be the misleading effect of rapid or saturation mutation of Cytb;

(4) Morphological analysis shows there are significant differences between N. sacer and N. bukit in external and skull indices.

Interestingly, N. sacer distributed in the east hills and southwest mountains in Shandong is divided into two independent small lineages. However, N. confucianus is only distributed in the southwest mountains (Mount Lu and Mount Meng) in the Shandong area and forms the same clade with N. confucianus from adjacent areas. Therefore, N. sacer and N. confucianus show sympatry characteristics in the southwest mountains area, but their genetic distance (5.6%) is higher than 5% and they have significantly different karyotypes and morphologies. The sympatry of N. sacer and N. confucianus in the southwest mountains indicates that they are distinct species rather than subspecies without hybridization.

Morphological analysis showed that N. sacer, N. bukit, N. confucianus, and N. lotipes had significantly different morphological characteristics, which are reflected in the larger skull, but have similar skull shapes and characteristics. The morphological characteristics of the tail are important traits for distinguishing different species and are probably associated with adaptations for an arboreal lifestyle in different forest types. Moreover, tails may play important roles in the recognition of conspecifics (Siegel 1970, Ge et al. 2018b). In a comparison of the morphological characteristics of N. sacer and N. confucianus, we found a significant difference in tail color: the upper tail color of N. sacer is brownish black, while the lower color is white, and approximately one-third of the tip was all white, which has been observed in holotype specimens found by Thomas (1908) in Yantai, Shandong; in contrast, the tail of most N. confucianus are brown-black with only a white tail tip.

The karyotype of N. sacer in this study is 2n = 46, FN = 55, 8m+4st+32t+X(sm)Y(t), which is consistent with that of N. confucianus from Shandong as described by Wang et al. (1997, 2003). The karyotype of N. confucianus in this study is 2n = 46, FN = 59, 6m+4sm+2st+32t+X(sm)Y(t), indicating that N. sacer differs from N. confucianus in karyotype. Cytogenetic evidence also supports N. sacer as a distinct and valid species.

The phylogenetic tree results showed that N. confucianus species complex is mainly divided into four clades. The first clade is N. confucianus, which is found in central China, extending from the northeast to southwest of China; The second clade is N. sacer, which is endemic to Shandong; The third clade is N. lotipes, which is distributed in the southeast of China; The four clade is distributed in southwestern Yunnan and southeastern Tibet, which has been considered as a new combination, N. pianmaensis Ge et al. 2018b. The range of the four species are consistent with Ge et al. (2018a, b). Ge et al. (2018a, b) conducted molecular phylogeny analysis of a species complex of common wild rat species in China and observed the historical dynamics of N. confucianus based on coalescence models. It was predicted that N. c. sacer was not a subspecies of N. confucianus, but rather N. bukit, and no N. confucianus was found in Shandong. The present study demonstrates that N. sacer is a valid species distributed only in Shandong.

The minimum genetic distance (K2P) between the four clades in this study was 0.053. The results of ABGD species delimitation showed that N. confucianus are divided into seven groups which are consistent with the phylogenetic tree.

In this study, N. sacer was found to be a relatively recent divergence from N. confucianus, which differs from the results of Zhang et al. (2016) and Ge et al. (2018a, b). We used more specimens from the east hills branch and southwest mountains branch of N. sacer, which also indicated that different datasets had significant effects on the systematic evolution relationship analysis within the genus Niviventer.

In addition, the phylogenetic findings in this study are similar to those of Lu et al. (2015) and Zhang et al. (2016), who observed paraphyly among N. andersoni and N. excelsior, and N. fulvescens and N. huang. Niviventer fulvescens and N. huang form a clade in the phylogenetic tree; genetic distance within the clade is 0.019 and these species formed a group according to the ABGD method. Therefore, N. huang may require synonymization with N. fulvescens. The relationships and taxonomic status of these species require further investigation.