Ancient Divergence Driven by Geographic Isolation and Ecological Adaptation in Forest Dependent Sundaland Tree Squirrels

Introduction

The small brown tree squirrels of the genus Sundasciurus are endemic to Sundaland and the Philippines. This genus was described by Moore (1958) based on the presence of a Y-shaped trans-bullar bony septum and includes species previously assigned to Sciurus and Tomeutes (Thomas, 1915; Chasen, 1940). Most species were widespread across multiple landmasses in Sundaland. This genus has received more attention than other genera in the subfamily (Callosciurinae/Nannosciurinae). Deep genetic divergences associated with habitat (such as highland or lowland forest) or geography (such as land mass) were obscured by conservative morphology and have resulted in multiple revisions of the genus in the last few years, increasing the number of species from 11 to 17 (Heaney, 1979; Den Tex et al., 2010; Thorington et al., 2012; Hawkins et al., 2016a). These latest systematic revisions have largely focused on high elevation ‘populations’ of widespread taxa which in all cases turned out to be very divergent and merited specific status. The recognition of these species was supported by both genetic and morphological lines of evidence. The genus was divided into two subgenera by Moore (1958). These subgenera, Aletesciurus and Sundasciurus, were defined based on morphological criteria (presence/absence of a sagittal crest, an inconspicuous/inflated antero-mesial lobe, and skull size), but these were not found to be reciprocally monophyletic with molecular data (Den Tex et al., 2010; Hawkins et al., 2016a). At that time, taxonomic rearrangements were not performed because of limited sampling and low support of key nodes in the tree, suggesting that additional data were necessary. The subgenus Aletesciurus contains one widely distributed species (Sundasciurus hippurus of the Malay Peninsula, Sumatra, and Borneo), a narrowly distributed Bornean mountain endemic (S. everetti), and nine species distributed in the Philippines (S. samarensis, S. philippinensis, S. mindanensis, S. rabori, S. steerii, S. moellendorffi, S. davensis, S. juvencus, and S. hoogstraali) (Chasen, 1937; Den Tex et al., 2010; Hawkins et al., 2016a). The subgenus Sundasciurus includes seven species restricted to Sundaland in two clades: a mainly highland clade on a long branch (S. altitudinis, S. brookei, S. jentinki, S. tahan, and S. tenuis) and a widespread lowland clade sister to the other two (S. fraterculus and S. lowii) (Thorington and Hoffmann, 2005; Den Tex et al., 2010; Thorington et al., 2012; Hawkins et al., 2016a). Many of these widespread species contain multiple subspecies described from a limited number of specimens that were generally identified by differences in pelage coloration (Corbet and Hill, 1992; Thorington et al., 2012). In order to evaluate the diversity in the genus and uncover other unidentified lineages, we used populations as the unit of analysis.

Despite recent research on the phylogenetics of some of the members of the genus Sundasciurus, less attention has been paid to lowland Sunda squirrels. Widespread across Sundaland, lowland species of rodents have been shown to have more and stronger genetic structure associated with land masses, as compared to larger mammals but these divergences have not yet led to taxonomic revisions. Given the forest-dependent nature of these squirrels, an improved knowledge of their evolutionary history is key for our understanding of the dynamics of Sundaland rainforests (Den Tex et al., 2010). Deep divergences and speciation events may be the consequence of past isolation of landmasses or forest pockets while lack of genetic structure could shed light on recent episodes of land or forest connection. Conflicting patterns have been shown across different phylogeographic studies on mammals distributed across Sundaland (Mason et al., 2011; Leonard et al., 2015). Some species of murids, carnivores, suids, greater mouse deer and pangolins have weak or no genetic structure across Sundaland and/or within Borneo (Gorog et al., 2004; Achmadi et al., 2013; Patel et al., 2017; Mason et al., 2019; Veron et al., 2020) while others, such as Sunda rats, viverrids, colugos or the lesser mouse deer, show much deeper divergences some of which predate to the Pleistocene (Den Tex et al., 2010; Mason et al., 2011; Camacho-Sánchez, 2017; Veron et al., 2019). Recently, Husson et al. (2019) provided substantial evidence for a new paleogeographic scenario in which prior to 400 thousand years ago (ka) the Sunda shelf would have been continuously exposed. According to this study, dispersal across this exposed shelf would have been possible through this period for terrestrial species. Nevertheless, geological connectivity does not necessarily imply population connectivity, different dates of divergences between populations on different land masses might reflect different dispersal capabilities across non-forested ecological barriers such as the possible “savanna corridor” in interior Sundaland (Heaney, 1991; Sheldon et al., 2015).

Within the genus Sundasciurus, previous phylogenies based on mitochondrial cytochrome b sequences have identified deep divergences between lineages within Low’s squirrel (S. lowii), and a lack of reciprocal monophyly between Low’s squirrel and the Fraternal squirrel (S. fraterculus) from the Mentawai islands (Den Tex et al., 2010; Hawkins et al., 2016a). Low’s squirrel is a good example of a widespread, lowland forest dwelling tree squirrel. It is currently listed as a Species of Least Concern by the IUCN, largely based on its wide historic distribution (Meijaard, 2016). However, the lowland tropical forest is rapidly disappearing in this region due to human induced habitat modifications and its current status and distribution need to be re-evaluated. Furthermore, it is possible that the widely distributed Low’s squirrel is actually a complex of multiple species each with a more restricted range. Low’s squirrel contains seven nominal subspecies distributed across Borneo (S. l. lowii), northern Borneo’s offshore islets Banguey and Balambangan (S. l. bangueyae), the Malay Peninsula (S. l. robinsoni), southern Malay Peninsula Riau offshore islands (S. l. seimundi), Sumatra (S. l. humilis), west Sumatra Batu offshore islands (S. l. balae), and the Natuna islands (S. l. natunensis) (Thorington and Hoffmann, 2005; Thorington et al., 2012) all based on morphology (Figure 1).

In this study, we use genetic data including complete mitochondrial genomes and twelve nuclear loci from across the Sundasciurus lowii and S. fraterculus distribution to construct a well-resolved phylogeny of the group, and to date divergences. These results are placed in the context of the morphological variation and distribution of the different taxa in order to re-evaluate and resolve the number of species in this group. Additionally, a comprehensive taxonomic evaluation of this species complex will have important conservation implications that could lead to the re-assessment of the IUCN Red List category for some of these squirrels based on their updated distributional ranges. It will also provide a better understanding of the ecological and evolutionary processes that led to the high diversification rates of these forest dependent squirrels in this region.

Materials and Methods

Materials

A total of 11 samples of Sundasciurus lowii and S. fraterculus from the major landmasses and isolated archipelago populations (Malay Peninsula, Borneo, Mentawai Islands and Natuna Islands) and 29 outgroups were included in the molecular part of the study (Supplementary Table S1). Populations from Sumatra and satellite islands (Batu and Riau) were not sampled, because only type specimens are available (S. l. humilis, S. l. balae, S. l. seimundi). In any case, different studies have suggested recent connections among Sumatra and Peninsular Malaysia (Ruedi and Fumagalli, 1996; Mason et al., 2011; Camacho-Sánchez et al., 2017) and among these satellite islands and the nearest landmasses (Mason et al., 2019) so our geographic coverage could potentially represent most of the genetic structure of this complex. Twelve samples were derived from tissue collected from historical specimens housed in three museums, 10 samples were modern tissue coming from museum loans and 14 sample sequences were downloaded from Genbank. Four samples were obtained from specimens collected during a previous expedition to Borneo (Camacho-Sánchez et al., 2019) and these were all taken according to the guidelines of the American Society of Mammalogists (Sikes, 2016) and as approved by Institutional Animal Care and Use Committees (Smithsonian Institution, National Museum of Natural History, proposal number 2012-04 and Estación Biológica de Doñana proposal number CGL2010-21424). Field research was approved by Sabah Parks [TS/PTD/5/4 Jld. 45 (33) and TS/PTD/5/4 Jld. 47 (25)], the Economic Planning Unit (100-24/1/299) and the Sabah Biodiversity Council (Ref:TK/PP:8/8Jld.2). All recognized species and subspecies of this complex were included in the morphological analyses (Figure 1 and Supplementary Table S2).

DNA Extraction, Amplification, Enrichment and Sequencing

DNA was extracted using phenol-chloroform with ethanol precipitation. Amicon Ultra 0.5 mL filters were used for the final clean instead of ethanol precipitation for historic samples. All museum samples were processed in an isolated ancient DNA laboratory following strict protocols to control for contamination.

Twelve nuclear loci (introns) previously found to be present in single copies and informative in closely related mammalian species (Igea et al., 2010) were targeted. Primer sets for four of those introns developed for Tamiasciurus squirrels by Chavez et al. (2014) worked in our taxa as a multiplex. For the other eight loci, we designed primers following the steps of Forcina et al. (2020). We developed 15 primer sets to amplify smaller fragments (ranging from 88-233 base pairs) of these same loci for the low quality specimens (Supplementary Table S1). PCR conditions for the two modern DNA and four historic DNA multiplexes are in Supplementary Table S3. PCR products were purified using streptavidin magnetic ‘SPRI’ beads (Rohland and Reich, 2012) and visualized on agarose gels. Given the degraded nature of historic DNA, three replicate PCRs were run per sample. Amplicon libraries were cleaned and checked on a gel with reference standards for quantification (Image Lab 5.2.1 software), and were pooled at equimolar concentrations. Illumina TruSeq adapters with an individual index were added to the purified PCR products in a second PCR as performed in Camacho-Sánchez et al. (2017).

We used target capture to recover complete mitochondrial genomes. We used an Illumina modification to the protocol based on Maricic et al. (2010) to enrich complete mitogenomes from historic samples. Enrichment conditions followed kit protocols as described in Hawkins et al. (2016b). Briefly, samples were denatured, incubated for ∼24 h, recovered the baits, and amplified from the p5/p7 priming site on the Illumina adapters for a limited number of PCR cycles to generate concentrations of mitochondrial DNA which could be sequenced. PCR products were purified as above.

Dual-indexed shotgun libraries were sequenced on an Illumina MiSeq using 2 × 300 (600 cycle) PE at the Center for Conservation Genomics, Smithsonian Conservation Biology Institute.

Preprocessing and Quality Scanning of Sequencing Data

Adaptor removal and quality trimming was performed with Trimmomatic (Bolger et al., 2014) with the sliding window parameter set to 5:20, read minimum length parameter to 50 bp and leading and trailing to 5. Single quality scans were run on the raw fastq files before and after trimming with FastQC (Andrews, 2010) and MultiQC (Ewels et al., 2016) was run to create a single report visualizing the output from all the samples of each run, enabling global trends and biases to be quickly identified. Mitochondrial genome reads were preprocessed following the protocol described in Hawkins et al. (2016b).

Assembly and Alignment of Mitochondrial Genomes

Quality trimmed reads were mapped to a reference mitochondrial genome with BWA-MEM algorithm or with Geneious Mapper iterative algorithm (up to 10 iterations and medium-sensitivity) for some divergent samples. The output of BWA was converted to BAM files with SAMtools (Li et al., 2009). These were later sorted, merged and PCR duplicates removed. Finally, BAM files were imported to Geneious where consensus sequences were called (minimum 5x and 75% threshold). Samples mapped in Geneious were inspected for duplicate removal with Dedupe plugin of the BBTools package v. 35.82 (Bushnell, 2015). We aligned our sequences with Callosciurinae whole mitochondrial sequences available in Genbank (Supplementary Table S1) with MAFFT v7.450 Geneious plugin automatic algorithm (Katoh and Standley, 2013) under default parameters (scoring matrix:200PAM/K = 2, gap open penalty = 1.53, offset value = 0.123). The control region was removed from the mitogenome assemblies because it was poorly assembled in many historic samples and has been shown to provide low phylogenetic resolution and overestimation of divergence times (Duchêne et al., 2011). ND6 is present in the light strand and was thereby reverse complemented. Protein-coding regions were translated and inspected for frameshift mutations and for the presence of unexpected stop codons to prevent inclusion of NUMTs. We computed summary statistics of the alignment with AMAS (Borowiec, 2016).

Genotyping and Alignment of Nuclear Data

For nuclear sequences, we imported the trimmed reads to Geneious and mapped these to the closest homologous sequences found in Genbank (Marmota marmota (NW_015351277.1, NW_015351233, XM_027934759, NW_015351257, NW_015351271, NW_01535121) and Tamasciurus hudsonicus (KF883544, KF883925, KF885067) with Geneious Mapper (low-medium sensitivity and up to 5 iterations). We called consensus sequences in Geneious with a minimum 5× coverage and 75% consensus threshold. Each locus was aligned independently with MAFFT v7.450 Geneious plugin automatic algorithm (Katoh and Standley, 2013) under default parameters (scoring matrix:200PAM/K = 2, gap open penalty = 1.53, offset value = 0.123). We translated the alignments to amino acids and inspected them for insertions, deletions and premature stop codons to prevent inclusion of paralogous sequences. Amplicon sizes exceeded the length generated by the Illumina sequencing for two loci (GDAP1, 795 bp and P4HA2, 758 bp). This generated missing data in the center of those alignments (203 and 295 bp), which were trimmed.

All nuclear sequences were resolved into statistically probable haplotypes using PHASE 2.1.1 (Stephens, 2005) with an acceptance threshold of 0.90. The online application SeqPHASE (Flot, 2010) was used to convert FASTA files to PHASE input files, as well as to convert PHASE output back to FASTA format. Due to the different size of modern and historic sample sequences, the initial gene alignments contained sections of missing data for the shorter historic samples. Since PHASE can not deal with missing data in heterozygous sites, two alignments (historical and modern samples) were generated for each gene to remove most missing data and overcome this problem. Genotypes that were still unresolved following PHASE were phased manually based on the original BAM alignment (Camacho-Sánchez, 2017).

Phylogenetic Analyses and Haplotype Networks

We performed a likelihood ratio test (LRT) with MEGA 7 to test the clock-like behavior of the mitochondrial sequences. A strict molecular clock was rejected. The optimal partitioning scheme and substitution model was selected by PartitionFinder 2.1.1.2 (Lanfear et al., 2016). We followed Mackiewicz et al. (2019) considering 3 codon positions for each individual protein coding gene and separate partitions for each of the RNA genes, but retained contiguous t-RNAs because individual t-RNA alignments would be too short to be informative (Nicolas et al., 2019). This reduced the number of initial partitions from 67 to 53. We used the “greedy” algorithm with branch length estimated as “unlinked”, AICc criteria and “–raxml” command line option to search for the best-fit partitioning scheme (11 partitions) and substitution model (GTR + GAMMA). We tested for saturation plotting transitions and transversions of each alignment partition against each other and against raw/uncorrected pairwise genetic distances with the APE R package and “dist.dna” command. We also performed Xia substitution saturation tests in Dambe 7 (Xia, 2018) for the 3rd codon of some fast-evolving genes (e.g., Cytochrome b) as a complementary line of evidence. We did not find saturation so all partitions were kept in downstream analyses. We performed phylogenetic inference through a maximum likelihood framework with RAxML (Stamatakis, 2014) on the mitochondrial genome dataset. We specified the output of PartitionFinder2 as input for RAxML (Stamatakis, 2014). We selected the rapid bootstrapping algorithm, GTR GAMMA substitution model and followed the extended majority rule stopping criterion (-autoMRE). The maximum likelihood phylogeny was used for the BEAST analyses.

Divergence times were estimated with our mitochondrial dataset in BEAST2 (Bouckaert et al., 2014). The mitochondrial genome alignment was split into the best partition schemes as identified by Partitionfinder2 with AMAS. We performed site modeling through a Bayesian approach with bmodeltest, selecting mutation rate estimation and empirical frequencies priors (Bouckaert and Drummond, 2017). Following MEGA7 tests, a lognormal (uncorrelated) relaxed clock was used under a Yule speciation tree, with linked partitions and default operators. Only one sequence per putative species was kept to comply with Yule speciation prior assumptions. Three independent runs of Markov chains (MCMC) for Monte Carlo simulations were run for 50,000,000 generations, with parameters and trees sampled every 5000 generations. Convergence was checked using Tracer 1.7.1 (Rambaut et al., 2018). For each run, the first 10% of sampled trees were discarded as burn-in. Two fossil constraints were selected to calibrate the phylogeny. A fossil of Callosciurus sp. (specimen number: YGSP 21682; left upper third molar that conforms to the proportions of living Callosciurus) from locality Y589 (Chinji formation, Siwaliks, Pakistan) dated by paleomagnetic correlation to 14 Ma (Flynn, 2003) for the Callosciurus and Sundasciurus clade; and a stem form of Tamiops (specimen numbers: Z2438, Z2439, Z2440, Z2441) from Z122 (Vihowa Formation) dated by paleomagnetic correlation to 18.7 Ma (L. Flynn, pers. comm.) as the oldest known fossil for the Callosciurus, Sundasciurus, Tamiops, and Dremomys clade. All constraints were set using hard minimum bounds and soft upper bounds using a lognormal prior, as suggested by Parham et al. (2012).

We performed TCS haplotype networks with popART (Leigh and Bryant, 2015) separately, for the modern and historic alignments of the mitochondrial genomes and nuclear intron sequences. We computed uncorrected pairwise genetic distances on the Cytochrome b gene with the APE R package and dist. dna command (Paradis et al., 2004). Results were visualized with ade4:table.paint. Finally, a multidimensional scaling analysis was performed for both datasets with stats:cmdscale and ggplot2. We defined the a priori groups for the species tree analysis within our ingroup following four criteria: Cytochrome b uncorrected genetic distances > 5% (Baker and Bradley, 2006) reciprocal monophyly in mitochondrial genome tree, clusters in MDS, presence of private alleles in the intron haplotype networks and phenotype differentiation. Based on these criteria we considered six groups within our ingroup as input for the ^∗BEAST species tree: Sundasciurus fraterculus, S. lowii robinsoni, S. l. natunensis, S. l. lowii from Sabah, S. l. lowii from Sarawak, and S. l. lowii from East Kalimantan.

We included all intron phased alleles available per group in the ^∗BEAST analysis, and unlinked partitions, selected an uncorrected lognormal clock, bmodeltest with empirical frequencies and mutation rate estimation, analytical population size integration model and Yule speciation prior. Three independent runs were conducted for 50000000 generations Markov chains (MCMC) for Monte Carlo simulations, with parameters and trees sampled every 5000 generations.

Species delimitation was conducted using the program BPP X1.2.2 (Yang, 2015) details regarding this analysis are specified in the Supplementary Material File 1.

Morphological Data Collection

The specimens examined here are housed in the following natural history repositories: American Museum of Natural History, New York, United States (AMNH); Natural History Museum, London, United Kingdom (BMNH); Estación Biológica de Doñana, Seville, Spain (EBD); Field Museum of Natural History, Chicago, Illinois (FMNH); Museum of Comparative Zoology, Harvard University, Cambridge, Massachusetts (MCZ); and National Museum of Natural History, Smithsonian Institution, Washington, DC, United States (USNM) Supplementary Table S2. Skins were carefully examined to avoid any misidentifications and a geographically balanced subset was photographed for further evaluation. During our initial surveys of the squirrel collections, we found that a high proportion of the skulls were broken and therefore decided to undertake traditional morphometric analyses instead of geometric morphometrics. Dental eruption and wear patterns were checked to age each specimen. Fusion of the presphenoid-basisphenoid and basisphenoid-basioccipital sutures were not informative to age individuals in this taxon since some large-sized adults with very worn molars had non-fused sutures. Only adults were included in the morphometric analyses. We combined data from males and females in the morphometric analyses because sexual dimorphism is uncommon in tree squirrels (Moore and Tate, 1965; Tenzin et al., 2013; C̀anády et al., 2015). However, if certain characters were sexually dimorphic, results would be interpreted in the context of phenotypic variation among putative species inclusive of potential patterns of sexual dimorphism, thereby providing a more conservative assessment of divergence (Meik et al., 2018). Seventeen cranial measurements were taken with a Fowler High Precision electronic digital caliper to the nearest 0.01 mm by AH (Supplementary Table S2). These were occipitonasal length (ONL), zygomatic breadth (ZB), interorbital breadth (IB), length of nasals (LN), anterior nasal breadth (ANB), length of rostrum (LR), height of rostrum (HR), mastoid breadth (MB), height of braincase (HBC), length of orbit (LO), length of diastema (LD), length of bony palate (LBP), breadth of bony palate at fourth premolar (BBP), length of auditory bulla (LB), crown length of maxillary cheek teeth (CLPMM), median palatal length (MPL) and least breadth of caudal point of zygomatic process of frontal bone (LBC). All measurements are as defined in Musser et al. (2010) except median palatal length (MPL) and least breadth of caudal point of zygomatic process of frontal bone (LBC), which are defined in Hayashida et al. (2007). In addition, given the curved rostrum of our target taxa, the height of rostrum (HR) was modified after conducting preliminary error measuring tests. We define it as the distance from dorsal surface of the nasals to ventral surface of the premaxillary at the premaxillary-maxillary suture and taken perpendicularly to the molar row axis. External measurements, head-body length (HB), tail (T), hindfoot (HF), ear (E) and weight (W), were taken directly from specimen labels (Supplementary Table S2). We examined, measured and photographed all of the types of our target taxa.

Morphometric Statistical Analyses

We performed standard univariate descriptive statistics (mean, standard deviation and observed range) for each of the subspecies (Thorington and Hoffmann, 2005) and also for each of the S. lowii Borneo populations in a separate analysis (Supplementary Table S4). We log-transformed each measurement to its standard deviation prior to computing the principal component analyses (PCA) so that the data was analyzed on the basis of correlations instead of covariances. PCA was implemented by the prcomp command in R studio (Allaire, 2012; R Core Team, 2018). Results were extracted and visualized with the following functions of the factoextra package (Kassambara and Mundt, 2017): fviz_pca_ind, plots PCA results; fviz_pca_biplot, biplot of individuals and variables; get_eigenvalue, extracts eigenvalues, variance percentage and cumulative variance percentage; get_pca_var, outputs each variable’s contribution to variance.

We conducted a discriminant analysis of the principal components (DAPC) with the dapc function of adegenet (Jombart and Ahmed, 2011). Identification of the clusters was made with find.clusters. This function first transforms the data using PCA, asking the user to specify the number of retained principal components (PCs) interactively. Then, it runs k-means algorithm with increasing values of k and computes associated summary statistics (by default, BIC). We performed a cross-validation with the xvalDapc() function with default parameters to decide the number of PCs to retain. This function runs 30 replicates of cross-validation for a number of PCs less than the total number of variables. After having an approximate idea of the PCA axes range to be retained, we re-ran the cross-validation with 1000 replicates. Different number of clusters had similar BIC values so a DAPC was run for each scenario (3-10 clusters). Different DAPCs (retaining 13-3 PCs and 6-3 discriminant functions) with a priori grouping for each of the subspecies (Wiley, 1978; Thorington and Hoffmann, 2005) and for each of the S. lowii Borneo populations (in a separate analysis) were also performed. We interpreted group memberships with the functions summary, compoplot and assignplot to test how clear-cut clusters were. These functions indicate the proportions of successful reassignment (based on the discriminant functions) of individuals to their original clusters. To increase the number of samples (from 55 to 152) and include the types, we made bivariate plots of craniodental variables based on the output of the DAPC and standard external measurements.

Registration of Nomenclature

The electronic version of this article in Portable Document Format (PDF) will represent a published work according to the International Commission on Zoological Nomenclature (ICZN), and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser at https://www.urn:lsid:zoobank.org:pub:727A0DA1-2B92-4400-A3C2-CE816BE00B1F

Results

Whole mitochondrial genomes were constructed from 24 individuals, and their alignment (without control region) was 15735 bp, with 5557 parsimony informative sites and 8.97% missing data. Missing data was mainly in the sequences from the eleven historic samples (15–47%). Mitochondrial genomes of all ingroup taxa were successfully assembled, although only 53% of one Sundasciurus lowii natunensis from two attempted samples could be reconstructed. Similarly, Sundasciurus lowii natunensis was the only taxon from which we failed to amplify any introns. Three intron loci were not informative or did not amplify correctly (GAD2, SLC17A9, Ulk1b) and were excluded from downstream analyses.

Mitochondrial Genome Phylogenetic Relationships and Divergence Dating

All nodes were strongly supported (Bootstrap support (BP) > 84/Posterior probability (PP) > 0.98) in the mitochondrial genome maximum likelihood tree (ML) and Bayesian maximum clade credibility tree (BMCC) and the topology was consistent between approaches (Figures 2, 3). Only the relationships among S. l. lowii had weak support in the ML tree (BP = 0.66) but strong support in the BMCC (PP = 0.99). These tree topologies were consistent with previous studies (Den Tex et al., 2010; Hawkins et al., 2016a). The Low’s squirrel complex contained two divergent clades from east and west Sundaland. The former included S. l. lowii from Borneo and S. l. natunensis from the Natunas, while the latter had S. l. robinsoni from peninsular Malaysia and S. fraterculus from Mentawai Islands. The sister relationship of S. fraterculus with S. l. robinsoni makes S. lowii sensu lato paraphyletic as in Den Tex et al. (2010). Within S. l. lowii, the population from Sarawak was sister to the clade containing the ones from Sabah and East Kalimantan. All genera (Tamiops, Dremomys, Callosciurus and Sundasciurus) were recovered as monophyletic. Previously suggested subgenera within Sundasciurus (Aletesciurus and Sundasciurus; Moore, 1958) were not found to be reciprocally monophyletic. Previous studies either could not resolve the lower nodes or found them not to be monophyletic (Den Tex et al., 2010; Hawkins et al., 2016a) and thus had suggested suspending that nomenclature (Hawkins et al., 2016a).

The time to the most recent common ancestor (TMRCA) of Callosciurinae and Sundasciurus was estimated at 26.86 (95% HPD: 20.00–33.85) and 12.25 (95% HPD: 10.67–13.67) million years ago (MYA). These dates are older than those previously estimated by Mercer and Roth (2003) who associated the diversification of Sundaland genera with the lowest pre-Pleistocene sea stand in the Cenozoic about 11 MYA. TRMCA for the subgenus Sundasciurus and Sundasciurus lowii complex were estimated at 11.65 (95% HPD: 10.09–13.11) and 7.01 (95% HPD: 5.74–8.33) MYA. At 4.41 (95% HPD: 3.39–5.53) MYA one of Sundasciurus lowii complex lineages split into the ancestors of S. l. lowii clade and S. l. natunensis. At 3.20 (95% HPD: 2.30–4.12) MYA the other Sundasciurus lowii complex lineage split into the ancestors of S. fraterculus and S. robinsoni. The lowland Sundasciurus tenuis and highland S. altitudinis exhibited remarkable deep levels of intraspecific divergence dated at 3.96 (95% HPD: 2.82–5.05) and 2.85 (95% HPD: 2.04–3.79) MYA. These dates are older to what was previously estimated from Cytochrome b by Den Tex et al. (2010).

Divergence and Private Alleles

Many subspecies had private alleles (Figure 4 and Supplementary Figure S1). The taxa S. l. robinsoni with S. fraterculus and S. l. lowii had private alleles for all introns. Within Borneo S. l. lowii of Sabah, S. l. lowii of Kalimantan and S. l. lowii of Sarawak had private alleles for seven out of nine introns. The narrow endemic S. fraterculus had private alleles for three of nine introns.

Cytochrome b uncorrected genetic distances among some of the taxa (S. l. robinsoni, S. l. natunensis, S. l. lowii of Sabah and S. l. lowii of Sarawak) were above 5% (Supplementary Figure S2), a scale of divergence which may suggest a revision of the status of those taxa could be of interest (Baker and Bradley, 2006). Genetic distances among these taxa were higher than that observed among other recognized species (Supplementary Figure S2). These results further highlight that while Cytochrome b distances might indicate potential unrecognized genetic species, it should never be used alone to define species. We note that high percent sequence divergence levels were also observed between S. t. tenuis and S. t. parvus, and among populations within S. hippurus, S. mindanensis, and S. altitudinis (Supplementary Figure S2), suggesting that additional work on these groups in the future could shed more light on the evolution of this genus.

Multilocus Species Tree and Species Delimitation

Phylogenetic relationships inferred by the intron species tree constructed in ^∗BEAST were concordant with the mitochondrial genome tree topology. Only paraphyly of the subgenus Sundasciurus in the nuclear species tree was inconsistent with the mitochondrial topology, but concordant with Hawkins et al. (2016a)Figure 3. However, deep nodes such as the MRCA of S. everetti and S. hippurus or shallow ones such as the relationships among Bornean S. lowii were weakly supported (PP = 0.62-0.63) while the other nodes were strongly supported (PP > 0.99).

Our analyses of nuclear data were consistent across the independent BPP runs and among the five prior combinations. We inferred a posterior probability of 1.0 for all species postulated by the ^∗BEAST guide tree except S. fraterculus and S. l. lowii from E. Kalimantan. The former was not supported (PP = 0.07–0.26) and the latter went from weakly to moderately supported (PP = 0.51–0.95). The most likely scenario in every run supported all species except S. fraterculus (PP = 0.41–0.88).

Morphometric Results

A high degree of morphometric differentiation among the different populations was observed (Figures 5, 6) with a large part (82%) of the variance explained by the first principal component. PC1 exhibited different loadings for most variables, with crown length of maxillary cheek teeth (17%), length of nasals (14%) and breadth of bony palate (10%) driving most of the variation. PC1 discriminated the larger-sized subspecies S. l. lowii and S. l. bangueyae from the smaller-sized subspecies S. l. robinsoni, S. l. humilis, S. l. balae and S. fraterculus. PC2 explained 5.7% of the variance and was mainly correlated with shape robustness/broadness, bulla and molar row length. The variables with a higher contribution were crown length of maxillary cheek teeth (26%), interorbital breadth (21%), anterior nasal breadth (20%), length of bulla (11%), and least breadth of caudal point of zygomatic process of frontal bone (7%). It discriminated the relatively gracile and large molar S. l. robinsoni, southern Natuna S. l. natunensis, east and west Kalimantan S. l. lowii and S. l. bangueyae, from the relatively robust/broad shaped and small molar S. fraterculus, S. l. balae, S. humilis, Sarawak and most Sabah samples of S. l. lowii, and northern Natuna S. l. natunensis. PC3 explained 3.3% of the variance and was correlated with anterior nasal breadth (35%), crown length of maxillary cheek teeth (14%), length of bulla (10%) and mastoid breadth (9%). It discriminated the narrow-nasal shaped S. l. robinsoni, S. l. balae, S. humilis, northern Natuna S. l. natunensis and east and west Kalimantan S. l. lowii from broad-nasal shaped S. fraterculus, S. l. bangueyae, southern Natuna S. l. natunensis and Sarawak and east-central Kalimantan S. l. lowii. Sabah S. l. lowii fell in between.

The DAPC discriminated most subspecies (Supplementary Figure S3). Some overlap was observed between S. l. balae and S. l. robinsoni, and between S. l. lowii and S. l. bangueyae. The Mount Kerinci sample of S. l. humilis (belonging to the junior synonym S. l. vanakeni) had an intermediate position between S. fraterculus and S. l. robinsoni, while the only west Bornean sample of S. l. lowii was clustered between other Borneo S. l. lowii and S. l. natunensis. Only samples from East Kalimantan could be discriminated from the remaining Bornean populations. This analysis identified crown length of maxillary cheek teeth, bony palate breadth, orbit length, rostrum length, interorbital breadth, and bulla length, as the measurements most useful for distinguishing among all taxa (Supplementary File S1).

Taxonomic Revision

Sundasciurus robinsoni (Bonhote, 1903) – Robinson’s squirrel

Sciurus robinsoni Bonhote (1903):24.—Type locality. Bukit Besar, Nawngchik, Patani, Peninsular Thailand, 2500 ft.

Holotype. BMNH 3.2.6.55, skull, skin, adult female, collected 30 August 1901 by C. H. Robinson and N. Annandale

Sciurus balae Miller, 1903:14.—Type locality. Tana Bala, Batu Islands, Sumatra

Sciurus piniensis Miller, 1903:14.—Type locality. Pulo Pinie, Batu Islands, Sumatra

Sciurus seimundi Thomas and Wroughton, 1909:440.—Type locality. Kundur Island, Riau Islands

Sciurus humilis Miller, 1913:24.—Type locality. Kateman river, East Sumatra

Sciurus vanakeni Robinson and Kloss, 1916:270.—Type locality. Barong Bharu, Korinchi, Sumatra, 4000 ft.

Subspecific Taxonomy

Sundasciurus robinsoni (Bonhote, 1903)

Sciurus robinsoni alacris Thomas, 1908e: 306.—Type locality. Selangor-Pahang border, Malaya, 3000ft.

Sciurus seimundi Thomas and Wroughton, 1909:440.—Type locality. Kundur Island, Riau Islands

Sundasciurus robinsoni balae (Miller, 1903)

Sciurus balae Miller, 1903:14.—Type locality. Tana Bala, Batu Islands, Sumatra

Sciurus piniensis Miller, 1903:14.—Type locality. Pulo Pinie, Batu Islands, Sumatra

Sciurus humilis Miller, 1913:24.—Type locality. Kateman river, East Sumatra

Sundasciurus robinsoni vanakeni (Robinson and Kloss, 1916)

Sciurus vanakeni Robinson and Kloss, 1916:270.—Type locality. Barong Bharu, Korinchi, Sumatra, 4000 ft.

Emended Comparison With Other Sundasciurus Species

The species Sundasciurus robinsoni has a dorsum that ranges from medium brown with orange agouti to dark brown (S. r. vanakeni), and its venter ranges from white to pale yellow/buff white, with a reduction in the extent of this pale coloration and a lack of distinct margins in the case of S. r. vanakeni. Some populations (S. r. balae and S. r. vanakeni) have a grayish ventral coloration in limbs while others do not (S. r. robinsoni). It can be easily distinguished from other medium-sized western Sundaland Sundasciurus based on its ventral coloration and tail. All populations of S. fraterculus except Siberut, S. tahan, and S. altitudinis have a venter fur coloration homogeneously admixed with gray. The only other medium sized squirrel found in syntopy, S. tenuis, is also usually ventrally darker (admixed with gray) and dorsally lighter, with reddish-brown coloration on the shoulders and hips, white/pallid yellow hair tips present on tail, and a relatively thinner and longer tail (85–95 % of head-body length; Corbet and Hill, 1992) than S. robinsoni (56–84% of head-body length). Males of S. fraterculus, S. tahan and S. tenuis have a darker orange wash in the scrotal area than S. robinsoni, which is peach colored. Bonhote (1903) pointed out S. robinsoni has a “paler tinge” than S. lowii and S. natunensis, who have a more reddish brown agouti fur. The taxa also differ in variation of fur length on the tail, and body size and proportions. Moderate variation of hair length along the tail is observed in S. robinsoni, S. fraterculus, S. natunensis (and S. lowii Northwest Borneo populations) and S. tenuis. Extensive variation with short hair close to the base and longer hair at the tips giving a tapered appearance is observed in S. lowii and S. tahan. A tapered tail, but thickest at the base is found in S. fraterculus (Supplementary Figure S4). S. robinsoni is smaller (head-body range: 112–140 mm; body mass range: 49–77 g) than S. lowii (head-body range: 127–169 mm; body mass range: 60–125 g), somewhat smaller but overlapping in size with S. natunensis (head-body range: 133–152 mm) and slightly larger but quite overlapping in size with S. fraterculus (head-body range: 109–125 mm). Although the absolute tail length of S. robinsoni is about the same as S. lowii, S. fraterculus, and S. natunensis, it is relatively longer (56–85% of head-body length as opposed to 53–76%, 62 – 72%, and 57 – 67%, respectively). Within S. robinsoni, S. r. balae has a relatively shorter tail than the other two subspecies. S. robinsoni also has a relatively larger ear length (9–15% of head-body length) than S. lowii (7–13% of head-body length) and absolutely and relatively larger ear length than S. fraterculus (9% of head-body length) and S. natunensis (8% of head-body length) (Supplementary Table S4).

As with overall body size, the skull of S. robinsoni is smaller (33.60–36.43 mm) than S. lowii (36.93–41.03 mm); somewhat smaller but overlapping in size with S. natunensis (35.35–36.55 mm), and somewhat larger but overlapping in size with S. fraterculus (32.46–34.51 mm) (Supplementary Table S4). The skull of all of the small/medium-sized Sundasciurus are quite similar and conserved in proportion. Not many discrete characters have been found to distinguish the different species. Bonhote’s statement “the bony palate extends well back the last molar which is not the case with S. lowii” can be considered accurate for S. robinsoni, and be extended to S. fraterculus and S. natunensis as well. However, the bony palate does not always end in line with the last molars in S. lowii, and does sometimes end after, as in the other species, but quite closer. Another feature described by Bonhote, “the bullae are more flattened and rounded in S. robinsoni than S. lowii and do not project so far downwards” does not apply when additional specimens are examined. Incisor size and projection, however, do seem to be quite diagnostic. Bonhote (1903) pointed out “the molar series of S. robinsoni is much shorter and smaller than S. lowii, but incisors are about the same size.” In line with this statement, the height of the molars of S. robinsoni is clearly shorter than that of S. lowii and S. natunensis. While the incisors of S. robinsoni are less orthodont than S. tenuis, they are more orthodont than S. lowii, S. fraterculus, S. natunensis and S. tahan, which are conspicuously proodont. Finally, S. robinsoni has chisel-shaped incisors with sharp-angled beveled edges as in S. fraterculus, S. natunensis, S. lowii lingungensis and some specimens from west Borneo, but in contrast to S. lowii lowii from the remainder of Borneo (Figure 6). S. robinsoni was characterized by the PCA analysis as having the most gracile skull, with a relatively short interorbital breadth, anterior nasal breadth and least frontal bone breadth (Figure 5). It also has a shorter molar row (5.57–6.16 mm), palatal length (16.03–17.83 mm) and breadth (6.51–7.5 mm) than S. lowii (6.63–7.81 mm, 18.15–20.65 mm, 7.99–8.89 mm), but a relatively larger braincase.

Sundasciurus robinsoni has a shorter muzzle than S. natunensis, so it can be differentiated based on a shorter rostrum length (13.13–14.86/14.84–15.36 mm), palatal length (16.03–17.83/17.48–18.44 mm), and relative and absolute narrower bony palate (6.51–7.5/7.81–8.07 mm). S. robinsoni has a relative and absolute shorter molar row (5.57–6.16 mm) and smaller bullae (6.17–7.34 mm) than S. natunensis (6.72–7.15 mm; 7.26–7.65 mm) but a relatively larger proportion of bony palate of the total palatal length, relatively larger height of rostrum and relatively larger braincase (Supplementary Figure S5). Postorbital processes are relatively shorter and less pointed in S. robinsoni than S. natunensis.

There are major distinctions between almost all populations of S. robinsoni and S. fraterculus analyzed here. However, the intermediate position of the incertae sedis Kerinci population, that we have tentatively considered a subspecies (Sundasciurus robinsoni vanakeni) until molecular lines of evidence are available, generates an overlap among S. robinsoni and S. fraterculus for most diagnostic characters. The nominal S. r. robinsoni and S. r. balae have larger tympanic bullae (6.46–7.34 mm) than S. fraterculus (5.88–6.43 mm) but there is overlap with the inclusion of S. r. vanakeni (6.17–7.34 mm) in the comparison. The larger maxillary cheek teeth of S. robinsoni (5.83–6.16 mm) than S. fraterculus (5.32–5.79 mm) becomes another invalid character once S. r. vanakeni (5.6–6.16 mm) is included in the comparison. Finally, the shorter anterior nasal breadth (4.27–4.88 mm) of S. robinsoni as compared to S. fraterculus (4.79–5.31 mm) also loses its diagnostic power with the inclusion of S. r. vanakeni (4.27–5.06 mm).

Although there are clear distinctions among these two species there are no discrete craniodental measurements to distinguish them. However, the first and third components of the PCA effectively separated these, especially the latter. PC3 loadings confirm that despite the previously mentioned overlapping measurements (anterior nasal breadth, maxillary cheek teeth and tympanic bullae length), their combination can still efficiently differentiate these taxa. Mastoid breadth, diastema length and interorbital breadth also seem to be diagnostic to a lesser degree with moderate contribution in the PC3 loadings (Supplementary File S1). S. robinsoni seems to have a somewhat broader mastoid breadth and narrower interorbital breadth than S. fraterculus. S. robinsoni also exhibits a relatively narrower least frontal bone breadth than S. fraterculus. Finally, S. robinsoni has a relatively shorter muzzle than S. fraterculus, which is reflected in a relatively shorter length and height of rostrum (Supplementary Figure S5).

Distribution: Malay Peninsula, Sumatra, and Batu and Rhio Archipelagos. The subspecies S. r. robinsoni and S. r. balae have been recorded from 0–610 m, on the Malay Peninsula and Rhio and Sumatra and Batu, respectively, while S. r. vanakeni has been recorded between 900–1372 m on Mount Kerinci.

Sundasciurus natunensis (Thomas, 1895)- Natuna squirrel.

Sciurus lowii natunensis Thomas, 1895:26.—Type locality. “Sirhassen Island” (Serasan), Natuna Islands.

Holotype. BMNH 94.9.28.40, skull, skin, adult male, collected 23 September 1893 by A. H. Everett.

Emended comparison with other Sundasciurus species.

Comparisons of S. natunensis with S. robinsoni are detailed in the previous section.

The slight brown pinkish flush described by Thomas (1895) was only found to be present in the type series and is likely an effect of the preservation of these specimens. As Thomas reasoned, these were “skinned out of spirit” and certain yellows turned into red in other specimens of that same expedition. Therefore, this species exhibits the same coloration as S. lowii, and can be distinguished from S. fraterculus based on paler ventral coloration. A major external distinction is that S. natunensis (together with west Borneo S. l. lowii) lack that tapered appearance in its tail (described in the previous section) that characterizes S. lowii (from remaining Borneo and northern Natunas), with a thinner tail base, and S. fraterculus, with a thicker tail base (Supplementary Figure 4). Thomas (1895) described S. natunensis, in relation to S. lowii, as having “smaller size, shorter feet and longer ears, which in the typical variety is a mere low rim, while the Natuna subspecies has a distinct upstanding conch”. The data presented here show smaller but largely overlapping head-body size (133.3–152. mm) and shorter but overlapping feet (33–35 mm with claw) than S. lowii (127–169 mm; 32–40 mm with claw). However, the only ear measurement available for S. natunensis (provided in the original description) is smaller (12 mm) than S. lowii (12–17 mm), not larger. Nevertheless, in relation to this, perhaps the most distinctive characteristic of S. natunensis is its relatively larger bullae than S. lowii. Besides this, the skulls of S. natunensis and S. lowii are quite similar and conserved in proportion, and only differ on the relatively thinner and less deep rostrum and narrower palate of S. natunensis (Figure 6missing and Supplementary Figure S5). The skull of S. natunensis (33.60–36.43 mm) is smaller than that of S. lowii (36.93–41.03 mm), as Thomas (1895) pointed out. The other non-overlapping characters are braincase height (12.34–13.01 mm) and orbit length (8.33–9.05 mm), which are also larger in S. lowii (13.07–14.30 mm; 9.03 – 10.06 mm). Our expanded dataset also supports Thomas’ statement that S. natunensis “postorbital processes are longer and slenderer” than those of S. lowii (except West Borneo populations), but also than S. fraterculus and S. robinsoni. The presence of chisel-shaped incisors with sharp-angled beveled edges are an additional diagnostic feature to differentiate S. natunensis from S. lowii.

The species S. natunensis has a larger skull (35.35–36.55 mm), rostrum (14.84–15.36 mm), mastoid breadth (15.66 – 16.56 mm), diastema length (8.91–9.68 mm), bony palate breadth (7.81–8.03 mm), upper molar row length (6.72–7.15 mm) and bulla length (7.26–7.65) than S. fraterculus (32.46–34.51 mm; 13.550 – 14.83 mm; 14.48–15.43 mm; 7.84–8.67 mm; 6.40–6.99 mm; 5.32–5.79 mm; 5.88–6.43 mm). In addition, these taxa also differ in skull proportions, with S. natunensis having a relatively more gracile skull than S. fraterculus, which exhibits a broader skull (IB and LBC) and thicker rostrum (HR) relative to skull length.

Distribution: Southern Natuna islands (only recorded in Sirhassen Island), and possibly in west Borneo. Despite discrete characters such as postorbital processes and tail shape of the latter populations closely resemble S. natunensis, we currently consider the taxonomic status of west Borneo populations as incertae sedis given the intermediate phenotypic position of these among S. natunensis and S. lowii shown in the PCA and DAPC and the lack of genetic evidence. Genetic data from these, southern Sumatra and northern Natuna populations is needed to clarify the taxonomy and distribution of this group.

Sundasciurus lowii (Thomas, 1895) -Low’s squirrel.

Sciurus lowii Thomas, 1895: 253.—Type locality. “Lumbidan, on the mainland opposite Labuan”, Sarawak, Borneo.

Holotype. BMNH 76. 5. 2. 14., skull, skin, adult male, collected by H. Low (unknown date).

Sciurus lingungensis Miller, 1901:123.—Type locality. Lingung Island, Near Bunguran Island, Natuna Islands.

Subspecific Taxonomy

Sundasciurus lowii lowii (Thomas, 1895)

Sundasciurus lowii bangueyae (Thomas, 1908)

Sciurus lowii bangueyae Thomas, 1908:387.— Type locality. “Banguey Island” (Banggi)

Sundasciurus lowii lingungensis (Miller, 1901)

Sciurus lingungensis Miller, 1901:123.—Type locality. Lingung Island, Near Bunguran Island, Natuna Islands

Distribution: S. l. lowii is present in Borneo, S. l. bangueyae in Banggi and Balambangan islands, (and possibly Malawali island as well) and S. l. lingungensis in the Northern Natunas (recorded in Bunguran, Lingung and Laut islands).

Taxonomic Notes

Corbet and Hill (1992) pointed out that geographical variation in Sundasciurus lowii sensu lato is slight. Many of the currently recognized subspecies are not morphologically distinct from others. We consider Sundasciurus robinsoni balae and Sundasciurus robinsoni vanakeni valid subspecies based on their ventral coloration and craniodental differentiation (Figure 5 and Supplementary Figure S3). We synonymize former Sundasciurus lowii humilis with Sundasciurus robinsoni balae due to external resemblance (ventral grayish coloration in limbs). We also synonymize based on external resemblance (lack of ventral grayish coloration in limbs) Sundasciurus lowii seimundi with Sundasciurus robinsoni robinsoni. Finally, we consider Sundasciurus lowii lingungensis and Sundasciurus lowii bangueyae valid subspecies given their craniodental differentiation (Figure 5 and Supplementary Figure S3).

Molecular species delimitation analyses suggest the presence of three species in Borneo within S. lowii (in Sabah, Sarawak and East Kalimantan). However, phenotypic divergence seems to be slight and the fine scale distribution of the morphs is unknown. The only Sarawak sample included in the PCA was clustered separately but close to Sabah + east Kalimantan, that largely overlapped. Regarding the DAPC, only east Kalimantan samples were differentiated from the remaining overlapping populations. Relative size of interorbital breadth seems to differentiate a small number of Sarawak specimens from the remaining S. lowii, but these samples are clustered very close to the other S. lowii. Sarawak specimens also seem to have darker fur than Sabah and East Kalimantan S. lowii and a lack of blond tips on the tail, which are only observed in Sabah populations. We consider these populations unconfirmed candidate species (Padial et al., 2010) until a better molecular and morphological sampling is performed in terms of geographic coverage and number of specimens included.

As pointed out by Hawkins et al. (2016a) the two subgenera, Aletesciurus and Sundasciurus, proposed by Moore (1958) are not supported as reciprocally monophyletic groups. The discrete characters (presence of sagittal crest, skull size, and shape of anterior-mesial lobe of auditory bullae) described as distinguishing subgenera were not valid for all species. For instance, the Palawan mountain squirrel (Sundasciurus rabori) belongs to Aletesciurus but lacks a sagittal crest and has an intermediate antero-mesial lobe and skull size (40.9-43.6 mm) among both subgenera (Heaney, 1979; Hawkins et al., 2016a). The northern Palawan tree squirrel (Sundasciurus juvencus) shows an antero-mesial lobe that resembles that of the subgenus Sundasciurus despite being assigned to Aletesciurus. Finally, different species of the subgenus Sundasciurus such as S. lowii, S. tahan or S. altitudinis have skulls that reach sizes (41–43 mm) that overlap with those of S. rabori (authors unpublished data). Although these discrete characters do not support the subgenera Aletesciurus and Sundasciurus, Heaney (1979) found support for them in his morphometric analyses. Molecular phylogenies suggest that there are 5–6 major phylogenetic lineages within Sundasciurus, which correspond to monophyly in the subgenus Aletesciurus (which contains 3 divergent clades), and paraphyly in the subgenus Sundasciurus (which also contains 2–3 divergent clades) (Den Tex et al., 2010; Hawkins et al., 2016a). The lack of universally valid diagnostic features and the phylogenetic evidence from this and previous studies demonstrate that current Sundasciurus subgeneric classification is invalid, so we synonymize Aletesciurus with Sundasciurus.

Discussion

Deep genetic and moderate phenotypic divergences between the lowland Sunda tree squirrels formerly recognized under the species Sundasciurus lowii were confirmed. The consistency between these two lines of evidence along with geography required the revalidation of Sundasciurus robinsoni and species level recognition of S. natunensis.

The patterns of divergence suggest ancient diversification events between these taxa that predate the Pleistocene. This does not support the Pleistocene speciation-pump hypothesis as the main driver behind the high levels of diversity found in Sundaland (Brown et al., 2013; Husson et al., 2019). The ancient patterns of divergence that we found in Sunda squirrels are consistent with those found in other small vertebrates and invertebrates (Klaus et al., 2013; Grismer et al., 2017; Karin et al., 2017; Manawatthana et al., 2017). This is particularly interesting in light of recent geological work suggesting that the Sunda shelf was continuously exposed until 400 ka (Husson et al., 2019). This suggests that dispersal should have been fluid across the entire shelf until this point. The much older times of divergence suggest that there was substantial structure, and likely speciation prior to the inundation of the Sunda shelf in the Middle Pleistocene.

The dated phylogeny also suggests that Sunda squirrels experienced a rapid radiation, given the almost simultaneous diversification event during the middle-late Miocene transition among major lineages of Sundasciurus. The mito-nuclear discordant relationship among these major lineages supports this hypothesis. There are two potential processes that might have driven this gene discordance, incomplete lineage sorting and hybridization, both of which have been viewed as byproducts of rapid radiations (Degnan and Rosenberg, 2006; Wiens et al., 2006; McLean et al., 2016). The latter has proved to be a widespread process across Sciuridae (Garroway et al., 2010; Chang et al., 2011; Chavez et al., 2014; McLean et al., 2016; Hawkins et al., 2016b), and may be the result of speciation outpacing the evolution of reproductive incompatibilities (Wiens et al., 2006). Interestingly, similar diversification patterns have also been shown during this narrow timeframe in stream toads, halfbeaks and colugos (De Bruyn et al., 2013; Grismer et al., 2017; Mason et al., 2019).

Land bridges that were formed after the inundation of the Sunda shelf may have allowed disparate dispersal rates between the different land masses. Whether or not a species successfully dispersed during these periods was likely due to a number of factors including the available habitat on the exposed shelf, distance between land masses, and presence or absence of close relatives/competitors for the same resources. Previous studies have suggested that during this time, despite shelf exposure, many rainforest taxa remained isolated (Den Tex et al., 2010; Leonard et al., 2015; Mason et al., 2019).

The Mentawai Islands are geographically very close to Sumatra, but were not connected by land bridges to Sumatra when other land masses were connected due to a deep-sea trench between them (Voris, 2000; Roberts et al., 2011). The endemic fraternal squirrel became isolated from Robinson’s squirrel, which was distributed across Sumatra and the Malay Peninsula in the Pliocene. This level of deep divergence between taxa on Mentawai and sister taxa elsewhere in Sundaland has also been identified in other taxa such as gibbons, civets and skinks (Chan et al., 2010; Patou et al., 2010; Karin et al., 2017) and even deeper divergence was found in macaques (Evans et al., 2017). This sister relationship between Mentawai and Sumatra-Malay Peninsula populations (Meijaard and Groves, 2004; Chan et al., 2010; Patou et al., 2010; Arifin et al., 2018) however, is not the most common in mammals (Wilting et al., 2012). For example, the Mentawai treeshrew (Tupaia chrysogaster) is sister to the endemic Bornean Long-footed treeshrew (T. longipes, Roberts et al., 2011) and the Mentawai endemic Koopman’s pencil-tailed tree mouse (Chiropodomys karlkoopmani) is sister to the Bornean endemic greater pencil-tailed tree mouse (C. major/C. calamianensis) (Musser and Carleton, 2005).

Although bathymetry suggests that the Mentawai islands have not been directly connected to Sumatra, it is possible that they have been connected to the Batu Islands which have been connected to Sumatra, creating a peninsula parallel to the latter (Voris, 2000; Roberts et al., 2011). Potential episodes of gene flow across this past connection in both directions could explain the lack of differentiation in nuclear introns but high mitochondrial divergence and important phenotypic differentiation among populations. The intermediate phenotype of the Siberut population, the closest Mentawai island to the Batu archipelago, which externally resembles Robinson’s squirrel but is clustered within the craniodental variation of the Fraternal squirrel would also support this. In accordance with this peninsular dispersal corridor, Karin et al. (2017) suggest a divergence of ∼2.60 Mya among skink populations of Nias (Batu Islands) and Siberut (Mentawai Islands) which would be slightly after to the divergence among Mentawai and Malay Peninsula squirrels. Given that western Sumatra and Mentawai have been suggested to be a forest refuge (Morley, 2000; Lim et al., 2010) it is not straight-forward to assess the direction of the potential dispersal events. Intron haplotype networks suggest that S. fraterculus is within the diversity of S. robinsoni, which could suggest a Malay Peninsula to Sumatra to Batu to Mentawai colonization. However, within Mentawai, just South Pagi was sampled, so data from other populations from Mentawai, Sumatra and Batu islands are needed to confirm this.

There is a deep, Late Miocene divergence between the western Robinson’s squirrel and Fraternal squirrel and the eastern Low’s squirrel and Natuna squirrel. At about the same time the western high elevation Sunda squirrels (S. tahan and S. altitudinis) diverged from the western slender and Jentink’s squirrels. This period was also characterized by diversification across Sundaland in semiarboreal rough-skinned skinks, bulbuls, barbets and Sundaic freshwater crabs (Den Tex et al., 2010; Klaus et al., 2013; Karin et al., 2017; Manawatthana et al., 2017). These common patterns could be due to dispersal barriers following the development of the Bornean highlands and river systems into their more modern representation by ∼8–10 Mya (Hall, 2013; Mason et al., 2011; Yang and Rannala, 2010). Alternatively, vicariance during fragmentation of the Sunda shelf/oceanic dispersal between islands during periods of high sea-levels by ∼5–6 Mya (Karin et al., 2017) if there were early inundations, could also explain this pattern.

The Natuna Islands are located in the South China Sea between Borneo and mainland Asia, a key location for understanding the biogeography of the region. Few studies to date (Den Tex et al., 2010; Mason et al., 2011, 2019) have been able to gather any information from historic specimen samples from this archipelago possibly due to poor specimen and DNA preservation. A deep divergence between populations of colugos from the northern and southern Natuna islands has been identified, with the former populations closely related to northeastern Borneo populations and the latter populations to the rest of Borneo (Mason et al., 2019). Here we show that the Natuna squirrels from the southern islands are deeply divergent from Low’s squirrels of central and east Borneo (Figure 3). These squirrels from the southern Natuna islands are also morphologically distinct from the northern Natuna squirrels, which are more similar to central and east Borneo specimens (Figure 5). Finally, populations from west Borneo have an intermediate phenotype (Supplementary Figure S3). These squirrels might have a similar pattern to the colugos, with one species on the southern Natuna islands and possibly extreme west Borneo (the Natuna tree squirrels, S. natunensis) and the other on the northern islands and rest of Borneo (S. lowii). However, additional genetic and morphologic data from the Natuna islands and Borneo is needed in order to confirm this hypothesis and the taxonomic status of west Borneo and north Natuna populations.

While divergence in allopatry seems to be the main process driving the speciation of Low’s squirrels, ecomorphological and behavioral adaptations in this clade suggest an important role of post-vicariance/colonization niche divergence, and a lack of niche conservatism. Previous studies have shown a significant correlation between locomotion modes and morphological traits in squirrels and other rodents. Tail length is thought to be associated with arboreality in squirrels (Hayssen, 2008) and other small mammals such as deer mice (Horner, 1954; Kingsley et al., 2017) and treeshrews (Martin, 1968). Nations et al. (2019) found a positive relationship of tail length with climbing in Philippine murines, but also a large variance and low predictive power for this trait across all locomotor categories. Thorington and Heaney (1981) provided a comparative table with different ratios that also shows this trend (e.g., the most terrestrial Dremomys rufigenis and S. lowii exhibit the relative shortest tail among Callosciurinae). Arboreal Sigmondontines and tree squirrels are known to have a more rounded and inflated vault and shorter nasals than terrestrial species (Lu et al., 2014; Camargo et al., 2019). Arboreal squirrels and Procyonids also seem to have shorter nasals or slightly shorter olfactory bulbs compared to terrestrial ones (Lu et al., 2014; Ahrens, 2014; Bertrand et al., 2019). A longer rostrum could imply a better sense of smell, which might be useful when living on the ground (Lu et al., 2014) while brain enlargement is an adaptation in arboreal species that use the space in three dimensions and need to integrate information of complex environments (Camargo et al., 2019). The relatively long tail, larger braincase and shorter rostrum of S. robinsoni as compared to the terrestrial S. lowii and S. natunensis might suggest arboreality. Ecological field studies (Kemper and Bell, 1985; Sharma, 1992; Abdullah et al., 2001; Yasuda et al., 2003; Saiful and Nordin, 2004; Jayaraj et al., 2013) and diet data finding bark, fruit and insects (Harrison, 1962; Abdullah et al., 2001) are also consistent with an arboreal life style for S. robinsoni. On the other hand, S. lowii on Borneo seems to be much more terrestrial (Davis, 1962; Wells et al., 2004; Camacho-Sánchez et al., 2019) with a diet including small arthropods, fruits and fungus (Davis, 1962). The Fraternal squirrel seems to have intermediate habits, since it uses the ground extensively, but also climbs showing a preference for heights up to five meters, but avoids the upper canopy (Whitten, 1981). It has only been observed eating bark, moss and lichen, but stomach content analyses also found annelid remains, in accordance with more terrestrial behavior.

The relatively longer rostrum of S. lowii, S. fraterculus and S. natunensis with respect to S. robinsoni might be due to its terrestrial foraging behavior, reflecting a more developed olfactory bulb as an adaptation to a more insectivorous diet. In addition, the more robust skull of S. fraterculus with regard to S. robinsoni might reflect a possible adaptation to the harder dietary elements such as an increased consumption of chitin from insects, lichens and bark, and an absence of soft food items such as fruit (Whitten, 1981). These findings suggest that not only has there been a high degree of genetic and phenotypic divergence, but also that morphological changes have possibly been the result of natural selection driven by adaptation to different niches, and are not only due to drift caused by geographic isolation (Mayr, 1963). It is difficult to interpret if the different degrees of scansoriality found among these species are due to different levels of interspecific competition, predator avoidance, food limitation or a combination of some of these. In any case, differences in the squirrel and treeshrew communities among islands are notable. On the Mentawai islands, S. fraterculus has to share its habitat with just one diurnal treeshrew (Tupaia chrysogaster), one terrestrial squirrel (Lariscus obscurus) and one arboreal diurnal squirrel (Callosciurus melanogaster), allowing it a greater range of vertical habitat than its sister species (Whitten, 1981). On Borneo, S. lowii has to share its habitat with up to six diurnal treeshrews and seven species of diurnal squirrels. Two of the treeshrews (Tupaia minor and Tupaia gracilis) and most of the squirrels are arboreal, and three of the latter are bark gouging specialists (Sundasciurus tenuis, Nannosciurus melanotis and Exilisciurus exilis) (Emmons, 2000; Zelditch et al., 2017). However, the only two terrestrial squirrels (Lariscus insignis or L. hosei and Rhinosciurus laticaudatus) are in very low densities or possibly absent where Low’s squirrel is present (Wells et al., 2007; Yasuda and Tsuyuki, 2012; Fitzmaurice, 2014; Chapman et al., 2018; Camacho-Sánchez et al., 2019; authors unpublished data) and are specialized on fruits and insects. On Sumatra and the Malay peninsula S. robinsoni is co-distributed with diurnal arboreal and terrestrial treeshrews (Tupaia minor and T. glis/T. ferruginea), and up to ten diurnal squirrels, three of which are terrestrial (Lariscus insignis, Menetes berdmorei and Rhinosciurus laticaudatus) (Lekagul and McNeely, 1977). These terrestrial squirrels seem to be more common in mainland field-studies than in Borneo (Kemper and Bell, 1985; Yasuda et al., 2000; Syakirah et al., 2001; Yasuda et al., 2003; Jayaraj et al., 2013). On the other hand, the bark gouging niche seems to be less exploited in Peninsular Malaysia, with just S. tenuis (Zelditch et al., 2017). This could have potentially added evolutionary selective pressures on the ancestor of S. robinsoni to alter its diet to avoid competition with other species. The relatively larger and more orthodont incisors of this species with regard to its closest relatives are also in line with this hypothesis. Convergence towards a bark gouging diet seems to have arisen at least seven times in Sciuridae and it exhibits higher rates of convergence relative to the number of living species (10 species) within that dietary niche (Zelditch et al., 2017).

Within Borneo, deep genetic divergences were detected across Sundasciurus lowii. Bornean populations exhibit a high level of genetic and some morphological diversity, but the fine scale distribution of this diversity is not yet characterized. Molecular species delimitation analyses suggest candidate species found in Sabah, Central Sarawak and East Kalimantan. However, much more sampling will be required to clarify their taxonomy. As in previous studies on Sunda squirrels, our results also point towards a possible underestimation of highland species diversity given the high genetic divergence observed among populations of S. altitudinis and high phenotypic divergence found between Mount Kerinci mid-elevation (900–1400 m) populations (S. r. vanakeni) and the remaining Sumatra and Malay Peninsula lowland populations (S. r. robinsoni and S. r. balae). The lack of molecular evidence for this pelage and craniodentally distinct population suggests to consider S. r. vanakeni as an unconfirmed candidate species (Padial et al., 2010). Given that a previous attempt failed to amplify DNA from the only known series of specimens of this taxon (Den Tex et al., 2010) future collection of new specimens along an elevational gradient in Kerinci will be needed to confirm the specific status of this population.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by Smithsonian Institution, National Museum of Natural History, proposal number 2012-04 and Estación Biológica de Doñana proposal number CGL2010-21424.

Author Contributions

AH, MH, JM, and JL designed the study. AH, MH, and AA performed field and museum data collection. AH and MH performed laboratory work and statistical analyses. AH and JL wrote the manuscript with input from MH, AA, and JM. All authors approved the submitted version.

Funding

AH was supported by an Ernst Mayr Travel grant and AMNH Collection Study Grant during museum data collection. This research received support from the SYNTHESYS Project financed by European Community Research Infrastructure Action under the FP7 “Capacities” Program. AH was also supported by a Spanish Ministry of Economy contract CGL2014-58793-P. The Spanish Ministry of Science and Innovation grant CGL2017-86068-P to JL supported this study. MH was supported by discretionary research funding from Humboldt State University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

L. Heaney provided important advice and guidance during specimen examination and morphometric data collection. G. Forcina provided valuable insights during marker development. L. Flynn and J. Barry provided expertise and guidance regarding fossil calibrations and A. Jaramillo and A. Sánchez regarding taxonomical nomenclature. M. Camacho, P. Carro, M. Lopez and different field assistants participated in fieldwork. Sabah Parks, Sabah Forestry Department, and Sabah Biodiversity Council provided research permits and logistic support, especially FTY Yuh. R. Portela, L. Heaney, D. Lunde, R. Voss, C. Urdiales, and M. Omura provided assistance at NHMUK, FMNH, USNM, AMNH, EBD-CSIC, and Harvard MCZ collections. L. Heaney, J. Esselstyn, A. Ferguson, L. Smith, B. Lim, and D. Lunde granted access to several tissue samples hosted at FMNH, ROM, and USNM. Nancy McInerney for help with MiSeq runs at the CCG. B. Lim, O. Johnson, M. La, R. J. Backhouse, and C. Lee shared pictures. Logistical support was provided by Laboratorio de Ecología Molecular, Estación Biológica de Doñana, CSIC (LEM-EBD) and Doñana ICTS-RBD.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00208/full#supplementary-material

References