A total of 197 Myotis specimens covering 21 species were analyzed for 5 external and 13 craniodental character measurements according to the Chiroptera Morphometric Standard from Yang et al. (2007) using MNT-150 vernier calipers (Meinert Industrial Co., Ltd., China, accuracy 0.01 mm). External measurements (Table S2) include: FA forearm length; TIB tibia length; HF hindfoot length including claws; EH ear length and TL tail length. Craniodental measurements (Table S3) include: GTL greatest length of skull; CCL condylo-canine length; ZB zygomatic breadth; MAW mastoid width; BCW braincase width; BCH height of braincase; PBL palatal bridge length; C¹M³L maxillary toothrow length; C¹C¹W upper canine width; M³M³W width across the upper molars; C₁M₃L mandibular toothrow length; ML mandible length including incisors and MH mandibular height. In addition, relevant morphological data were supplemented by review of the literature to better reflect morphological differences in Myotis species.

Although morphometric approaches effectively quantify overall body size differentiation—an important criterion in species delimitation—they may disregard discrete morphological characters with considerable diagnostic utility. Acknowledging that the rigorous coding and analysis of such qualitative traits significantly enhances species recognition accuracy, we conducted a comprehensive review of taxonomic literature on Myotis from China and neighboring regions (Smith and Xie 2009; Kruskop 2013; Ruedi et al. 2015, 2021; Wilson and Mittermeier 2019). Based on this evaluation, we selected 9 stable morphological characters (defined as traits exhibiting minimal intraspecific variation) for systematic categorization and coding. The characters and rules are as follows, VH ventral hair color – 1 to 5: nearly black, tan, gray or off-white, golden, orange-red; DH dorsal hair color – 1 to 5: nearly black, tan, grayish brown, golden, orange-red; WM location of wing membrane attachment – attached near the base of the toes vs attached at the ankle; SU structure of the uropatagium – the tip of tail is free from membrane, vs tail is totally included into uropatagium; SC sagittal crest – distinct vs inconspicuous or absent; LC lambdoid crest – distinct vs inconspicuous or absent; PB profile of braincase – flattened vs arched; P² location of the first upper premolar – in the dental formula vs intruding laterally from the dental formula; P³ location of the second upper premolar – in the dental formula vs intruding laterally from the dental formula. The coding rules diagram and coding results for each species could be found in Figure S1 and Table S4.

The maximum likelihood (ML) method was employed for phylogenetic reconstruction based on the mitochondrial cyt b gene and combined cyt b-Rag2 genes. Additionally, preliminary analyses revealed that the Rag2 gene evolves much slower than the cyt b gene, resulting in insufficient phylogenetic signal. Therefore, we did not analyze Rag2 independently and retained it solely as a corroborative nuclear marker in the concatenated matrix.

Maximum likelihood method is sensitive to nucleotide substitution models, we selected the optimal models using Bayesian Information Criterion (BIC) in ModelFinder (Kalyaanamoorthy et al. 2017). The optimal model for both gene datasets was TIM2+F+I+R3. Subsequently, maximum likelihood phylogenetic trees were constructed based on the above model using IQ-TREE v.2 (Minh et al. 2020), with Kerivoula furva as the outgroup. The branching support was obtained after 1000 non-parametric bootstrap resampling. To assess divergent levels among species, interspecific genetic distances were calculated using the Kimura 2-parameter (K2P) model in MEGA v.11 (Tamura et al. 2021). In addition, substitution saturation in the cyt b datasets was evaluated using DAMBE v.7 (Xia et al. 2018) to confirm the reliability of phylogenetic analyses.

Species delimitation analyses were conducted using the SPdel pipeline (Ramirez et al. 2023), which is designed to integrate multiple species delimitation methods and generate consensus molecular operational taxonomic units (MOTUs). Following its default workflow, we applied 5 methods to both the cyt b dataset and the combined dataset: automatic partitioning (ASAP, Puillandre et al. 2021), generalized mixed yule coalescent (GMYC, Monaghan et al. 2009), Poisson tree process (PTP, Zhang et al. 2013), Bayesian Poisson tree process (bPTP) and multi-rate Poisson tree process (mPTP, Kapli et al. 2017). The pipeline requires only the corresponding sequence alignment and the original maximum likelihood tree file as input, and automatically handles all method-specific processing requirements. After obtaining delimitation results from the 5 methods, consensus MOTUs were generated based on majoritily consistent units: A group is retained as the consensus MOTU if it is supported by more than half of the methods. All analyses were performed using default parameters as implemented in the SPdel pipeline (for details, see: https://github.com/jolobito/SPdel).

To investigate the distinctions and potential benefits between metric-only and combination of metric and character-coded data in species classification, two distinct data matrices were meticulously constructed. Matrix 1 comprises 18 external and craniodental metric measurements, whereas Matrix 2 integrates both 9 phenotypic encoded data and the complete set of metric measurements. All specimens had complete measurements and character codings, with no missing data present in either matrix. Therefore, no imputation was required, and the full dataset was retained for all subsequent analyses.

Initially, Principal Component Analyses (PCA) were conducted independently on the external and craniodental measurements, as well as on Matrix 1. The first two principal components were extracted and visualized as scatterplots. Morphometric pairwise distances were subsequently estimated by calculating the Euclidean distance between the centroids of each species in the PCA scatterplot, based on Matrix 1. Additionally, the correlation between interspecific K2P distance (derived from cyt b gene) and morphometric distance was analyzed. To reconstruct similarity relationships from a morphometric standpoint, hierarchical clustering (HC) was employed to generate the morphometric dendrogram of Myotis.

While the above methods effectively visualize morphometric relationships, limitations exist in identifying optimal trait combinations for species identification. To address this challenge, we incorporated decision tree and random forest analyses. The former provides an interpretable framework for feature exploration and generates explicit classification pathways that can directly assist in constructing identification keys; the latter captures complex trait interactions by aggregating numerous decision trees, overcoming the reliance of single decision trees on single-node thresholds and achieving higher predictive accuracy. For both algorithms, the response variable was species identity (21 Myotis species), and the predictor variables were the two data matrices described above.

In the decision tree model, to mitigate the issue of small-sample bias, which often leads to classification models favoring the majority class and neglecting the minority class as noise (He and Garcia 2009; Buda et al. 2018; Chen et al. 2024), the synthetic minority over-sampling technique (SMOTE, Chawla et al. 2002) was utilized to improve the data distribution pattern and enhance the model generalization. Model analysis was conducted on both the original dataset and the SMOTE-augmented dataset. For the random forest model, we used 10-fold cross-validation with 200 trees per forest to obtain robust classification accuracy estimates while maintaining computational efficiency. The same tree number was applied across all folds to ensure reproducibility. Subsequently, feature importance was assessed by the Mean Decrease in Gini (MDG), which was then used to identify trait combinations with high classification value. Finally, we integrated the classification pathways generated by decision trees with the key features identified by random forests to update the identification key for Chinese Myotis.

All of the above morphological analyses were implemented using the R packages: caret (Kuhn 2008), cluster (Maechler et al. 2021), DMwR (Torgo 2010), factoextra (Kassambara and Mundt 2020), FactoMineR (Le et al. 2008), ggplot2 (Wickham 2016), ggpubr (Kassambara 2025), igraph (Csardi and Nepusz 2006), randomForest (Liaw and Wiener 2007) and rpart (Therneau and Atkinson 2023) in R v.4.3.2.

To facilitate future species determination and academic collaboration, we used a Rexcan DS3 Silver laser scanner (maximum resolution 0.01 mm; Solutionix, Korea) to construct high-resolution 3D digital models of 21 representative Myotis skulls (see Figs 5 and S2 for examples). Although micro-computed tomography (μCT) scanners would capture more accurate details with better resolution, the considerably smaller size of our files than those produced by a μCT scanner (e.g., ~40 MB vs. ~600 MB in the case of a skull scan), and sufficient accuracy suggest that laser 3D scanners can be used as an alternative for shape analyses and morphological studies (Yu et al. 2021). All 3D model files are publicly available in File S3. These digital resources were generated as independent data products of this study, aiming to establish a lasting foundation for subsequent taxonomic, ecological, and conservation research on this challenging genus.

K2P genetic distances based on the cyt b gene ranged from 8.3% to 21.9% among the 30 Myotis species (Fig. 2C). Except for three species pairs (M. alticraniatus vs. M. laniger, M. bombinus vs. M. pequinius, M. fimbriatus vs. M. pilosus), all interspecific distances exceeded 10%. Analysis of the combined cyt b-Rag2 dataset yielded largely consistent results (Fig. S5).

To quantify morphological differences, interspecific Euclidean distances of the PCA scatterplots for each species were calculated (Fig. 2D). It was found that M. altarium, M. blythii, M. chinensis, and M. pilosus were the most morphologically distinct, primarily due to size differences. In contrast, the remaining small and medium-bodied species exhibited high morphological similarity. Correlation analysis between the genetic distance matrix (constructed based on the cyt b gene) and the morphological distance matrix detected no significant relationship between morphological and genetic distance matrices among the 21 Myotis species (r = -0.073, p = 0.29; Fig. 2E). A substitution saturation test conducted in DAMBE confirmed that sequence substitutions were not saturated (p < 0.001), suggesting that morphological variation in Chinese Myotis is shaped by local ecological adaptation and convergence rather than genetic divergence alone.

The decision tree models trained on two separate matrices initially distinguished only 7 and 8 taxonomic groups, with overlap observed across all sizes of species (Fig. S6). After augmenting Matrix 2 (metric + encoded data) via SMOTE, model performance improved in a data-dependent manner. A threefold increase in data volume enabled the model to distinguish up to 16 taxonomic units (Fig. 3A); further increases did not improve this limit, suggesting that single-node decision thresholds are inadequate for capturing multi-trait synergies. We note, however, that synthetic data may not fully represent biological variations, potentially influencing model generalization.

Figure 3. 

A Decision tree classification model constructed on SMOTE-augmented data (3-fold). B Characteristicimportance analysis based on the Random Forest model (sorted by Mean Decrease Gini).

In random forest models, 10-fold cross-validation showed that incorporating categorical phenotypic data improved classification accuracy of Chinese Myotis. The model trained on Matrix 1 (metrics only) achieved a mean accuracy of 77.9%, with misclassifications spanning multiple size classes (e.g., M. chinensis as M. blythii; M. formosus as M. indochinensis or M. pequinius; M. laniger as M. muricola, etc.). The model using Matrix 2 (metrics + categorical phenotypic data) reached 90.5% accuracy, with most errors involving medium and small species such as M. formosus and M. laniger.

To identify key diagnostic traits, we evaluated feature importance using random forest models based on the Mean Decrease Gini. The most informative characters included FA, MAW, HF (incl. claws), VH, ZB, ML, CCL, TIB, DH and GTL (sorted by classification contribution, Fig. 3B). These comprise both stable qualitative traits (e.g., pelage color: VH, DH) and key metric thresholds (e.g., FA, CCL). Subsequently, by integrating the random forest importance analysis, the classification pathways revealed by decision trees, and a comprehensive review of taxonomic literature (Smith and Xie 2009; Kruskop 2013; Ruedi et al. 2015, 2021; Wilson and Mittermeier 2019), we compiled a set of diagnostically informative morphological characters (Table S5) and gathered taxonomic notes for poorly known species. Based on these results, we propose a comprehensively revised identification key for Chinese Myotis (see Appendix and File S2).

The taxonomy of Chinese Myotis taxa remains fraught with uncertainties, which hamper our understanding of their true diversity and evolutionary relationships. Based on integrative evidence, including molecular species delimitation, detailed morphological comparisons, and reference to the latest taxonomic revisions (Ruedi et al. 2015, 2021; Liu et al. 2023), we re-evaluated several species and species complexes to establish a foundation for taxonomic research within this genus. For the taxonomic treatment and revisions of these controversial species, please refer to Fig. S7 for a quick overview of their specific taxonomic transitions.

The taxonomic status of M. davidii (Peters, 1869) has been obscured by nomenclatural complexity and historical misidentification. Benda et al. (2016) synonymized M. nipalensis (Dobson, 1871) and M. aurascens Kuzyakin, 1935 with M. davidii, but Ruedi et al. (2021) reinstated M. nipalensis as a distinct species based on phylogenetic evidence and noted that many GenBank sequences labeled as M. davidii from Southern China actually represent M. alticraniatus (Osgood, 1932) or M. badius Tiunov, Kruskop & Feng, 2011. Our results support this revision: Southern Chinese populations formerly attributed to M. davidii are genetically divergent from the true M. davidii occurring in Russia and elsewhere (Fig. 1). Their morphology also aligns with the type specimens of M. alticraniatus and M. badius (Tiunov 2011; Ruedi et al. 2021; Figs 4, 5), particularly in the nyctalodont configuration of the first lower molar (Fig. S8A, B). The proposed elevation of M. aurascens to species rank by Yang et al. (2023) is also controversial. Their cited M. davidii sequence has been reassigned to M. alticraniatus (Ruedi et al. 2021; Liu et al. 2023), and their newly submitted M. aurascens sequence (OK053029) falls within a highly supported clade (BS = 100, Fig. 1) containing M. davidii sensu Benda et al. (2016), which was itself revised from M. aurascens. Whether M. davidii and M. aurascens are conspecific requires re-evaluation using topotype material of M. davidii from Beijing. Currently, it is evident that all sequences from southern China identified as M. davidii pertain to M. alticraniatus, which is distinct from M. aurascens found in Eastern Europe and West Asia. We therefore recommend treating M. badius as a junior synonym of M. alticraniatus and reassigning the southern Chinese populations as M. alticraniatus.

The M. siligorensis group presents particularly challenging taxonomic issues. Myotis alticraniatus was formerly considered a subspecies of M. siligorensis (Horsfield, 1855) but was elevated to species level by Ruedi et al. (2021) based on morphological and phylogenetic distinctions. Due to limited sampling and minimal subsequent research (Xiao et al. 2017; Chen et al. 2025; Ding et al. 2025), this revision has not been widely adopted in China, where many records of M. siligorensis likely represent misidentified M. alticraniatus. For instance, a sequence from China accessioned as M. siligorensis (FJ215679) clusters within the M. alticraniatus clade in our phylogeny (Fig. 1). We speculate that many Chinese populations reported as M. siligorensis may in fact be M. alticraniatus given the species’ wide geographical distribution in China. However, in the absence of true M. siligorensis specimens for direct comparison, we advise caution in using these contested sequences. Clarifying the taxonomic identity of these populations will require meticulous morphological and genetic analyses.

The M. frater complex also necessitates taxonomic re-evaluation. Northern populations (including M. longicaudatus Ognev, 1927, M. kaguyae Imaizumi, 1956 and M. eniseensis Tsytsulina & Strelkov, 2001) were long considered subspecies of M. frater Allen, 1923 until Ruedi et al. (2015) consolidated them into a single species distinct from true M. frater of southeastern China. Our phylogenetic analysis supports this taxonomic revision: Samples identified as M. frater (but likely M. longicaudatus) from northeastern China (Liu et al. 2023) form a clade distinct from true M. frater (Fig. 1), with a cyt b genetic distance of 13.5% (Fig. 2C). The species M. soror Ruedi, Csorba, Lin & Chou, 2015 also belongs to this complex, as described from Taiwan. Our specimens from Sichuan initially identified as M. frater were re-identified as M. soror following re-examination (Figs 4, 5). The two species are morphologically distinguishable by ear morphology: Myotis frater exhibits a less distinct notch on the posterior margin of the concha and multiple convex folds on its interior, whereas M. soror shows the opposite condition (Fig. S8C, D).

Mislabeled sequences in public databases pose serious obstacles to species identification in Myotis. For example, several sequences assigned to M. montivagus (Dobson, 1874) are likely erroneous, while Wilson and Mittermeier (2019) suggested that these sequences are closer to M. indochinensis Son et al., 2013. We provisionally refer to these as M. cf. montivagus (Figs 4, 5). In our phylogeny, this group forms a well-supported clade sister to the M. indochinensis clade (Fig. 1). Strikingly, the forearm length of these specimens (35.3–36.8 mm) is considerably shorter than that of true M. montivagus (39.2–41.5 mm) and also differs from that of M. indochinensis (43.7–45.6 mm). Resolving the identity of this entity will require integrative taxonomic investigations using verified reference specimens.

Finally, we provide new evidence regarding the occurrence of M. longipes (Dobson, 1873) in China. Liu et al. (2023) re-identified most historical Chinese records of M. longipes as M. laniger (Peters, 1870) following Ruedi et al. (2021), while Wei et al. (2025) consequently treated M. longipes as M. csorbai Topál, 1998 in their checklist. This indicates that the current consensus denies the existence of M. longipes in China. However, samples collected from Xizang in the present study are genetically affiliated with topotypic M. longipes from India (MW054878/79) and are clearly differentiated from M. laniger (including misidentified M. longipes; Fig. 1), representing the first confirmed record of M. longipes in China (Figs 4, 5). Morphologically, M. laniger and M. longipes differ in dentition: In M. laniger, the second upper premolar (P³) is situated within the toothrow, the length ratio of the first upper premolar (P²) to the upper canine (C¹) is approximately 1/2, and the area of P³ is about one-third that of C¹ (Fig. S8E). In M. longipes, P³ is displaced laterally from the toothrow, the P²/C¹ ratio is approximately 1/3, and the area of P³ is only about one-quarter that of C¹ (Fig. S8F). Notably, Ruedi et al. (2021) suggested that M. longipes and M. csorbai are conspecific, contradicting those from Wei et al. (2025). We thus retain M. longipes as a valid species present in China.

Against the backdrop of the rapid global decline in biodiversity, establishing an accurate and practical taxonomic framework is crucial for species recognition and effective conservation (Khater et al. 2021; Samayoa et al. 2022; Sol et al. 2023; Hu et al. 2025). However, the revision and refinement of the Myotis taxonomic system in China has long relied on relatively singular species delimitation works, such as isolated reports of new species or distribution records (Xiao et al. 2017; Yang et al. 2023; Ding et al. 2025), lacking systematic integration. Meanwhile, existing identification keys (e.g., Smith and Xie 2009; Kruskop 2013) suffer from limitations including incomplete geographic coverage, insufficient specimens of key taxa, and inconsistent data sources. Consequently, they fail to reflect current species diversity and incorporate recent taxonomic revisions, making it imperative to develop a new taxonomic framework and identification system grounded in broader sampling and multidimensional evidence.

To achieve accurate identification and construct a reliable taxonomic tool for Chinese Myotis species, this study employed multiple statistical classification methods (PCA, HC, DT, and RF) to systematically analyze morphological trait variation (Fig. 3B). In the light of key taxonomic references (e.g., Smith and Xie 2009; Ruedi et al. 2015, 2021; Wilson and Mittermeier 2019), we filtered a suite of diagnostic characteristics (Table S5) and subsequently compiled an updated identification key encompassing all 30 Myotis species in China (Table 1). This key incorporates three major refinements: 1) emphasis on readily observable and measurable external traits to support preliminary field identification; 2) preferential use of encoded traits that remain stable across individuals, particularly for rarely collected species (e.g., M. formosus, M. frater, etc.) to minimize the effect of limited sample size; and 3) enhanced discrimination of morphologically conserved groups through the inclusion of subtle craniodental characters (Figs 5, S1, S8). The key not only reflects the most current understanding of Myotis diversity in China but also offers a practical tool for future surveys, taxonomic research and conservation actions.

It is important to note that the proposed taxonomic framework may be susceptible to overfitting, given the high character-to-specimen ratio. We mitigated this risk by employing an ensemble strategy based on the random forest algorithm and 10-fold cross-validation. However, potential biases may persist due to the use of synthetic SMOTE-augmented data, which may not fully represent true morphological variation. Concurrently, this study primarily focused on 21 Myotis species from eastern China (Figs 4, 5). For the remaining 9 species from western regions, their morphological descriptions relied primarily on literature records. Therefore, we provided two versions of the identification key: one encompassing all 30 species (see below for details), and another specific to the 21 eastern species (File S2). We suggest that future studies should expand specimen sampling across species and geography to improve the key’s robustness and generalizability. We anticipate that the taxonomic framework (Table 1) and open data resources presented here will support further research and conservation efforts.