Abstract
The Himalayan Shrew belongs to the genus Soriculus , which is a
monotypic genus. We investigated the phylogeographic patterns,
speciation, divergence time, and demographic history ofSoriculus nigrescens in
southwestern China. A total of 128
samples from 29 regions were analyzed for both mitochondrial and nuclear
genes. Phylogenetic analyses of the
concatenated mtDNA and nuDNA data revealed three highly divergent
lineages within the S.
nigrescens . One group represented
subspecies S. n. minors , another group represented subspeciesS. n. nigrescens , which is
made up of two Clades A and B. The
species delimitation analyses, based on two methods, supported the
species status of the two Clades ofS. n. nigrescens . In addtion,
it was found that individuals at different altitudes in Motuo were
divided into two Clades. Bayesian skyline plotting analyses and
ecological niche modeling also supported demographic and range
expansions during the LGM for S. n. nigrescens .
We propose that S. n.
nigrescens appears to be composed of two putative species, and S.
n. minor should be elevated to species status. Our study also
suggested that climate change since the Miocene periods and the uplift
of the QTP may have resulted in the diversification and speciation ofS. nigrescens .
Keywords: The Himalayan Shrew;
Population demography; Species delimitation; Population expansion
Introduction
It seems that the speciation driving
patterns, genetic variation and population structure of many species are
all affected by paleoclimate fluctuations and paleogeological events
(Avise, 2009; He and Jiang, 2014; Feng et al., 2016; Ye et al.,
2016a, 2016b). Both current and
historical environmental factors, such as topography, climate and
habitat, were strongly affected evolution of small mammals (Kent et al.,
2014). The shrews have limited dispersal abilities due to their small
body, high metabolism, and low vagility (Churchfield, 1990; Churchfield
et al., 2012). As a consequence, the evolution of terrestrial shrews
were strongly affected by Climatic and Topographic situations (Vega et
al., 2010; Wierzbicki et al., 2011; Willows-Munro and Matthee, 2011).
Shrews could provide excellent models for exploring the genetic
consequences of climatic changes and topographic isolation (Chen et al.,
2015).
The
distribution pattern of biodiversity and its causes are one of the core
contents of biogeography research, which is of great significance for
understanding the evolutionary laws of biodiversity, the process of
species formation, and formulating the protection and management of
species diversity (Chen et al.,
2015). There are a total of 25 ecological hotspots in the world, all of
which have high biodiversity (Myers
et al., 2000). Southwest China is a
hotspot of biodiversity in East Asia due to the high level of endemism
and extensive biodiversity (He and Jiang, 2014).
Geological and climatic events
affect the morphology, spatial structure and genetic structure of
animals, especially for thermostatic
vertebrates
(Wu et al., 2013). The uplift of the
Qinghai-Tibet Plateau (QTP) and the Quaternary glacial period, orogeny,
the formation of the Asian monsoon and the violent Quaternary climate
oscillations have profoundly affected the distribution and genetic
structure of existing organisms (An
et al., 2001; He and Jiang, 2014; Lei et al., 2014; Qu et al., 2015; Ye
et al., 2016a,b). During the ice
age, although the ice sheet was not completely covered, it had a huge
impact on species.
Research by
Reumer
(1987) indicated that temperature and humidity are the key factors
affecting the distribution, number and species of shrews.
The Himalayan Shrew
(Soriculus nigrescens ) Gray,
1842 belongs to the Soricinae, which
is distributed in southwest China,
Nepal, India and Bhutan (Smith and Xie, 2010). And it is a monotypic
genus. The classification of Soriculus nigrescens was confusing
before. Based on morphological data, Motokawa (2003) dividedSoriculus nigrescens into two subspecies, S. n. nigrescens(Gray, 1842) and S. n. minor Dobson, 1890 (Motokawa, 2003).Soriculus nigrescens minor was considered as senior synonym ofS. radulus Thomas, 1922, and S. n .nigrescens was regarded
as senior synonyms of S. aterrimus , S. sikimensis ,S. n. pahari , S. n. caurinus , and S. n. centralis(Motokawa, 2003). But currently there is a lack of molecular
phylogenetic research data on Soriculus nigrescens .
In this study, we analyzed the genetic structure of the population,
calculated the time of differentiation and the historical dynamics of
the population, reshaped the niche model of Soriculus nigrescens ,
and studied the relationship between its lineages and system evolution,
in order to explore the genetic diversity and the mechanism of species
diversity formation in the Himalayas in China and provide data for
differentiation mechanism of small mammals in the Himalayas.
Materials and methods
Sampling, sequencing and alignment
From 2009 to 2019, a total of 124
samples of theS. n. nigrescens were
collected from 28 regions in Tibet of China. Three samples of theS. nigrescens from Yunnan were downloaded from the Genbank.
Additional sequence data of Soricidae
were downloaded from the GenBank
and included in the analysis (Table S1). Detailed specimen and locality
data are presented in Fig 1A and
Table S2.
Total genomic DNA was extracted from muscles or livers using a
chloroform procedure (Sambrook and Russell, 2006). The
mitochondrial
Cyt-B gene, nuclear genes
APOB,
BRCA-1 and RAG-2 were amplified by highly conserved primers, detailed
primer data are presented in
Appendix Table S3. The same PCR
primers and primer walking were used for sequencing, detailed PCR
amplification and sequencing methods followed Chen et al., (2015). The
Cyt-B, APOB, BRCA-1 and RAG-2
sequences were aligned and examined by eye in MEGA v5 (Tamura, et al.,
2011).
Genetic diversity and structure
In this study, DnaSP 5.10 was used
to calculate the haplotype number (H), haplotype polymorphism, haplotype
diversity (HD), and nucleotide diversity (π) of the sample, to analyze
the polymorphic data of DNA (Librado and Rozas, 2009). Estimated genetic
divergence among the major mtDNA lineages (see below) through the paired
Kimura two-parameter (K2P) model implemented in MAGA v5.0 (Tamura et
al., 2011). Each nuclear gene calls the PHASE typing algorithm (Stephens
and Donnelly, 2003) in DnaSP v5.10 (Librado and Rozas, 2009), and uses
the MP method to remove redundant connections. Haplotype median-joining
networks (MJN) of genealogical relationships for nuDNA were also
constructed with the software Network v4.6 (Bandelt et al., 1999).
Phylogenetic analyses
We used BEAST v1.6.1 for Bayesian phylogenetic reconstructions (Drummond
and Rambaut, 2007). Bayesian gene trees were reconstructed independently
for concatenated Cyt-B and each nuDNA gene. The jModeltest v2
(Darriba et al., 2012) was used to
determine the model of evolution for mtDNA and nuDNA genes by Bayesian
Information Criterion (BIC)(Table S4). The prior value of the parameter
uses the default value, the starting tree uses a random spanning tree,
and each data set is run at least 4 times. Each operation is 100 million
generations, and every 5,000
generations are sampled. We used Tracer v1.5 to analyze the log file in
the calculation result to check the convergence of the calculation. When
the posterior value of the data set is greater than 200 (ESS≥200), the
analysis is considered stable and the result is credible (Drummond and
Rambaut, 2007). After that, according to the posterior distribution,
Tree Annotator is used to determine the burn-in. In addition, when using
Figtree v1.4.3 to observe Bayesian trees, only when the DNA sequences of
certain groups are clustered on nodes with a posterior probability value
greater than 0.95 (PP> 0.95), it is proved that the group
is highly supported (Huelsenbeck and Rannala, 2004).
Species delimitation
First, we
constructed the species tree by
using the *BEAST model (Heled and Drummond, 2008). The species tree was
reconstructed for S. nigrescens based on a coalescent-based
method implemented in *BEAST. Chodsigoa hypsibius was selected as
outgroup. Only individuals with both mtDNA and nuDNA data were included
in the species delimitation analysis. In order to meet the construction
requirements of the *BEAST model, the model we adopted is the same as
the phylogenetic analysis. *BEAST program setting: The main parameter
setting is consistent with the analysis of the divergence time, the
running generation is modified to 100 million generations, and every
5,000 generations are sampled.
Second, we used Bayesian species-delimitation method to delineate
species boundaries (Yang and Rannala, 2014). We tested the validity of
our assignment of 2 putative species (Clade A and Clade B) in BPP v3.1.
Only dataset 1 (nuDNA) and dataset 2 (mtDNA + nuDNA) were included in
the analysis. Two alternative rjMCMC algorithms (algorithms 0 and 1)
implemented in BPP. We used algorithm 0 with fine tuning parameters for
subsequent analysis as the original analysis showed that algorithms 0
and 1 produce similar results. We used Gamma-distributed priors (G) to
specify the ancestral population size (θ) and root age (τ). The three
combinations were modeled to allow a range of speciation histories: G
(2, 2000) for θ and G (2, 2000) for τ ,
G (1, 10) for θ and G (2, 2000) for
τ and G (1, 10) for θ and G (1, 10) for τ. Each rjMCMC was run for
100,000 generations and sampled each 100 generations after discarding
10,000 generations as pre-burn-in.
Divergence time estimates
Cyt-B,
APOB and BRCA-1 were used to calculate the divergence time ofS. nigrescens population
through the BEAST v1.6.1 (Drummond,
et al., 2012). Following He et al., (2018), we set the mean to 36 Ma,
with a standard deviation of 0.135, so that the median of the prior was
35.7 Ma and the 95% CI was 28.6 - 44.5 Ma. BEAST analyses used unlinked
substitute model, linked clock models, linked tree, a random starting
tree, a birth-death process tree prior, a relaxed lognormal clock model,
and the programs default prior distributions of model parameters. Each
analysis was run for 100 million generations and sampled every 5, 000
generations. Convergence was
accessed using Tracer v1.4 to confirm the effective sample sizes (ESS)
of all parameters were higher than 200 (Rambaut and Drummond, 2007). The
first 10% of the trees were discarded as the burn-in.
Population demography
We used Bayesian skyline plots (BSP) implemented in BEAST v1.6.1
(Drummond et al., 2012) to examine the historical demographics ofS. n. nigrescens with the
Cyt-B and nuDNA data. The HKY + G and
GTR + I + G nucleotide substitution models selected by jModeltest 0.1
were used. The BSP analysis was
executed iterations with 100 million and sampled at every 5000th tree,
with the first 10% discarded as the burn-in.
Convergence
was accessed using Tracer v1.4 to confirm the effective sample sizes
(ESS) of all parameters were higher than 200 (Rambaut and Drummond,
2007).
Ecological niche analysis
We used ecological niche modeling (ENM) to predict the distribution ofS. n. nigrescens (Clades A, B in Fig. 1B) at the present time and
during the last glacial maximum (LGM).The maximum entropy algorithm
implemented in MAXENT v3.3.3 (Phillips et al., 2006; Phillips and Dudík,
2008) was used to model the ecological niches. A total of 28
distribution localities were collected. We used three environmental
layers in different periods: Last Interglacial (LIG), Last Glacial
Maximum (LGM) and Present. The environment layer comes from WorldClim
(Hijmans et al., 2005) (website
http://www.Worldclim.Org). The main climatic elements considered and
reflected are monthly precipitation and average precipitation, minimum
and maximum temperature, seasonal variation and limit level. The Maxent
analysis was used with 80% of the species records for training and 20%
for testing the model, 10 subsample replicates and a maximum of 2000
iterations. The area under the curve (AUC) of the receiver operating
characteristic (ROC) plot was used for model evaluation.
RUSULT
Sequence characteristics and genetic diversity
The 127 individuals used in this study were from 28 sampling sites, and
3,093 bp of gene sequence was obtained after sequencing, including 1,953
bp of nuclear gene sequence (APOB [513 bp], BRCA-1 [765 bp],
RAG-2 [675 bp]) and 1,140 bp mitochondrial sequence
(Cyt-B [1,140 bp]). All coding
genes fragments were tested in MAGA v5
(Tamura et al., 2011), and the
results showed that there was no premature termination of coding genes.
The mitochondrial gene Cyt-B has the most mutation sites compared with
nuclear genes, and mitochondrial genes have higher genetic
polymorphisms. For all the Soriculus nigrescens samples, the
haplotype diversity (h) and nucleotide diversity (π) of mitochondria are
relatively high (h = 0.97141, π = 0.07291). Among all amplified genes,
mitochondrial Cyt-B had the highest
haplotype, and 51 haplotypes were
identified in Cyt-B. Among nuclear genes, RAG-2 has the highest
haplotype number, and 25 haplotypes were identified in RAG-2. The
nucleotide diversity and haplotype diversity of mitochondrial genes of
Clade B are significantly higher than those of Clade A
(Table S5). The average distances
between species were 7.0% (Clade A and Clade B), 12.4% (Clade A and
Clade YN), and 11.3% (Clade B and Clade YN)
(Table 1).
Phylogenetic relationship and population genetic structure
The phylogenetic tree results of the mitochondrial data set, nuclear
gene data set, and the combined data set showed that S. n.
nigrencens can be divided into three clearly Clades (Clade A, Clade B
and Clade YN), with higher posterior probability (Fig 2). But the
phylogenetic status of Clade YN is not completely consistent in the
three data sets. But Clade A and Clade B are always the closest
relatives. Clade A is from Yadong County, Dingri County, Nielamu County,
and Motuo County in Tibet. Clade B are from Bomi County, Milin County,
Gongbujiangda County, Bayi District and Motuo County in Tibet. In Motuo
County, Clade A and Clade B exist at the same time, Clade A is in the
low-altitude area, and Clade B is in the high-altitude area. The MJN for
the nuDNA data showed a similar structure to the phylogenetic tree,
Clade A, Clade B and Clade YN have no crossover phenomenon and can be
completely divided into three Clades.
We used the nuclear genes APOB, BRCA-1, and RAG-2 of Clade A and Clade B
to analyze the mediation network diagram to obtain 3 MP trees. The
results showed that the MP tree results of all nuclear genes of S.
n. nigrescens are consistent with the previous mitochondrial
phylogenetic tree results, and Clade A and Clade B are clearly divided
into two clades (Fig S1). In all nuclear gene mediation network
diagrams, Clade A and Clade B have no crossover phenomenon, and can be
completely divided into two branches, which is consistent with the
results of the previous nuclear gene phylogeny diagram.
Species delimitation
The *BEAST species-tree and BPP
analyses gave identical results. The results of the *BEAST species-tree
delimitation analysis of individuals with both the nuDNA (APOB, BRCA-1
and RAG-2) and Cyt-B data are shown in Fig 3A and Table S6. The results
of the *BEAST species tree supported
the S. nigrescens as three
species, andS. n. nigrescens divided into
two species. When using either dataset 1 (nuDNA) and dataset 2 (mtDNA +
nuDNA) of S. n. nigrescens , BPP consistently supported the
validity of 2 hypothetical species with high posterior probabilities (PP
≥ 0.96)(Table S6).
Divergence time estimation
Divergence dating estimated that S. n. nigrescens diverged withS.n.minor 11.11 MYA during the early
late Miocene.
The divergence time of S. n.
nigrescens further differentiated into two clades of Clade A and Clade
B was at 5.69 MYA during the late Miocene (PP ≥ 0.95). The results
showed that it was affected by QTP change. The first uplift of QTP
promoted the differentiation of S. n. nigrescens and S. n.
minor . And the second uplift of QTP promoted the differentiation of
Clade A and Clade B. It is consistent with the uplift of QTP in the Late
Miocene. Since then, Clade A has differentiated several times in the
Pliocene epoch (PP ≥ 0.95); Clade A has an earlier differentiation time
than Clade B, which is also in the Pliocene epoch (PP ≥ 0.95)(Fig 3B).
Historical demography
In the combined analysis of Clade A and Clade B, the results of the
BSP combined with nuclear gene and
mitochondrial gene analysis show that the curve gradually rises from the
present to the LGM (26.5 - 19.0 ka), indicating that the population
expansion of S. n. nigrescens has occurred (Fig 4A). In the
analysis of mismatch distribution, the mitochondrial gene had multiple
jagged peaks, and all of nuclear gene mismatch distribution showed
unimodal distributions (Fig 4B).
Ecological niche modeling
In this study, the ENM of S. n. nigrescens in three different
climatic environments (current, LGM and LIG) have achieved good
prediction results (Fig 1B). The results of divergence time show that
the history of differentiation between the two branches of Clade A and
Clade B can be traced back to the
Late Miocene. The results of ENM
show that the population size of each branch of the S. n.
nigrescens migrated and changed during the severe fluctuations of the
last glacial maximum and the last interglacial period. Since the LIG,
the distribution range of S.
n. nigrescens has gradually expanded. Compared with the BSP results,
the ENM results show that the population expansion time may be earlier.
Discussion
Biogeographic inferences
First, our study supports allopatric occurrence of species within theS. n. nigrescens . In this study, these species have obvious
systematic geographic characteristics, and their distribution ranges do
not overlap. Clade A is mainly distributed in the altitude areas of
Yadong County, Nielamu County, Dingri County and
Motuo County in Tibet. In the Motuo
area, the collected specimens are deciduous broad-leaved forest areas
below 1200m altitude. And there are also specimens of the S. n.
nigrescens in the residential areas of Motuo
County. The specimens of Clade B
were collected from Milin, Bomi,
Gongbujiangda, Bayi and Motuo. And
these samples from the Motuo area came from the high-altitude area of
Motuo, which is bordering the Milin. Clade YN was collected from Pianma,
Yunnan Province (He and Jiang, 2014). The divergence of these lineages
may be to the uplift of the QTP and the accompanying climate
fluctuations. The QTP experienced four phases (12 - 14 Ma
(Turner et al., 1993; Coleman and
Hodges, 1995; Dettman et al., 2003; Lu et al., 2004; Sun et al., 2005),
6 - 8 Ma (Molnar et al., 1993; Molnar, 2005; Song et al., 2007), 3.6 - 4
Ma (Li and Fang, 1999; Song et al., 2007; Nie et al., 2008) and 1.2 -
0.6 Ma (Li and Fang, 1999). Each phase of a geological event can lead to
habitat fragmentation, which may prevent spread and promote vicariant
speciation in many taxa (He et al.,
2010; Lei et al., 2014).
The first uplift of QTP promoted
the differentiation of S. n. nigrescens and S. n. minor .
And the second uplift of QTP promoted the differentiation of Clade A and
Clade B. The rapid rise of QTP and the evolution of the Asian monsoon
have reshaped the landscape of East Asia.
Second, we found that climate change
in the Miocene affected the genetic pattern and population history ofS. n. nigrescens .
Climate
change and cyclical Quaternary glaciations have been revealed as the
major factors affecting the demographic history of many extant species
in Europe and North America (Hewitt, 1996, 2004; Ursenbacher et al.,
2008; Fijarczyk et al., 2011). BSP demographic analyses indicated
population expansion for S. n. nigrescens . The atypical star-like
radiation from the MJN analysis implied the presence of a bottleneck for
the two clades in the past. Our ENMs analyses revealed the predicted
range has changed between the LGM and the present for S. n.
nigrescens (Fig. 1B). Furthermore, the BSP analyses showed that the
expansion events began from approximately 0.02 MYA, suggesting the role
of Pleistocene climatic changes in driving the demographic expansion of
the S. n. nigrescens .
Nuclear and mtDNA genetic structure
The analyses based on the Cyt-B and nuDNA data revealed a pattern which
is consistent with the distribution for the S. n. nigrescens . The
results of Cyt-B and nuDNA both support the existence of S. n.
nigrescens as two clades.
The inconsistency between mitochondrial and nuclear data may be caused
by many factors. The incomplete lineages ordering of ancestral
polymorphism and the different evolution rate of different markers may
cause inconsistencies between mitochondrial and nuclear gene data (Toews
and Brelsford, 2012; Lin et al., 2014; Ye et al., 2016b). Incomplete
lineages division is often used to explain the inconsistency between
mitochondrial and nuclear results. If incomplete lineages sorting
results in the observed pattern of mitochondrial DNA, it can be expected
that the divergence pattern of nuclear DNA will be shallow, as a
relatively large effective population size will result in a relatively
long sorting time (Ye et al., 2016b). We found less nuclear genetic
divergence compared to the mtDNA, and this is consistent with incomplete
lineage sorting. The nuDNA appeared to be conservative and have lower
evolutionary rate than mtDNA, which underestimated the genetic structure
of S. n. nigrescens . Moreover, the rate of evolution may also
have a strong impact on the ability to detect differentiation. The slow
rate of evolution may not produce enough polymorphisms, and the existing
genetic structure cannot be well recognized, and the excessively high
rate of evolution may even mask the signals of high differentiation
(Hedrick, 1999; Brito, 2007). The nuDNA appeared to be conservative and
have lower evolutionary rate than mtDNA, which underestimated the
genetic structure of S. n. nigrescens . Thus, the different
evolution rates of mtDNA and nuDNA markers may lead to the inconsistency
of nuclear mitochondria revealed in this study.
Taxonomic consideration
Traditional taxonomy is based on the morphological characteristics of
species for identification and description, which may underestimate the
species diversity due to the phenotypic plasticity and local
morphological adaptations to the similar environmental habitats (Wan et
al., 2013). There are cryptic species of mammals that would likely not
be recognized based solely on classical studies of morphology of voucher
specimens housed in museums (Kaliontzopoulou et al., 2010; Barley et
al., 2013). Motokawa (2003) divided S. nigrescens into two
subspecies, S. n. minor andS.
n. nigrescens based on morphological judgment. In this study, we
further divide S. n.nigrescens into two clades, Clade A and Clade B. Each of them
formed a well-supported clade on both the Cyt-B and nuclear gene trees.
And species delimitation based on the method BPP confirmed their species
status with high supports. In addition, by extensive sampling, our mtDNA
genealogy and MJN results all supported the existence of two highly
supported genetic clades (Clade A, B) for S. n. nigrescens .
The pairwise distance among the
three clades (0.070, 0.113, 0.124) is obviously higher than the
divergence threshold for identifying mammal cryptic species.
Since their non-overlapping
geographic distributions, it suggested the extensive lineage
diversification in the S. n. nigrescens . The BPP species
delimitation analysis indicated the two clades were candidate species.
Furthermore, Clade YN always clustered separately in molecular analysis,
and it does not overlap with other two clades. Therefore, the species
delimitation results in this study reveal that there are two species in
the S. n. nigrescens , and S. n. minor should be elevated
to species status.
The type specimen of S. n. nigrescens was collected from
India, West Bengal, Darjeeling, and
kept in the Natural History Museum (London). Among the collection areas
in this study, Yadong County is the area where is closest to the type
specimen production area of S. n. nigrescens . The samples
collected in Yadong County were identified as Clade A. In addition, the
samples collected from Nielamu County were also identified as Clade A.
In this study, we found that samples collected from central Nepal were
clustered with Clade A. Therefore, we hypothesized that the specimens of
Clade A and the type specimens were of the same species.
Conclusions
In the present study, we obtained molecular sequences ofS. nigrescens in Tibet and
Yunnan province. We analyzed the mtDNA and nuDNA data to study the
diversification, phylogeographic structure, and evolutionary history ofS. nigrescens . We detected
three highly supported Clades within the S. nigrescens based on
mtDNA and nuDNA genealogy analysis.
We
propose that S. n. nigrescensappears to be composed of two putative species (Clade A and Clade B),
and S. n. minor should be elevated to species status.
The urgent morphological description
needs to further verify our hypothesis. Our study also suggested that
climate change since the Miocene periods and the uplift of the QTP may
have resulted in the diversification and speciation of S.
nigrescens .