Molecular phylogeny and taxonomy of the genus Nectogale (Mammalia: Eulipotyphla: Soricidae)

Abstract The elegant water shrew, Nectogale elegans, is one of the small mammal species most adapted to a semi‐aquatic lifestyle. The taxonomy of the genus Nectogale has received little attention due to difficulties in specimen collection. In this study, we sequenced one mitochondrial and eight nuclear genes to infer the phylogenetic relationship of Nectogale. Phylogenetic analyses revealed two large clades within Nectogale. One clade represented N. elegans, and the other was regarded as N. sikhimensis. The split between N. elegans and N. sikhimensis dated back to the early Pleistocene (2.15 million years ago [Ma]), which might be relevant to the Qinghai‐Tibet Plateau (QTP) uplift. The morphological comparison showed several distinguishing characters within Nectogale: the shape of the mastoids, the first lower unicuspid (a1), and the second upper molar (M2). Overall, the molecular and the morphological evidences supported that the genus Nectogale consists of two valid species: N. elegans and N. sikhimensis.

Nectogale from Sikkim and nominated them N. sikhimensis (Winton, 1899) based on morphological differences with N. elegans. He considered the two specimens browner in color without a clear boundary between dorsal and ventral pelage, and the shorter cusps of the first upper incisor. However, Allen believed that this should be considered a subspecies (Allen, 1938), and this was followed by some scholars (Ellerman & Morrison-Scott, 1951;Smith & Xie, 2009;Wang, 2003). On the other hand, N. sikhimensis has been treated by other scholars as conspecific or synonym of N. elegans (Hoffmann, 1987;Hutterer, 1993Hutterer, , 2005Paradiso, 1975;Repenning, 1967).
The phylogenetic relationship of this genus remains to be studied due to the lack of specimens of N. sikhimensis in previous studies based on morphological (Hoffmann, 1987) and phylogenetic analyses (He et al., 2010;Ohdachi et al., 2006). Is the N. sikhimensis a taxonomically valid species? In this study, we collected specimens of Nectogale species from Sichuan, Yunnan, Qinghai, and Tibet in China, and used one mitochondrial and eight nuclear genes to infer the phylogenetic relationship, estimated the divergence time, and explored the evolutionary history of Nectogale.

| Sampling and DNA sequencing
A total of 44 specimens were collected from Sichuan, Yunnan, and Tibet, including 42 individuals of Nectogale and two of Chimarrogale ( Figure 1 and Table 1). Two individuals from C. styani (Winton, 1899) were used as outgroups for phylogenetic analyses because of the close phylogenetic relationship between the genera Chimarrogale and Nectogale (He et al., 2010;Ohdachi et al., 2006). All specimens were identified based on their morphology and distributions following de Winton and Styan (1899), Smith andXie (2009), andWang (2003).
In addition, some of the other sequences used in molecular analyses were downloaded from GenBank (Tables S1 and S2). The voucher and museum number, location data, and elevation are provided in Table S1.
Total genomic DNA was extracted from the muscle or the liver tissue preserved in 95% ethanol using the phenol/proteinase K/sodium dodecyl sulfate method (Sambrook et al., 1989). All muscle or liver tissues were stored in 100% ethanol at −70°C after DNA extraction for further analysis. We amplified one mitochondrial gene

| Molecular data processes and analyses
All sequences were edited using EditSeq (DNASTAR, Lasergene v7.1) and aligned using MEGA5 (Tamura et al., 2011). Two methods were used to infer phylogenetic relationships: Bayesian Inference (BI) and maximum likelihood (ML). All sequences were divided into two datasets: (1) an eight nuclear gene combined dataset (nDNA); (2) a CYT B gene dataset (mtDNA), including seven sequences of Nectogale and two of C. styani downloaded from GenBank (Table S1). The sequences of C. styani were used as outgroups. Bayesian gene trees were reconstructed independently for the CYT B gene and concatenated nDNA fragments. The Bayesian analysis was performed with BEAST v1.6.1 (Drummond et al., 2012). Each gene was treated as a partition. The analysis used unlinked substitute models, linked clock models, linked trees, an uncorrected lognormal, a relaxed molecular clock model, a birth-death tree prior, and default prior (including a random starting tree). Each analysis ran for 100 million generations and was sampled every 5000 generations. We repeated the analysis 10 times and evaluated the convergence using Tracer 1.5 (Drummond et al., 2012).
In addition, we calculated the average genetic distance between the clades. All calculation results were based on within and between group pairwise analysis using the Kimura 2-Parameter (K2P) model in MEGA5, with 1000 bootstrap replications (Tamura et al., 2011).
We estimated the divergence time using BEAST v1.7.5. Because mitochondrial genes may overestimate the true divergence time (Phillips, 2009;Zheng et al., 2011), the combined nDNA genes dataset was used for the estimation of divergence times. All calibration age constraints were treated as log-normal distributions.
We used three calibration points: (1) The division of Soricinae and Crocidurinae occurred at 36 million years ago (Springer et al., 2018).
The rest of the parameters were set as in the phylogenetic analyses.
We used Network v4.5 to construct phylogenetic trees of each nuclear gene and mitochondrial gene CYT B. Before building the tree, we phased our data in DnaSP v5.10 (Librado & Rozas, 2009;Stephens & Donnelly, 2003), then used the haplotype to construct the median-joining network (Bandelt et al., 1999), and finally used the maximum parsimony strategy to select the optimal phylogenetic tree (Polzin & Daneshmand, 2003).
We explored the population genetic structure in STRUCTURE 2.3.1 (Pritchard et al., 2000) based on nDNA data. The number of genetic clusters (K value) was estimated by the admixture model with correlated allele frequencies. We set the K value from 1 to 10 to find the best value for K. Structure was run with 400,000 after 100,000 runs as burn-in. Structure results were based on 10 independent runs. We got the results from Structure Harvester (https:// taylo r0.biolo gy.ucla.edu/struc tureH arves ter/).

F I G U R E 2
Results of Bayesian phylogenetic analyses of mtDNA dataset (left) and nDNA dataset (right). Numbers above branches refer to Bayesian posterior probabilities (PP). Branch lengths represent substitutions per site.

| Morphological analyses
We obtained 33 skulls from adult individuals. All specimens were deposited at Sichuan Normal University (SCNU), Sichuan Academy of Forestry (SAF), and Sichuan University (MSCU). Three external measurements-namely, head and body length (HB), tail length (TL), and hind foot length (HL) were measured in the field or recorded from original specimen labels. However, they were not used for morphological analyses because these measurements may show considerable inter-observer variation (Jiang et al., 2003). Fourteen craniomandibular variables were measured using a digital caliper graduated to 0.01 mm. These were condyloincisive length (  variables between groups were also tested using an independentsample t test. We used 32 of the 33 relatively complete skulls to analyze morphometric variation using craniomandibular variables by a principal component analysis (PCA). The overall variables were log 10 -transformed before conducting the PCA, and a few missing values were replaced by means. On the basis of the results of our molecular analyses, we assigned two specimens of N. elegans from Tibet (Markam) (see Section 3).

| Phylogenetic relationships and genetic distance
We The K2P genetic distances of CYT B among the 29 samples in the genus Nectogale ranged from 0% to 14.9%, with an average of 5.5%.
The genetic distance within each clade was ≤1.63% (

| Divergence time
The Bayesian analysis of the divergence time tree showed the same phylogenetic relationship of Nectogale as the mitochondrial and nuclear gene tree (Figure 3)

F I G U R E 4
Median-joining network based on mitochondrial and nuclear genes. Each circle represents a single haplotype scaled by its frequency; red dots represent missing or non-sampled haplotypes.

| Mitochondrial and nuclear gene networks and population genetic structure
Sequence characteristics of each gene used in network analyses are presented in Table 3. There were 17 haplotypes in the mitochondrial CYT B gene, and haplotype diversity (Hd) was 0.864. In nuclear genes, the haplotype diversity of VWF (0.919) and GHR (0.855) was higher than that in the other six nuclear genes (Table 3). The network analyses based on the nuclear gene (ADORA3, GHR, and ATP7A) showed that N. sikhimensis and N. elegans did not share haplotypes.
The plot for ∆K and K revealed peaks at K = 2, 4, and 9 ( Figure S2A). The L(K) mean values showed the highest probability of ln for K = 4, after which it stabilized ( Figure S2B). Combining the results of L(K) and ∆K, the best number of clusters in STRUCTURE was 2 and 4 ( Figure S2C). At K = 2, the structure analyses identified two clades in the genus Nectogale ( Figure S2), a result that was consistent with the phylogenetic analysis based on the nuclear gene ( Figure 2). At K = 4, the entire population was divided into four clusters corresponding to the four main sample localities (Zayü of Tibet, Medog of Tibet, Yunnan, and Sichuan).

| Morphological analyses
The skull measurements are given in Table 4. There were significant differences when all individuals in the two clades were subjected to an independent-sample t test for each variable (  (Table 5) Differences also exist in the mandible. In N. sikhimensis, the first half of the horizontal ramus of all specimens from Yunnan and most specimens from Sichuan is stouter than in N. elegans (Figure 6a, b). In N. sikhimensis, the tip of the a 1 is shorter than in N. elegans resulting in a smaller appearance (Figure 6g, h).

| DISCUSS IONS
Many studies have shown a high level of species diversity in the tribe Nectogalini. For example, two new species of the genus Chodsigoa (Kastschenko, 1907) were recently described (C. hoffmanni and C. dabieshanensis) . Some subspecies, such as Chodsigoa furva (Anthony, 1941), Chimarrogale leander (Thomas, 1902), and Neomys milleri (Mottaz, 1907) were recommended species status (Burgin & He, 2018;Chen et al., 2017;Igea et al., 2015;Yuan et al., 2013). New genera may exist in the tribe. For example, Episoriculus fumidus (Thomas, 1913)  Within the genus Nectogale, N. sikhimensis was first described from Sikkim based on differences in fur and teeth from N. elegans, but no molecular studies were performed to date to evaluate the genetic divergence of these two species. In this study, the results of phylogenetic trees based on mitochondrial and nuclear indicated that Nectogale was divided into two clades. Clade A is representing N.
sikhimensis, clade B is N. elegans. Given the large genetic distance and morphological differences between the two taxa, we recover species status of N. sikhimensis and support that the genus Nectogale consists of two species. It is worth noting that the results of the structure analyses revealed that the species of the genus Nectogale were divided into two clusters (K = 2, corresponding to N. sikhimensis and N. elegans) or four clusters (K = 4, corresponding to four main sample localities of specimens). The results are very interesting, especially because this could indicate that cryptic diversity might be present. The results may be important to provide additional insights for future work.
The type locality of N. elegans is in Baoxing, Sichuan, which was also found in Yunnan, Shaanxi, Gansu, and Qinghai, China. On the basis of our collection, it is also distributed in Markam, eastern Tibet of China (Figures 1 and 2). The type locality of N. sikhimensis is Sikkim (India). This species is mainly distributed in southeast QTP of China and was also found in Bhutan, Myanmar, and Nepal (Smith & Xie, 2009). It seems that Mt. Gaoligong and Mt. Boshulaling are the boundaries between the two species.
The role of geographic isolation in speciation and diversification has been demonstrated in previous studies (Qu et al., 2014;Xing & Ree, 2017). Although rivers are not barriers for water shrews (Yuan et al., 2013) Springer et al., 2018), instead of 20 million years ago (He et al., 2010;Reumer, 1994).
This discordance between mitochondrial and nuclear gene trees has been observed in our study. For example, the two subclades B1 and B2 were strongly supported in mtDNA gene trees but not in nDNA trees. This discordance often was attributed to several factors, i.e., mitochondrial capture (Dong et al., 2014), explosive speciation (Krause et al., 2008), incomplete lineage sorting (Edwards, 2009), and introgression (Funk & Omland, 2003;Yannic et al., 2010). The most likely reason for cyto-nuclear incongruence is that the mitochondrial gene CYT B was captured during the divergence of the subclades B1 and B2. This process has been observed in many cases (e.g., Bryson et al., 2010;Dong et al., 2014;Markova et al., 2013;Tang et al., 2012).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
New DNA sequences in this study were deposited in GenBank (Accession numbers ON160936-ON161124, ON219777-ON219792). (https://www.ncbi.nlm.nih.gov/).