Molecular relationships of the Israeli shrews (Eulipotyphla: Soricidae) based on cytochrome b sequences

2.1 Sample origin

Thirty-one samples were used in this work. We received eight C. leucodon, 12 C. suaveolens, six S. etruscus, and three C. ramona samples (all from Israel) from the Steinhardt National Collection of Natural History, Zoological Museum at Tel Aviv University (Israel). The location of the samples is indicated in Figure 1A. In addition, the Field Museum of Natural History (Chicago, USA) kindly provided two samples of Crocidura nana from Tanzania, which were also considered in order to examine the possible relationship of this species to C. ramona. Both ramona and nana share a small size and a greyish tinge to their fur (Dobson 1890; Ivanitskaya et al. 1996). Voucher number, collection location, and date of all samples are provided in Supplementary Table S1.

Figure 1:

Shrew sample locations. A. samples sequenced in this work. B. shrew specimen present in the Steinhardt National Collection of Natural History, Zoological Museum at Tel Aviv University (Israel). The two major biomes of Israel: the Mediterranean biome, and the desert, are indicated in green and yellow, respectively. Red squares, blues stars, green circles and black triangles indicate C. ramona, S. etruscus, C. suaveolens gueldenstaedtii and C. leucodon specimen respectively. The museum records encompass 79 C. leucodon, 594 C. suaveolens gueldenstaedtii, 13 C. ramona and 443 S. etruscus.

2.2 DNA isolation

Two types of museum samples were used. The first was that of tissue samples from specimens captured relatively recently (2002–2015) and preserved in 70–100% ethanol. The second type of sample was that of skin pieces of older specimens (collected 1969–1995) preserved as study skins (Supplementary Table S1). Each type of sample was extracted following a different protocol.

For the specimens preserved in alcohol, a sample of the leg muscle (about 0.5 gr), ear (about 0.5 gr/sample), or kidney/heart (about the same weight) was cut into small parts and digested using 0.5 mL of lysis buffer (1% SDS, 10 mMTris-HCl pH8, 125mMNaCl, 5 mM EDTA, and 0.5 mg/mL Proteinase K) at 55 °C. Following homogenization, the DNA was extracted using a standard phenol-chloroform protocol followed by ethanol/sodium-acetate precipitation (Sanbrook et al. 1989).

For the four specimens preserved as skins, we used the QiagenDNeasy Blood and Tissue Kit following the protocol of Iudica et al. (2001) for small bones from dried mammalian museum specimens. This protocol starts with washing the tissue 3 times in 250 µL of 1% PBS for 10 min at 55 °C. This step allows the removal of any inhibitors and residual fixatives present in the skin specimens. Iudica et al.’s (2001) protocol includes a long digestion time (in our case digestion was conducted for 15–36 h, depending on the sample) with repeated additions of proteinase K (specifically 10 μL every 6 h, overnight 20 μL). The long digestion was necessary since the skin tissues used were hard and difficult to digest.

2.3 Amplification and sequencing of the cytochrome b

Amplification of the complete cyt b gene was conducted differently for the ethanol-preserved samples and the skin samples. For the ethanol-preserved samples the DNA quality following extraction allowed us to amplify large fragments. Consequently, amplification of the cyt b gene was conducted with the primers D1 and R1, which are located in the tRNA-glu and tRNA-thr, respectively (primer sequences are presented in Supplementary Table S2). The PCR amplifications were conducted using the BIOTAQ™ DNA polymerase (Bioline, London). Following a denaturation at 94 °C for 2 min, the PCR was set for 40 cycles of denaturation at 94 °C for 40 s, annealing at 59 °C for 40 s, and elongation at 72 °C for 2.5 min. The 40 cycles were followed by a final elongation step at 72 °C for 10 min.

Because the yield of the amplification was insufficient to directly sequence the PCR fragments obtained, these products were re-amplified. The re-amplifications were conducted using the nested primers D2 and R2, which are located in the tRNA(Glu) and tRNA(Pro) genes, respectively. The PCR conditions were as described above.

For the four samples preserved as skins, the DNA was too degraded to enable the amplification of fragments longer than 400 bp. The sequencing of the complete cyt b gene was conducted by amplifying small overlapping fragments of ∼200–350 bp. Different primer sets were used for each species (all primer sequences are given in Supplementary Table S2). Following denaturation at 94 °C for 3 min, the PCR was set for 34 cycles of denaturation at 94 °C for 45 s, annealing at 51 °C for 45 s, and elongation at 72 °C for 45 s. The 34 cycles were followed by a final elongation step at 72 °C for 10 min.

PCR products were directly sequenced using Big Dye Terminator v1.1 (Applied Biosystems) on an ABI 310 sequencer by the DNA Sequencing Unit at the G.S. Wise Faculty of Life Sciences, Tel Aviv University.

2.4 Phylogenetic reconstructions

Four datasets were used to reconstruct the phylogenetic position of the obtained Israeli sequences. The first dataset comprised representatives of the Crocidura diversity. Specifically, sequence searches were conducted on the National Center for Biotechnology Information (NCBI) nucleotide collection database, with the queries ‘Crocidura cyt b’ and ‘Crocidura cytochrome b’, and all sequences longer than 500 bp were downloaded. We also downloaded all sequences from Dubey et al.’s (2008) dataset, and sequences from the same individuals were concatenated. From the list of sequences obtained from the NCBI searches and the Dubey dataset, we selected the longest sequence for each Crocidura species. In addition, two C. leucodon, two C. suaveolens, three C. ramona, two S. etruscus (all from Israel), and two C. nana (from Tanzania), obtained as described above, were added to the dataset. We also added S. etruscus (LC126597 and DQ630397), S. dayi (DQ630389 and DQ630432), two S. montanus (GQ290374 and DQ630388), and three S. murinus (LC126460, LC126565, and LC126577) sequences from different countries as outgroups (Dubey et al. 2008) (Supplementary Table S3).

The C. leucodon dataset was created by downloading C. leucodon cyt b sequences from NCBI using the keywords “C. leucodon cyt b” and selecting only nucleotide results. A total of 56 sequences longer than 500 bp were downloaded. Eight newly-obtained C. leucodon sequences from Israel were added to this dataset. The tree was rooted with cyt b sequences of Crocidura musseri (FJ813927 and FJ813929) and Crocidura obscurior (KC684154, KC684155, and KC684158) downloaded from NCBI. Again, the outgroup choice was based on Dubey et al. (2008) (Supplementary Table S4).

The C. suaveolens complex dataset was created by downloading C. suaveolens cyt b sequences from NCBI using the keywords “C. suaveolens cyt b” and selecting only nucleotide results. A total of 203 sequences longer than 500 bp were downloaded. Four Crocidura shantungensis cyt b sequences (KF144163, AB077278, KF144159, and AY843447) and one Crocidura zarudnyi sequence (AY925211) were also added since they are known to be closely related to C. suaveolens (Dubey et al. 2008). A total of 12 C. suaveolens sequences from Israel were added to this dataset. The tree was rooted with cyt b sequences from Crocidura brunnea (DQ059025), Crocidura lasiura (AY843503), and Crocidura nigripes (DQ059024), since these species are closely related to C. suaveolens (Dubey et al. 2008) (Supplementary Table S5).

The S. etruscus dataset was prepared by first downloading Suncus cyt b sequences from NCBI using the search “Suncus cyt b” and selecting only nucleotide results. A total of 250 sequences was downloaded and sequences shorter than 500 bp were removed. From this dataset we selected all representatives of S. etruscus, Suncus fellowesgordoni, S. madagascariensis, S. malaynus, and S. dayi. In addition, we selected three S. montanus, five S. murinus, and two S. stoliczkanus sequences as outgroups. This choice was based on previous phylogenetic works (Dubey et al. 2008; Meegaskumbura et al. 2012a, b; Omar et al. 2011). Finally, we added six S. etruscus sequences from Israel to this dataset (Supplementary Table S6).

The nucleotide sequences were aligned using a codon-based alignment as implemented in Geneious 7.1.9, with MAFFT version 7.0.1.7 under the L-ins-i algorithm (Katoh and Standley 2013). Following alignment, columns with at least 50% gaps were excluded. The alignment files are provided in Nexus format in Supplementary Datasets S1–S4. Phylogenetic trees were reconstructed for each dataset separately under the maximum likelihood (ML) and Bayesian criteria. ML analyses were performed with the program RaxML version 8.1.2 (Stamatakis 2014) as implemented in RAxMLGUI 1.5b2 beta (Silvestro and Michalak 2012). The analyses were run with the options “ML + thorough bootstrap”, “20 runs”, “1000 bootstrap repetitions”, and “GTRGAMMA”. Bayesian reconstructions were performed with MrBayes version 3.2.6 (Ronquist et al. 2012) under the GTR + Gamma model. The parameters of the analyses were: two runs with four chains each, sampling every 100 generations, and Burninfrac set to 0.25. The analyses were run for 15,000,000 generations for the C. leucodon and S. etruscus datasets and for 30,000,000 generations for the C. suaveolens complex and the entire Crocidura dataset. For each dataset, we verified that the average standard deviation of split frequencies was below 0.01 before the burnin threshold was reached. We also verified that the Potential Scale Reduction Factor (PSRF) parameters were always close to 1.0 at the end of the run.

2.5 Network analyses

Median-joining networks (Bandelt et al. 1999) were inferred using cyt b haplotypes and the program NETWORK v 10.0.0.0 under default parameters (available at http://www.fluxes-engineering com/sharenet.htm). The datasets used in the network analyses were constructed from the matrices used in the phylogenetic analyses by removing outgroups and shorter sequences as well as excluding all positions containing ambiguous or missing data. The C. leucodon dataset comprised 51 taxa and 1041 positions, of which 155 were variable. The C. suaveolens complex dataset comprised 75 taxa and 901 positions, of which 76 were variable. The S. etruscus dataset comprised 23 taxa and 516 positions, of which 98 were variable.

3.1 Phylogenetic relationships among Crocidura species

The phylogenetic tree reconstructed based on the Crocidura dataset, and which encompassed representatives of the Crocidura species diversity, agreed overall with the findings of Dubey et al. (2008). Because our analyses were based on a single gene (1140 bp), whereas Dubey et al. (2008) had used four genes (3306 bp), branch support is much lower in our case. It should be noted that the topological differences between Dubey’s trees and the trees in Figure 2 and Supplementary Figure S1 only relate to the lowly supported nodes. Specifically, our results recover the monophyly of the Afrotropical clade and the Asian clade described by Dubey et al., but with no support (ML Bootstrap Percentage, BP < 50%; Bayesian Posterior Probabilities, PP = 0.65 for the Asian clade and PP = 0.98 for the Afrotropical clade). The monophyly of the Old-World clade could, however, not be recovered.

Figure 2:

Maximum likelihood tree of Crocidura cytb sequences with emphasis on the Asian and Old World species. Phylogenetic relationships inferred from a matrix of 1,128 nucleotide positions for 131 individuals. Maximum likelihood bootstrap supports above 50% and Bayesian posterior probabilities above 0.70 are indicated near the corresponding node separated with a slash. Sequences obtained in this work are indicated in bold. Authors of sequence data in Supplementary Table S3.

Regarding the Israeli samples, the phylogenetic results confirmed the presence of four distinct species: S. etruscus, C. leucodon, C. suaveolens gueldenstaedtii, and C. ramona (Figure 2 and Supplementary Figure S1). The first three species strongly clustered (BP = 100, PP = 1) with sequences of the same species. In agreement with Dubey et al. (2008), C. ramona sequences formed a distinct lineage whose phylogenetic position could not be determined since the monophyly of the Old-World clade was not supported in our analyses (Figure 2). The C. ramona sequences are clearly unrelated to the C. nana sequences that we obtained from Tanzania, which are nested within an Afrotropical clade with Crocidura jouvenetae, Crocidura crossei, and Crocidura lusitania (BP = 86, PP = 1.0; Supplementary Figure S1).

3.2 Phylogenetic relationships within the C. leucodon clade

The C. leucodon tree recovered the monophyly of the main known C. leucodon clades (Dubey et al. 2007b; Mahmoudi et al. 2019). We found a western clade (with samples from western Europe and western Turkey; BP = 71; PP = 0.99), an eastern clade (with samples from Bulgaria, Romania, Georgia, and eastern Turkey; BP = 99%; PP = 1.0), and an Iranian clade (with samples from the Hyrcanian region of Iran; BP = 85%; PP = 1.0) (Figure 3; Dubey et al. 2007b; Mahmoudi et al. 2019). Within this tree, the Israeli C. leucodon cyt b sequences formed an isolated and strongly supported clade (BP = 91; PP = 1.0) (Figure 3). This Israeli clade is thus placed as a sister clade to the eastern clade, with (BP = 79; PP = 0.99). The average p-distance between sequences from Israel and sequences of the eastern clade was 0.023 (minimum = 0.019; maximum = 0.030). The results of the network analyses support the division into the four distinct C. leucodon lineages observed in the phylogenetic analysis (Supplementary Figure S2).

Figure 3:

Maximum likelihood tree of Crocidura leucodon cytb sequences. Phylogenetic relationships inferred from a matrix of 1,077 nucleotide positions for 63 individuals. Maximum likelihood bootstrap supports above 50% and Bayesian posterior probabilities above 0.70 are indicated near the corresponding node separated with a slash. Sequences obtained in this work are indicated in bold. Authors of sequence data in Supplementary Table S4.

3.3 Phylogenetic relationships within the C. suaveolens complex

The C. suaveolens complex is a highly diverse lineage. Dubey et al. (2007a) divided this complex, and its closest outgroup, C. shantungensis, into ten different clades, all of which were recovered in our phylogenetic reconstruction (Figure 4). All Israeli “C. suaveolens” sequences were found to be nested within clade V, ‘gueldenstaedtii’, of Dubey et al. (2007a) (BP = 100; PP = 1.0; Figure 4). This clade includes samples from western Asia, Corsica, Minorca, and Crete. The Israeli sequences did not form a distinct clade within clade V. The average p-distance between sequences from Israel and other sequences of the ‘gueldenstaedtii’ clade was 0.005 (minimum = 0; maximum = 0.018). Similarly the network analysis of the ‘gueldenstaedtii clade’ sequences did not separate Israeli haplotypes from Iranian, Georgian, Turkish, or Corsican haplotypes (Supplementary Figure S3).

Figure 4:

Maximum likelihood tree of Crocidura suaveolens complex cytb sequences. Phylogenetic relationships inferred from a matrix of 996 nucleotide positions for 215 individuals. Maximum likelihood bootstrap supports above 50% and Bayesian posterior probabilities above 0.70 are indicated near the corresponding node separated with a slash. Sequences obtained in this work are indicated in bold. Authors of sequence data in Supplementary Table S5.

3.4 Phylogenetic relationships within the S. etruscus clade

S. etruscus is a widespread Eurasian species whose distribution extends from France to Vietnam. Previous phylogeographic studies have shown that this species is probably paraphyletic and may be polyphyletic (Meegaskumbura et al. 2012a). Our phylogenetic tree (Figure 5) separated S. etruscus into two major clades: an eastern one and a western one. In agreement with Meegaskumbura et al. (2012a), the eastern clade included S. etruscus sequences from South Asia and S. madagascariensis sequences (BP = 77; PP = 0.99), supporting the view that S. madagascariensis is a junior synonym of S. etruscus. However unlike the findings of Meegaskumbura et al. (2012a), sequences from western European and the Near-East clustered together to form the Western clade (BP = 67; PP = 0.86). The two S. etruscus clades weakly cluster together (BP = 59; PP = 0.76) and are closely related to S. malayanus and S. fellowesgordoni sequences (BP = 100; PP = 1.0). One S. etruscus sequence from Vietnam (KF110756) formed a distant lineage, suggesting it could be a different species. Within our tree, the Israeli sequences did not form a distinct clade but, rather, clustered within the western clade (BP 67/0.86) together with specimens from France and Iran (Figure 5). The average p-distance between sequences from Israel and other sequences of the western clade was 0.003 (minimum = 0; maximum = 0.013). The results of the network analyses support the phylogenetic results (Supplementary Figure S4).

Figure 5:

Maximum likelihood tree of Suncus cytb sequences. Phylogenetic relationships inferred from a matrix of 1,140 nucleotide positions for 36 individuals. Maximum likelihood bootstrap supports above 50% and Bayesian posterior probabilities above 0.70 are indicated near the corresponding node separated with a slash. Sequences obtained in this work are indicated in bold. Authors of sequence data in Supplementary Table S6.

4.1 Crocidura ramona

Our work, based on sequences from type specimens, corroborates the view that C. ramona is a distinct species and probably endemic to Israel and the West Bank (Dubey et al. 2008; Ivanitskaya et al. 1996). C. ramona has been suggested to be related to the Palearctic “flat-headed rock-shrews” (i.e., Crocidura pergrisea, C.arispa, Crocidura armenica, Crocidura serezkyensis and C. zarudnyi) (Burgin et al. 2018b; Kryštufek and Vohralík 2001). Although the DNA of most of the latter has not yet been sequenced, our analysis indicates that C. ramona and C. zarudnyi are distantly related (Figure 2). Similarly, we have shown here that C. ramona is not closely related to another silvery-gray shrew – C. nana, as the two formed two distinct and distant clades in our phylogenetic analyses (Supplementary Figure S1). We note that there is some doubt regarding C. nana’s distribution. It is usually considered to be restricted to Somalia and Ethiopia (Cassola 2019; Hutterer 2005), whereas the samples sequenced in this work originated from Tanzania. We cannot be certain therefore that C. nana is the correct species assignment for the Tanzanian samples that we sequenced.

Because no morphological comparisons have been carried out between the skulls of C. katinka and C. ramona, the exact relationship between these two species remains to be determined. Although C. katinka has been suggested to be related to certain African species with cranial similarities (e.g., Crocidura bottegi, C. obscurior, Crocidura bottegoides) (Burgin et al. 2018b; Hutterer and Kock 2002), our analyses have demonstrated that C. ramona is unrelated to C. obscurior.

It is also possible that C. ramona is conspecific with C. portali. Thomas (1920) described C. portali as a small shrew (“though not excessively so”), with “pale drab-grey” pelage; and indicated that it has “clearly nothing to do with the C. russula group”. The gross morphology of the holotype (BMNH #19.4.11.9) and its size agree with C. ramona – though a detailed examination is still needed in order to confirm or refute this. Kryštufek and Vohralík (2001) suggested that C. portali may be a valid species, related to Crocidura arispa (and other members of the pergrisea group), and a senior synonym of C. ramona. Hutterer and Kock (2002) indicated that C. ramona and C. portali are similar in skull dimensions as well as in pelage, but note that they are also similar to C. gmelini, a species that they note requires “a better definition”. C. gmelini (ranging from Iran to Mongolia) has been synonymized with C. suaveolens, a species distantly related to C. ramona (Figure 2). C. arispa, C. ramona, and C. katinka are currently all considered valid species, while C. portali is not (Burgin et al. 2018b; IUCN 2020). Clearly, a taxonomic revision of these taxa is warranted. Neither the DNA of the C. portali type, nor that of any shrews identified as C. arispa, C. pergrisea, or C. katinka, has been sequenced and, unfortunately, we failed to amplify any cyt b fragment from a tissue of a specimen identified as C. portali (BMNH ZD 1971.817). We thus tentatively ascribe the sequence we obtained to C. ramona, pending a taxonomic revision.

4.2 Crocidura leucodon

The phylogenetic reconstruction indicates that the Israeli C. leucodon forms a distinct clade related to the eastern clade (Figure 4, Dubey et al. 2007b). The different C. leucodon clades have been suggested to have diverged during the Pleistocene glaciations (Dubey et al. 2007b; Mahmoudi et al. 2019). It is unlikely however that the Israeli clade corresponds to a fourth refugium from the Ice Age. Rather, the observed mitochondrial genetic divergence is possibly the result of the edge position of the Israeli population at the southernmost part of the C. leucodon range. In agreement, Mendelssohn and Yom-Tov (1999) noted that C. leucodon is less abundant than C. s. gueldenstaedtii, basing their conclusion on the number of shrews deposited in museum collections. The current data available at the Steinhardt Museum of Natural History support the view that the range of C. leucodon is more limited than that of C. suaveolens in Israel, since C. leucodon samples are rarer in the collection (79 leucodon vs. 594 suaveolens), and with one exception, restricted to the Mediterranean biome (Figure 1B). Edge populations have been suggested to harbor adaptive traits and genetic variability that may be important when considering future conservation needs (Hampe and Petit 2005; Mátýas et al. 2009). We thus also recommend that further population and genomic studies be carried out on this species, since its population status is currently unknown (Shenbrot et al. 2016), and assessment from museum collections alone may provide a biased representation of the population status.

4.3 Crocidura suaveolens gueldenstaedtii

The phylogenetic analysis indicated that all Israeli “C. suaveolens” sequences are part of clade V – the ‘gueldenstaedtii’ clade (Dubey et al. 2007a), and close to the Turkish and Georgian sequences (Figure 5). While the Israeli populations represent the edge of this species complex range, they do not present a mitochondrial genetic difference from the rest of its clade, unlike the situation in C. leucodon. This may be due to their larger population size and wider range in Israel compared to that of C. leucodon (Figure 1B). However, it is also possible that genetic differences exist in the nuclear genome. Based on the high similarity between the Israeli and Balkan sequences we agree with the status of “least concern” listed by the IUCN (Palomo et al. 2016). Our molecular cyt b data unequivocally place Israeli ‘suaveolens’ specimens within the ‘gueldenstaedtii’ clade. The specific status of the form awaits a thorough taxonomic revision, including an examination of animals from the type localities, preferably alongside specimens from the type localities of other members of the complex (e.g., C. ‘gmelini’, C. suaveolens monacha, C. portali) as well as phylogenetic information from nuclear markers. Until such a study is carried out, we tentatively refer to Israeli specimens as members of C. suaveolens gueldenstaedtii.

4.4 Suncus etruscus

The Israeli S. etruscus sequences appear to be closely related to the European ones and one Iranian sequence. The differences between the Israeli sequences and those of other western clade members are very small and comparable to the distances observed within the C. suaveolens gueldenstaedtii clade. This suggests that the western clade encompasses individuals from France in the west to Iran in the east and Israel in the south. However, this species has been poorly sampled in molecular studies to date. As a case in point, only a few S. etruscus specimens have been sequenced from Iran where members of both the western and eastern clade are present (Darvish et al. 2017; Ohdachi et al. 2016). More data are needed in order to decipher the population structure of this species.