Presumptive First Record of Myotis aurascens (Chiroptera, Vespertilionidae) from China with a Phylogenetic Analysis

Simple Summary The taxonomic status of Myotis aurascens has been debated. Based on morphological and molecular data, our study showed that M. aurascens was closely related to M. ikonnikovi, M. alcathoe, and M. mystacinus and was not a synonym for M. davidii. Therefore, the integrated analysis demonstrated that M. aurascens should be considered a distinct species rather than a synonym of M. davidii. This study clarified the taxonomic status of M. aurascens, which should help enrich biodiversity and improve the species determination of bats in China. Abstract Bat groups have a high degree of species diversity, and the taxonomic status and phylogenetic relationships among bat species have always been research hotspots. Due to the fact that morphological characteristics do not always reflect the evolutionary relationships among species, mitochondrial DNA has been widely used in the study of species relationships due to its maternal inheritance pattern. Myotis aurascens has been suggested as a possible synonym for M. davidii. However, the status of this classification has been controversial. In this study, the morphological and molecular characteristics of a M. aurascens captured from Inner Mongolia, China, were analyzed to determine its taxonomic status. In terms of morphological features, the body weight was 6.33 g, the head and body length were 45.10 mm, the forearm length was 35.87 mm, and the tragus length was 7.51 mm. These values all fell within the species signature data range. Nucleotide skew analysis of the protein-coding genes (PCGs) suggested that only five PCGs (ND1, ND2, COX2, ATP8, and ND4) showed AT-skew value within the mitogenome of M. aurascens. Except for ND6, the GC-skew values of the other PCGs were negative, reflecting the preference for C and T bases compared to G and A bases. Molecular phylogenetic analyses based on mitochondrial PCGs indicated that M. aurascens was a distinct species from M. davidii and phylogenetically closer to M. ikonnikovi, M. alcathoe, and M. mystacinus. Genetic distance analysis also showed that M. aurascens and M. davidii were distantly related. Therefore, the integrated analysis demonstrated that M. aurascens should be considered a distinct species rather than a synonym of M. davidii. Our study could provide a reference for enriching species diversity and research on conservation in China.


Introduction
Myotis is one of the most diverse genera in Mammalia and includes about 90 species worldwide. Myotis species are widely distributed all over the world and are found on all continents except Antarctica [1,2]. Currently, 27 species of the genus and their distributions have been documented in China [3]. Most species in this genus have been well recognized on the basis of morphological and phylogenetic analyses. Nevertheless, the taxonomic statuses of several congeners remain unresolved, for instance, Myotis aurascens A single individual of M. aurascens was collected from Hulun Lake National Nature Reserve, Inner Mongolia, China, on 1 August 2018. The sampled individual was captured using mist nets during animal investigations. The morphological indexes were measured using an electronic balance (accurate to 0.01 g) and vernier calipers (accurate to 0.01 mm). The sample was frozen in an ultra-low temperature freezer. DNA was extracted from the tissue using the DNeasy Blood & Tissue Kit (QIAGEN, Beverly, MA, USA). All sample procedures and experimental methods were approved by the Qufu Normal University Institutional Animal Care and Use Committee (Permit Number: 2022002).

MtDNA Sequencing, Assembly, and Annotation
The library was constructed and sequenced using an Illumina MiSeq platform (Illumina, San Diego, CA, USA) with 150 paired ends. The raw reads were filtered and quality controlled to obtain clean reads for subsequent mitochondrial genome assembly. First, A5-MISeq (v. 20150522) [15] and SPAdes (v. 3.9.0) [16] were used to assemble the clean data into contig and scaffold sequences. Then, Blastn (v. 2.2.31; https://blast.ncbi.nlm.nih.gov/Blast. cgi?PROGRAM=tblastn&PAGE_TYPE=BlastSearch&BLAST_SPEC=&LINK_LOC=blasttab& LAST_PAGE=blastn, accessed on 15 May 2022) was used to compare sequences of high sequencing depth with the Nt database, and the mitochondrial sequences of each spliced result were selected. Finally, Pilon (v. 1.18) [17] software was used to correct the results to obtain the final mitochondrial sequence. After annotating the sequence using the online software Banklt, the genome was deposited in GenBank with the accession number OK053029.

Phylogenetic Analysis
Phylogenetic analysis was performed by comparing the mitogenome sequences of M. aurascens with those of other Vespertilionidae (including 43 Myotis, 1 Hypsugo, 5 Murina, 3 Nyctalus, 3 Pipistrellus, 3 Plecotus, and 2 Vespertilio), Miniopterus fuliginosus, and Tadarida latouchei as outgroups of the phylogeny. The models of evolution were evaluated using corrected Akaike information criteria (AICc) in ModelTest (v. 3.7) [20] to determine the best nucleotide substitution model [21]. Two methods, RAxML [22] and MrBayes [23], were used to construct a maximum likelihood (ML) tree based on three PCGs (ND1, Cytb, and COX1). The RAxML tree was evaluated with a bootstrap test with 1000 generations, and the Markov chain Monte Carlo (MCMC) analyses of the Bayes tree were run for 1,000,000 generations. The PAUP [24] software was used to construct an ML tree based on 12 PCGs. StarBeast (v. 2.6.7) software [25] was used to construct phylogenetic trees for species classification to prove our viewpoint. P-distance and the maximum composite likelihood method were used to calculate the genetic distance of nine species with 1000 bootstrap replications in MEGA7 [26].

Habitat and Morphology
The sample of M. aurascens was collected from Hulun Lake National Nature Reserve (48 • 54 28.41" N, 117 • 5 2.65" E), and the habitat type was lakeshore cliffs (Figure 1). The measurement results of the external morphology showed that the body mass (BM) was 6.33 g, the forearm length (FL) was 35.87 mm, the length of the head and body (LHB) was 45.10 mm, the ear length (EL) was 13.22 mm, and the ear width (EW) was 7.25 mm (Table 1). Compared with the report of Kim [6], the tibia length (TBL), tail length (TL), and tragus length (TRL) of the Inner Mongolian specimen were shorter, at 15.11 mm, 32.70 mm, and 7.51 mm, respectively. The TL of the Chinese population was also shorter than those in Oh's study (36.00-42.00 mm) [27]. Compared with M. davidii, M. aurascens is larger in BM, FL, LHB, EL, and EW values. Compared to M. davidii, the M. aurascens specimen was larger. PCA analysis results based on morphological data (FL, LHB, TL, and EL) supported that our study sample was M. aurascens ( Figure S1).

Genome Organization
The assembled sequence was deposited at NCBI in our previous study with the accession number OK053029 [14]. The sequence of M. davidii was downloaded from NCBI under accession number NC_025568.1. As shown in Figure 2, the mitogenome of M. aurascens was a circular DNA molecule of 16,771 bp in length, and that of M. davidii was 17,531 bp. Both genomes contain 37 genes, including 13 PCGs, 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs), and a D-loop. The length difference between the two species was mainly in the rRNA genes and D-loop. Most of the genes were encoded in the heavy strand (H-strand), while eight tRNAs and ND6 were encoded in the light strand (L-strand), which was similar to that of other species of the same genus [29].

Genome Organization
The assembled sequence was deposited at NCBI in our previous study with the accession number OK053029 [14]. The sequence of M. davidii was downloaded from NCBI under accession number NC_025568.1. As shown in Figure 2, the mitogenome of M. aurascens was a circular DNA molecule of 16,771 bp in length, and that of M. davidii was 17,531 bp. Both genomes contain 37 genes, including 13 PCGs, 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs), and a D-loop. The length difference between the two species was mainly in the rRNA genes and D-loop. Most of the genes were encoded in the heavy strand (H-strand), while eight tRNAs and ND6 were encoded in the light strand (L-strand), which was similar to that of other species of the same genus [29].
The nucleotide composition of M. aurascens is shown in Table 2. The percentage of A + T (64.8%) was significantly higher than that of C + G (35.2%), which was considered to play a vital role in the transcription and replication of the mitochondrial genome [30]. The percentages of A + T and C + G were 63.64% and 36.33%, respectively. The analysis of the complete mitogenome nucleotide skew showed that it had a positive AT-skew value (0.046) and a negative GC-skew value (−0.262) (Table S1). The nucleotide composition of M. aurascens is shown in Table 2. The percentage of A + T (64.8%) was significantly higher than that of C + G (35.2%), which was considered to play a vital role in the transcription and replication of the mitochondrial genome [30]. The percentages of A + T and C + G were 63.64% and 36.33%, respectively. The analysis of the complete mitogenome nucleotide skew showed that it had a positive AT-skew value (0.046) and a negative GC-skew value (−0.262) (Table S1).

Characterization of Coding Genes
The total length of the 13 PCGs was 11,405 bp, which accounted for 68.0% of the whole mitogenome (Table 2). Similar to the whole genome, the average A + T content of PCGs in the mitogenome was 65.92%, ranging from 62.46% (COX1) to 70.41% (ATP8), and was higher than G + C (34.08%) in 13 PCGs. Only five PCGs had positive AT-skew values in the mitogenome of M. aurascens (ND1, ND2, COX2, ATP8, and ND4), while the GC-skew value was negative except for ND6, reflecting the preference for C and T bases compared to G and A bases ( Table S1). All of the initiation codons of the PCGs were ATG or ATA codons, except for the ND2 gene, which started with ATT. All of the termination codons were TAA or truncated T residues, except for the Cytb gene, which stopped with AGA ( Table 2).
The count and RSCU values of M. aurascens are shown in Table S2 and Figure 3. Among them, 26 codons were used more frequently (RSCU > 1, Table S2), and the most commonly used codons were AUU-Ile (249), CUA-Leu (236), AUA-Ile (231), UUA-Leu (190), and CUG-Phe (159). The entire length of the 22 typical tRNA genes was 1514 bp and ranged from 59 bp (Ser1) to 74 bp (Leu2). Eight of these were located on the L-strand and fourteen on the H-strand, and all of these tRNA genes were able to form the classical cloverleaf secondary structure except for seryl-tRNA (Table 2 and Figure S2).

Phylogenetic Analysis within Myotis
Firstly, we determined if our sample was M. aurascens and if the sequence we downloaded was M. davidii. We downloaded the published sequences of the COI, Cytb, and ND1 genes of M. aurascens and the Cytb sequence of M. davidii from NCBI. The phylogenetic tree results confirmed that the sequence we selected was correct (Figures 4 and S3). To explore the taxonomic status of M. aurascens and its phylogenetic relationship with M. davidii and M. mystacinus, phylogenetic trees were constructed using all or some of the PCGs, respectively. Based on 12 PCGs (excluding ND6) of 45 species (Table S3), a phylogenetic tree was obtained using the ML method with 1000 replications, in which T. latouchei and M. fuliginosus were set as the outgroups ( Figure 5). As a result, there were two main clades within the genus Myotis (blue and rose); compared with M. davidii, M.

Phylogenetic Analysis within Myotis
Firstly, we determined if our sample was M. aurascens and if the sequence we downloaded was M. davidii. We downloaded the published sequences of the COI, Cytb, and ND1 genes of M. aurascens and the Cytb sequence of M. davidii from NCBI. The phylogenetic tree results confirmed that the sequence we selected was correct (Figures 4 and S3). To explore the taxonomic status of M. aurascens and its phylogenetic relationship with M. davidii and M. mystacinus, phylogenetic trees were constructed using all or some of the PCGs, respectively. Based on 12 PCGs (excluding ND6) of 45 species (Table S3) Since the mitochondrial genome of M. mystacinus was not available, phylogenetic trees of 64 species were constructed based on single (Table S4, Figure 4) and combined gene sequences of ND1, Cytb, and COX1.
Based on COX1 gene sequences, the phylogenetic tree showed that M. aurascens was closely related to M. nipalensis. Based on the combined sequences of the three genes, the two methods obtained different topological structures: the RAxML tree showed M. aurascens and M. ikonnikovi to be sister species and then clustered with M. alcathoe and M. mystacinus with moderate bootstrap support ( Figure S4). On the contrary, the Bayes tree showed that M. aurascens first clustered with M. alcathoe and M. mystacinus, and then clustered with M. ikonnikovi with strong bootstrap support. However, the Bayes tree obtained high support values, so the result was more convincing (Figure 6). Both topological structures indicated that M. mystacinus and M. davidii were distantly related and should belong to different species. Further, StarBeast was used to construct phylogenetic trees for species classification. Several species closely related to M. aurascens and M. davidii were selected, Since the mitochondrial genome of M. mystacinus was not available, phylogenetic trees of 64 species were constructed based on single (Table S4, Figure 4) and combined gene sequences of ND1, Cytb, and COX1.
Based on COX1 gene sequences, the phylogenetic tree showed that M. aurascens was closely related to M. nipalensis. Based on the combined sequences of the three genes, the two methods obtained different topological structures: the RAxML tree showed M. aurascens and M. ikonnikovi to be sister species and then clustered with M. alcathoe and M. mystacinus with moderate bootstrap support ( Figure S4). On the contrary, the Bayes tree showed that M. aurascens first clustered with M. alcathoe and M. mystacinus, and then clustered with M. ikonnikovi with strong bootstrap support. However, the Bayes tree obtained high support values, so the result was more convincing (Figure 6). Both topological structures indicated that M. mystacinus and M. davidii were distantly related and should belong to different species. Further, StarBeast was used to construct phylogenetic trees for species classification. Several species closely related to M. aurascens and M. davidii were selected, and the results support our viewpoint ( Figure S5).

Phylogenetic Analysis within Vespertilionidae
Further, we discussed the evolutionary relationships among the genera within Vespertilionidae. As shown in Figures 5 and 6, the evolutionary tree constructed using different gene sets obtained the same topology. Within Vespertilionidae, two main clades were recovered. One clade included the genera Myotis and Murina, and the other clade included the genera Plecotus, Vespertilio, Pipistrellus, and Nyctalus. These results indicated that Myotis is closely related to Murina. Pipistrellus was closely related to Nyctalus, and Hypsugo was closely related to Vespertilio, they clustered together, then with Plecotus, to form the third cluster. All of these structures had strong support values.

Phylogenetic Analysis within Vespertilionidae
Further, we discussed the evolutionary relationships among the genera within Vespertilionidae. As shown in Figures 5 and 6, the evolutionary tree constructed using different gene sets obtained the same topology. Within Vespertilionidae, two main clades were recovered. One clade included the genera Myotis and Murina, and the other clade included the genera Plecotus, Vespertilio, Pipistrellus, and Nyctalus. These results indicated that Myotis is closely related to Murina. Pipistrellus was closely related to Nyctalus, and Hypsugo was closely related to Vespertilio, they clustered together, then with Plecotus, to form the third cluster. All of these structures had strong support values.

Genetic Distance Analysis
To further confirm the genetic relationship between M. aurascens and M. davidii, the pairwise distances of several related species were calculated in the two phylogenetic trees. As is shown in Table 3, the pairwise distances between M. aurascens and M. davidii were 0.138 and 0.161, higher than those between any other species. The pairwise distance be-

Genetic Distance Analysis
To further confirm the genetic relationship between M. aurascens and M. davidii, the pairwise distances of several related species were calculated in the two phylogenetic trees. As is shown in Table 3, the pairwise distances between M. aurascens and M. davidii were 0.138 and 0.161, higher than those between any other species. The pairwise distance between M. aurascens and M. pilosus was the smallest (0.133 and 0.154), indicating that the two species are closely related.
The pairwise distances between M. aurascens and M. davidii calculated based on three genes were 0.134 and 0.156, higher than M. ikonnikovi, M. alcathoe, and M. mystacinus (Table 4). This research showed that M. aurascens and M. davidii were distantly related and clearly need to be classified as different species.

Discussion
At present, there is no report on the study of M. aurascens in China, which may be due to its unclear taxonomic status. Due to their similar appearance, many previous studies have confused it with M. mystacinus. However, through morphological and molecular studies, it was confirmed that M. aurascens and M. mystacinus are different species [7]. Our results also supported this view ( Figure 6).
However, the current classification system classifies M. aurascens as M. davidii, not as a separate species [2]. Based on the phylogenetic tree constructed using a single mitochondrial gene (COX1, Figure 4A), we found that the closest genetic relationship was between M. aurascens and M. nipalensis (Nepal Myotis). In China, M. nipalensis is mainly distributed in Qinghai, Gansu, Xinjiang, Hubei, Jiangsu, and Tibet [31,32]. Although the geographical locations of the two species is far apart, M. nipalensis was previously considered a subspecies of M. mystacinus [33] and was listed as an independent species later [34]. Therefore, it can be inferred that there are certain similarities between the two species in morphology and taxonomic status, which explains why they have the closest genetic relationship in the phylogenetic tree.
According In conclusion, we suggest that M. aurascens should be recognized as a separate species. In addition, we also briefly discussed the inter-genus relationships within the Vespertilionids, and the results showed that Myotis were the closest relatives to Murina.
However, a comparative analysis of molecular and morphological data with the type specimens of M. aurascens and M. davidii was lacking. Therefore, we could not conclude that our research subjects were M. aurascens and M. davidii. However, by combining our results with the previously published molecular data of M. aurascens and M. davidii, we have clear evidence that they are two different species. Next, we need to collect molecular and morphological data on the type specimen to prove our hypothesis.
M. aurascens is mainly distributed in grasslands and forests, and they roost mainly in rock crevices [34]. Using radio-tracking techniques, Korean researchers found that the main range and foraging sites of M. aurascens were concentrated near water bodies, while they mainly inhabit forests and nearby areas during the day [35]. In this study, the M. aurascens specimen was captured in crevices of the coastal cliffs of Hulun Lake (Figure 1), and this habitat type is consistent with previous descriptions. Combined with previous surveys, we found that M. aurascens was distributed near the Shuanmazhuang protection stations and Hulungu in Hulun Lake National Nature Reserve, Inner Mongolia, China.
Since M. aurascens is located in the Hulun Lake National Nature Reserve, its habitat and range are well protected. Although samples were collected several times in the vicinity, the range of M. aurascens is narrow and the number is small compared to the other locally distributed Vespertilio sinensis and Vespertilio murinus. This suggests that further comprehensive and in-depth surveys of the species should be carried out to discover more locations in this region. Habitat degradation due to human activities (e.g., grazing) and climate change are impacting the population of M. aurascens. Hence, it was urgent for us to understand its population status and suggest strategies for conservation of the species.

Conclusions
Bats not only play important ecological and economic roles in pest control and plant pollination but also have important scientific research value. Hence, it is urgent for us to understand their population status and establish a monitoring network to ensure sustainable bat populations in China. In this study, based on phenotypic characteristics, molecular phylogenetic trees, and genetic distances, we found that M. aurascens should be considered a distinct species rather than a synonym of M. davidii. The M. aurascens we discovered mainly lives in the cliff crevices of Hulunbuir Grassland in Inner Mongolia, China. Therefore, this region should strengthen its distribution area survey and expand research to neighboring provinces to find the true distribution of this species.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ani13101629/s1, Figure S1: PCA analysis between M. aurascens and M. davidii based on morphological data; Figure S2: The secondary structure of tRNA gene; Figure S3: The phylogenetic tree constructed from M. davidii data downloaded from NCBI; Figure S4: The ML analyses of phylogenetic relationship based on ND1, Cytb, and COX1 used two methods; Figure S5: The phylogenetic trees for species classification using StarBeast; Table S1: Base composition of the mitogenomes of M. aurascens; Table S2: Frequency and RSCU values of codon in PCGs in the mitogenome of M. aurascens; Table S3: Complete mitochondrial sequence used in molecular phylogenetic analyses in this study; Table S4: COI, ND1 and Cytb sequence used in molecular phylogenetic analyses in this study.
Author Contributions: X.Y., H.D. and H.Z. designed the study; X.Y., X.C., X.G., X.S. and G.S. performed the bioinformatics analysis; X.Y. prepared the draft of the manuscript. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: The animal study protocol was approved by the Ethics Committee of Qufu Normal University (protocol code: 2022002; 7 January 2022). : Tables S3 and S4 provide the NCBI accession numbers of the species used in this study.