\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{High gene flow across heterogeneous tropical montane environments in a
Bornean endemic small mammal}



\author[1]{Lillian Parker}%
\author[1]{Melissa Hawkins}%
\author[2]{Miguel Camacho-Sanchez}%
\author[1]{Michael Campana}%
\author[1]{Jacob West-Roberts}%
\author[1]{Tammy Wilbert}%
\author[3]{Haw Chuan Lim}%
\author[3]{Larry Rockwood}%
\author[2]{Jennifer Leonard}%
\author[1]{Jesús Maldonado}%
\affil[1]{Smithsonian Conservation Biology Institute}%
\affil[2]{Estación Biológica de Doñana}%
\affil[3]{George Mason University College of Science}%


\vspace{-1em}



  
  \date{June 08, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Rapid shifts in environmental variables associated with elevational
changes in montane ecosystems provide opportunities to test hypotheses
regarding the effects of environmental heterogeneity on gene flow and
genetic structure. In tropical mountains, spatial environmental
heterogeneity combined with seasonal environmental stability is
predicted to result in low dispersal across elevations. Few studies have
investigated the genetic consequences of elevational environmental
heterogeneity in tropical montane mammals. Here, we use a population
genomics approach to test the hypothesis that mountain treeshrews
(\emph{Tupaia montana}) exhibit limited gene flow across elevational
gradients and between two neighboring peaks within Kinabalu National
Park (KNP) in Borneo. We sampled 83 individuals across elevations on Mt.
Tambuyukon (MT) and Mt. Kinabalu (MK) and sequenced mitogenomes and
4,106 ultraconserved elements containing an average of 1.9 single
nucleotide polymorphisms per locus. We detected high gene flow across
elevations and between peaks. We found greater genetic differentiation
on MT than MK despite its lower elevation and associated environmental
variation. This implies that, contrary to our hypothesis, genetic
structure in this system is not primarily shaped by elevation. We
propose that this pattern may instead be the result of colonization
history combined with restricted upslope gene flow on MT due to unique
plant communities associated with its upper montane habitats. Our
results serve as a foundation to identify and mitigate future effects of
climate change on mountain treeshrews in KNP. Given predictions for 2100
CE, we predict that mountain treeshrews will maintain genetic
connectivity in KNP, making it an important conservation stronghold.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Montane ecosystems are highly valuable to the study of evolutionary and
genetic consequences of environmental heterogeneity due to the rapid
shifts in environmental variables (e.g. temperature, precipitation, and
solar radiation) over short distances. Relative to those in the
temperate zone, tropical mountains experience less seasonal temperature
variation which generates greater elevational stratification (Polato et
al., 2018). Janzen (1967) hypothesized that this temporal thermal
stability paired with spatial environmental heterogeneity should select
for narrow thermal tolerances, which in turn result in low effective
dispersal and population isolation across elevational gradients
(Ghalambor, Huey, Martin, Tewksbury, \& Wang, 2006; Gill et al., 2016).
Consistent with Janzen's (1967) hypothesis, many studies have documented
narrower elevational ranges among tropical montane species than
temperate species across diverse taxonomic groups (McCain, 2009;
Ghalambor et al., 2006). Fewer studies have tested the prediction
stemming from Janzen's (1967) hypothesis that restricted gene flow among
populations spanning elevational gradients results in genetic
divergence. The available data regarding this prediction are
contradictory - some studies have found significant population genetic
divergence across elevations, for example, in insects (Polato et al.,
2018; Gueuning et al., 2017) and birds (Gadek et al., 2018; DuBay \&
Witt, 2014; Linck, Freeman, \& Dumbacher, 2019). Others have found high
rates of gene flow, sometimes in the presence of adaptive phenotypic
divergence (Cheviron \& Brumfield, 2009; Gadek et al., 2018; Branch,
Jahner, Kozlovsky, Parchman, \& Pravosudov, 2017).

Very few studies have investigated the population genetic structure of
small mammals across elevational gradients in tropical montane
ecosystems (Muenchow, Dieker, Kluge, Kessler \& von Wehrden, 2018; but
see Yu, 1995). We address this knowledge gap by elucidating the
population genetic structure of the mountain treeshrew, \emph{Tupaia
montana} , across its full elevational range on two mountains in
Kinabalu National Park (KNP), Sabah, Borneo: Mt. Kinabalu and Mt.
Tambuyukon (Figure 1). The mountain treeshrew provides an interesting
system in which to study the effect of environmental gradients on
population structure because it has a broad elevational distribution
compared to other small mammals in KNP (Camacho-Sanchez, Hawkins, Tuh
Yit Yu, Maldonado, \& Leonard, 2019; Nor, 2001).

Understanding the population genetic structure of tropical montane taxa
like the mountain treeshrew is important for conservation because it
enables researchers to identify metapopulation dynamics and distinct
evolutionary units warranting protection and to predict and track
species' responses to changing environmental conditions (Moritz, 1994;
Castillo Vardaro, Epps, Frable, \& Ray, 2018; Camacho-Sanchez et al.
2018). This is critical given the vulnerability of tropical montane
ecosystems to the impacts of global climate change (GCC) (Feeley,
Stroud, \& Perez, 2017; Lenoir \& Svenning, 2015) and the paucity of
population genetic studies in Southeast Asia (Muenchow et al., 2018).

\emph{Kinabalu National Park}

Kinabalu National Park (6\selectlanguage{ngerman}°09\selectlanguage{english}'N 116deg39\selectlanguage{english}'E) is a UNESCO World Heritage Site
located in Sabah, Malaysia. The park comprises two major mountains, Mt
Kinabalu (MK) and Mt. Tambuyukon (MT), which contain a wide range of
habitats from lowland tropical forest to high elevation scrub (Kitayama,
1992). At 4,095 meters above sea level (masl), MK is the tallest
mountain in the Sundaland biogeographical region. It is relatively
young, having reached its present height ca. 1 million years ago (Mya)
(Hall et al. 2009). Eighteen kilometers to the north of MK, the
less-studied MT stands at 2,579 masl (Figure 1a). MT is older - its
major uplift occurred 7--11 Mya as part of the Crocker Range (Hall et
al., 2009). Although multiple studies have focused on species richness
across elevational gradients in KNP (e.g. Camacho-Sanchez et al., 2019;
Liew, Schilthuizen, \& bin Lakim, 2010; Nor, 2001), studies on
population genetic structure in this landscape are limited
(Camacho-Sanchez et al., 2018).

\emph{Mountain treeshrews}

Treeshrews are small-bodied mammals in the order Scandentia native to
South and Southeast Asia. The most comprehensive subordinal molecular
phylogeny of treeshrews to date showed that Borneo is a center of
treeshrew diversity, with 9 of 20 extant species including both
widespread and endemic lineages (Roberts, Lanier, Sargis, \& Olson,
2011).

The mountain treeshrew is the only montane \emph{Tupaia} species endemic
to Borneo. Its distribution includes montane regions of Sabah and
Sarawak (Figure 1 inset) (Payne and Phillipps 2016; Hawkins, 2019). The
species diverged from its sister lineage, the ruddy treeshrew (\emph{T.
splendidula} ), ca. 4--7 Mya (Roberts et al., 2011). The ruddy treeshrew
and other closely related species occupy lowland forest, potentially
indicating elevational niche partitioning (Han, Sheldon, \& Stuebing,
2000). On MK, the mountain treeshrew's range spans from ca. 900 to 3200
masl, encompassing four vegetation zones following Kitayama (1992):
lowland (\textless{}1200 masl), lower montane (1200--2350 masl), upper
montane (2350--2800 masl), and subalpine (2800--3400 masl). On MT, the
species spans from ca. 900 masl to the summit and three vegetation zones
(Camacho-Sanchez et al., 2019; Phillipps \& Phillipps, 2016). The
ecology, behavior, and evolution of this species are poorly known. The
most comprehensive ecological study to date was conducted by Emmons
(2000), who tracked individuals for one month at Poring Hot Spring, MK
(900 masl). Emmons (2000) suggests that the mountain treeshrew has small
home ranges compared to other \emph{Tupaia} species, at only 2.5
hectares. Despite the lack of ecological information, the mountain
treeshrew's broad elevational distribution permits hypothesis testing
regarding the effects of environmental heterogeneity on population
genetic structure.

\emph{Study objectives}

Our first objective is to determine if there is significant genetic
differentiation between mountain treeshrews on MT and MK, which would
have implications for conservation monitoring and management of the
species in KNP and beyond. Our second goal is to assess population
genetic structure across the species' elevational extent on both peaks.
We test the hypothesis that gene flow is restricted across the steep
ecological gradient that mountain treeshrews inhabit. We predict that
mountain treeshrews will exhibit significant differentiation in neutral
genetic markers 1) between mountain peaks, due to limited dispersal
across the lowland habitat that connects them, and 2) across elevations
-- with greater differentiation on MK due to its higher elevation and
associated environmental variability. Our final objective is to evaluate
the utility of ultraconserved element (UCE) loci for estimating
population genetic parameters and structure at fine geographic scales
within a species. Previous studies have shown that UCEs are sufficiently
variable to resolve phylogenies on a phylogeographic scale (e.g. Harvey,
Smith, Glenn, Faircloth, \& Brumfield, 2016; Mason, Olvera-Vital,
Lovette, \& Navarro- Siguenza, 2018; Smith, Harvey, Faircloth, Glenn, \&
Brumfield, 2014) and to answer questions regarding recently diverged
species (Oswald et al., 2016; Winker, Glenn, \& Faircloth, 2018).
However, based on a Web of Knowledge literature search (accessed April
12, 2020) using the keywords `population' AND `ultraconserved element',
ours is the first study to use these markers to study population
genetics at a landscape scale. We describe a modified UCE processing
pipeline to generate single nucleotide polymorphism (SNP) and phased
pseudo-haplotype sequence datasets.

\textbf{Materials and methods}

\emph{Sample collection}

We trapped small mammals on MK and MT during two field seasons in 2012
and 2013. Our trapping methodology and permitting information is
described in Camacho-Sanchez et al. (2019). Briefly, we set traps from
ca. 331 to 2,509 masl on MT and from ca. 503 to 3,466 masl on MK (Figure
1a). The mountain treeshrew was the most frequently caught species,
representing 37.5\% of all catches. For this study, we included
92\emph{Tupaia} individuals from our collection: 84 mountain treeshrews
and 8 outgroup individuals from 2 congeners, the pygmy treeshrew
(\emph{T. minor, n = 2} ) and the large treeshrew (\emph{T. tana, n = 6}
) (Table S1). We also received a loan of genomic DNA from two ruddy
treeshrew voucher specimens housed at the University of Michigan Museum
of Zoology (UMMZ174428 and UMMZ174429).

\emph{Laboratory methods}

All laboratory work was performed at the Center for Conservation
Genomics (CCG), Smithsonian Conservation Biology Institute, Washington,
DC. We extracted DNA from liver and ear punch samples using a DNeasy
Blood and Tissue Kit (Qiagen, Valencia CA) following the manufacturer's
protocol. We amplified whole mitogenomes in two fragments using long
range PCR, fragmented the PCR products to an average length of 500 base
pairs (bps) using a Qsonica Q800R sonicator (QSonica, Newtown, CT, USA),
and prepared single-indexed DNA libraries for sequencing using a Kapa
LTP Library Preparation kit (Kapa Biosystems, Wilmington, MA) following
Hawkins et al. (2016). We pooled libraries equimolarly and sequenced on
an Illumina MiSeq with 2 x 100 base pair (bp) reads (Illumina, Inc., San
Diego, CA).

We used in-solution DNA hybridization to enrich genomic DNA for UCEs
following Hawkins et al. (2016). We sheared DNA extracts and constructed
indexed libraries as above. We quantified libraries using a Qubit\selectlanguage{english}?
fluorometer (Life Technologies) with a 1x dsDNA HS assay kit and
multiplexed 4-8 samples equimolarly prior to enrichment. We used a
NimbleGen SeqCap EZ\selectlanguage{english}? kit (Roche, Basel, Switzerland) containing 54,689
unique 60-bp DNA probes representing 5,561 vertebrate UCE loci with an
average of 4x tiling per base per locus to enrich multiplexed libraries
following the manufacturer's protocol. Post-enrichment libraries were
amplified with 12--14 cycles of PCR using Kapa HiFi HotStart DNA
polymerase (Kapa Biosystems, Wilmington, MA) following the
manufacturer's protocol.

Following visualization on a Bioanalyzer 2100 (Agilent Technologies,
Santa Clara, CA) with High Sensitivity DNA kits, enriched libraries were
quantified via qPCR using the Kapa Biosystems Illumina Library
Quantification Kit (Kapa Biosystems, Wilmington, MA). Samples were
pooled equimolarly and sequenced with 2 x 150 bp reads on Illumina
HiSeq2000 (Semel Institute of Neurosciences, UCLA, \& University of
Copenhagen, Denmark) and MiSeq\selectlanguage{english}? platforms (CCG).

\emph{Data analysis: Mitogenomes}

\emph{Assembly and alignment}

Mitogenome amplicon reads were quality filtered with Trimmomatic v.0.33
(Bolger, Lohse, and Usadel 2014) with parameters SLIDINGWINDOW: 4:15 and
MINLEN: 36. Since the only publicly available mitogenome representing
any \emph{Tupaia} species (the northern treeshrew, \emph{T.
belangeri}NC\_002521; Schmitz, Ohme, \& Zischler, 2000) is highly
divergent from our study species (Roberts et al., 2011), we first
generated reference mitogenomes for the mountain treeshrew and 3 more
closely related outgroup species: the pygmy treeshrew, large treeshrew,
and ruddy treeshrew. For each species, we selected one individual with
the highest number of sequencing reads (pygmy treeshrew, BOR 443; large
treeeshrew, BOR 010, \& ruddy treeshrew, UMMZ174429) and assembled
sequences\emph{de novo} using the MIRA v1.0.1 plugin in Geneious v9.1.2
(Biomatters Ltd.). Quality filtered sequence reads were mapped to the
appropriate reference using BWA-MEM v0.7.10 (Li 2013) with default
parameters. We also assembled mitogenomes from UCE-enriched library
sequences (Supplemental Information). Consensus sequences were generated
with Geneious (lowest coverage to call a base 5x, and Highest Total
Quality parameters) and aligned with the MAFFT v7.450 plugin (Katoh et
al\emph{.} 2002). We transferred annotations from the northern treeshrew
reference to the consensus sequences. To rule out the presence of
nuclear copies of mitochondrial genes (NUMTs), we translated all
protein-coding genes to check for frame shifts or stop codons.

\emph{Genetic diversity and population structure}

Because the inclusion of close relatives can bias estimates of genetic
diversity and structure (Goldberg \& Waits, 2010), we removed sequences
representing first-order relatives identified by our SNP dataset
(methods described below) and performed all subsequent mitogenome
analyses with the reduced data, which we refer to as the ``unrelated
dataset''. We defined haplotypes and calculated haplotype diversity
(\emph{H} \textsubscript{d}), nucleotide diversity \((\pi\)),
and Tajima's \emph{D} using DNAsp (Librado \& Rozas 2009). We estimated
the differentiation between MK and MT and between high and low
elevations within each peak through analysis of molecular variance
(AMOVA) in Arlequin v3.5 (Excoffier \& Lischer, 2010) with a permutation
test of 10,000 replicates to assess statistical significance. We
visualized relationships among haplotypes by generating a median-joining
network with the PopART software (Leigh \& Bryant, 2015).

\emph{Phylogenetic analysis and modeling demographic history}

We performed phylogenetic analyses to place the mitochondrial lineages
detected in our mountain treeshrew samples within an evolutionary
framework with respect to other Bornean treeshrew species (i.e. the
large treeshrew, pygmy treeshrew, and ruddy treeshrew). We used
PartitionFinder v2.0 (Lanfear et al. 2012) to select partitions and
substitution models (Table S2) and estimated a phylogeny using MrBayes
v3.2.6 (Ronquist \& Huelsenbeck 2003). We then used BEAST v.1.8.4
(Drummond \& Rambaut 2007) to estimate the timing of divergence between
the mountain treeshrew mitochondrial lineages we identified.

To infer demographic history, we performed a Bayesian coalescent skyline
plot analysis using BEAST v2.0 (Bouckaert et al. 2014). We used a time
to most recent common ancestor (TMRCA) prior of 450,000 years before
present (lognormal distribution, \selectlanguage{greek}µ \selectlanguage{english}= 0.45, \selectlanguage{greek}σ \selectlanguage{english}= 0.2), the estimated date
of divergence between the two mitochondrial lineages as determined by
the dating analysis performed in BEAST. Analytical details are provided
in Supplemental Information.

\emph{SNPs within UCEs}

\emph{Genotyping SNPs}

To generate the SNP dataset, we followed the PHYLUCE v.1.5.0 pipeline
with default parameters (Faircloth, 2016) for sequence trimming,\emph{de
novo} assembly of contigs, identification of UCE loci, and sequence
alignment. We generated a pseudo-genomic reference by aligning each
locus with MAFFT v7.407 and trimming using Gblocks v0.91b with default
parameters (Castresana, 2000). We then used Geneious v9.1.2 to generate
a consensus sequence for each locus, replacing ambiguity codes with an
appropriate nucleotide at random. We used Picard v1.106
(http://broadinstitute.github.io/picard/), and SAMtools v1.9 (Li et al.,
2009) to create sequence dictionaries and reference indices from the
reference. We used the PHYLUCE script \emph{snps.py} to automate
alignment of trimmed reads from each sample to the reference with
BWA-MEM v0.7.17, and then called SNPs with the HaplotypeCaller tool of
the Genome Analysis Toolkit v3.7 (McKenna et al., 2010) following Giarla
and Esselstyn (2015). Using VCFtools v0.1.14 (Danecek et al., 2011), we
removed SNPs that failed to pass GATK quality filters (QD \textless{}
2.0 \textbar{}\textbar{} FS \textgreater{} 60.0 \textbar{}\textbar{} MQ
\textless{} 40.0 \textbar{}\textbar{} HaplotypeScore \textgreater{} 13.0
\textbar{}\textbar{} MappingQualityRankSum \textless{} -12.5
\textbar{}\textbar{} ReadPosRankSum \textless{} -8.0), and selected SNPs
with a minimum depth of coverage of 8 per individual and a minor allele
frequency [?] 5\%. We used \emph{HD\_plot.py} (McKinney, Waples, Seeb, \&
Seeb, 2017) to filter SNPs resulting from putative paralogs or wrongly
assembled contigs from the dataset by removing SNPs with heterozygosity
\textgreater{} 0.75 and a read-ratio deviation score D \textgreater{}10.
The D statistic is a measure of deviation from the expected allelic read
ratio of 1:1 when reads are summed over all heterozygous individuals.
This method more accurately identifies true SNP loci than methods
relying on read depth or heterozygote excess alone (McKinney et al.,
2017). After filtering with \emph{HD\_plot.py} , we removed SNPs that
were out of Hardy-Weinberg Equilibrium after Bonferroni correction for
multiple comparisons (\emph{p} \textless{} 10\textsuperscript{-5}) using
VCFtools v0.1.16\emph{.} To generate a set of unlinked SNPs, we selected
a single SNP per UCE using VCFtools v0.1.16 `-thin 2000'. We used the
unlinked SNP dataset with 10\% missing data for all SNP-based analyses
except calculation of genetic diversity and effective population size,
Principal Components Analysis (PCA), Discriminant Analysis of Principal
Components (DAPC) and Bayesian cluster analysis with STRUCTURE v2.3.4
(Pritchard, Stephens, \& Donnelly, 2000), for which we used a dataset
with no missing data.

\emph{Generating phased pseudo-haplotype sequences}

In order to include multiple SNPs per UCE locus as well as invariant
sites, we generated multiple sequence alignments of pseudo-haplotypes.
We did this by using the EMIT\_ALL\_SITES output mode of the GATK
HaplotypeCaller tool. We filtered the resulting VCF file to include only
UCE loci with at least one SNP with no more than 10\% missing data. We
then generated alignments from the VCF file with a custom Ruby
script,\emph{vcf2aln} v0.4.2 (https://github.com/campanam/vcf2aln,
Supplemental Information). This script utilizes phasing information
where present and randomly selects an allele where phase is unresolved.

We trimmed aligned, phased UCE sequences with Gblocks Version 0.91b
(Castresana, 2000) using default parameters and quantified informative
sites with the \emph{phyluce\_align\_get\_informative\_sites.py} PHYLUCE
script. For the final dataset, we retained loci with more than one and
fewer than 10 parsimony informative sites (PIS). We removed loci with
more than 10 PIS because these are likely the result of assembly or
alignment errors.

\emph{Genetic diversity and effective population size}

We removed individuals that were identified as first-degree relatives
(parent-offspring or full siblings) according to the KING v2.1.4
software (Manichaikul et al., 2010), i.e. those with kinship
coefficients [?] 0.18. Using the unlinked SNPs with 10\% missing data, we
first calculated pairwise kinship values and identified putative family
groups with the KING software and then ran a PC-AiR analysis with
GENESIS v2.2.2 (Conomos, Reiner, Weir \& Thornton, 2016) in R (R Core
Team 2017, applies to all subsequent use of R) to identify an
``unrelated'' subset of individuals. We used GenAlEx v6.503 to estimate
the power of the SNP dataset to differentiate individuals by
calculating\emph{P\textsubscript{IDsib}} , the probability of two
individuals having identical genotypes assuming siblings are present in
the data (Waits, Luikart, \& Taberlet, 2001).

To estimate genetic diversity, we calculated average expected and
observed heterozygosity (\emph{H\textsubscript{e}}
and\emph{H\textsubscript{o}} ) and the inbreeding coefficient
(\emph{F\textsubscript{IS}} ) with VCFtools v0.1.16 using our SNPs with
no missing data. We also concatenated FASTA alignments of UCE sequence
pseudo-haplotypes end-to-end for all individuals in the unrelated
dataset and used the maximum composite likelihood method to calculate
nucleotide diversity (\(\pi\)) in MEGA v7.0.26 (Kumar,
Stecher, \& Tamura, 2016).

We estimated effective population sizes (\emph{N\textsubscript{e}} )
using the linkage disequilibrium model with random mating (Waples \& Do,
2008) implemented in NeEstimator v2.1 (Do et al., 2014). We report
estimated \emph{N} \textsubscript{e} values using `Lowest Allele
Frequency Used' 5\% and 95\% confidence intervals generated by the
`Parametric method' for unrelated individuals for each population
cluster identified by STRUCTURE separately.

\emph{Population structure}

We characterized population genetic structure using the SNP dataset with
no missing data. We performed PCA and DAPC with the Adegenet 2.1.1
package (Jombart, 2008) in R. For the DAPC analysis, we first
conducted\emph{K} -means clustering and selected the number of clusters
based on the lowest Bayesian Information Criterion (BIC) value. We
performed cross-validation with the function xvalDapc to determine the
number of PCs to retain by calculating the lowest root mean squared
error. We then ran DAPC, retaining 20 PCs and 2 discriminant functions.

We used the Bayesian clustering algorithm in the program STRUCTURE
v2.3.4 to infer the number of population clusters (\emph{K} ) and the
proportion of individual membership assigned to each cluster
(\emph{q\textsubscript{k}} ). We used a burn-in of 500,000 steps
followed by 1,000,000 recorded steps, tracking the probability of the
data given\emph{K} (LnP(D)) to ensure that we ran the program long
enough for the values to stabilize. We used the admixture model, no
location priors, and assumed correlated allele frequencies (Falush,
Stephens, \& Pritchard, 2003). We performed a simulation with \emph{K}
from 1 to 7 with 10 replicates each and identified meaningful \emph{K}
values using the \selectlanguage{greek}Δ\selectlanguage{english}\emph{K} method (Evanno, Regnaut, \& Goutdet, 2005)
implemented in STRUCTURE HARVESTER v0.6.94 (Earl \& vonHoldt, 2012). We
combined replicate runs using CLUMPP v1.1.2 (Jakobsson \& Rosenberg,
2007).

Using the SNP dataset with 10\% missing data, we performed an AMOVA and
calculated pairwise \emph{F} \textsubscript{ST} values between
populations identified by STRUCTURE in GenAlEx v6.5 (Peakall \& Smouse,
2012), with 10,000 permutations to generate the null distribution. To
investigate local spatial structure, we performed a Mantel test using
ade4 v1.7 (Dray \& Dufour, 2007) in R on both the full and unrelated
dataset. We tested for a correlation between pairwise genotypic distance
and Euclidean geographic distance with 9,999 permutations to generate
the null distribution. We also generated a Mantel correlogram to test
for spatial autocorrelation between pairs of treeshrews at different
distance classes using GenAlEx. We first calculated pairwise linear
geographic and genotypic distances, and then used the `Spatial' option
with 9,999 permutations. We defined 7 distance classes (0.2, 1.0, 2.0,
5.0, 10.0, 15.0, and 18.0 km) based on Sturges's Rule (Sturges, 1926),
chosen to ensure sufficient comparisons within each class. Finally, we
calculated the average, median, and maximum geographic distances between
pairs of individuals in each kinship class corresponding to first,
second, third-order, and distant relatives (Table 1). To test for
significant differences between the means in each kinship class, we
performed a one-way ANOVA in R with a Tukey Honest Significant
Differences test and a Bonferroni correction for multiple comparisons
(Combs, Puckett, Richardson, Mims, \& Munshi-South, 2018).

\emph{Modeling migration}

We used the program MIGRATE-N v3.6 with its Bayesian implementation
(Beerli, 2005) to compare support for six different models of population
structure and migration (Figure 2). We used the phased pseudo-haplotype
sequence dataset for this analysis to take advantage of the higher
information content in DNA sequences relative to SNPs. Although our
workflow may generate chimeric sequences in instances where phase is
unresolved, we do not expect that this would affect model selection
(Andermann et al., 2018; Beerli \emph{pers. comm.} ). To ensure that
phase did not affect model inference, we ran the MIGRATE-N analysis
twice with different configurations of variants within haplotypes while
maintaining all other settings.

We compared the following models: 1) panmixia, 2) four populations (high
elevation MK, low elevation MK, low elevation MT and high elevation MT)
with bidirectional migration between adjacent pairs, 3) three
populations (high elevation MT, low elevation MT and all of MK) with
bidirectional migration between adjacent pairs, 4) three populations
with migration between all pairs, 5) two populations (high elevation MT
separate from all others) with bidirectional migration, and 6) two
populations with unidirectional migration from high elevation MT (Figure
2).

For model 2, we assigned individuals to populations based on their
sampling location: low elevation individuals \textless{} 2000 masl, and
high elevation [?] 2000 masl. For models 3, 4, 5, and 6, we assigned
individuals based on STRUCTURE output for \emph{K} = 3 and \emph{K} = 2,
respectively. We randomly selected five individuals from each population
cluster (\emph{n} = 10 haplotypes). We did not include all individuals
because for coalescent processes, increasing the sample size above this
does not necessarily improve accuracy, but substantially increases
computation time (Felsenstein, 2005). For each model, we ran two long
chains of 20,000,000 steps, sampled every 100 steps with 50,000 steps
per chain discarded as burn-in, and with four heated chains. To ensure
comparability across models, we ran the most complex model first and
used the same prior distributions and run parameters for all subsequent
models. We assessed chain mixing through acceptance ratios and ESS of
parameters and genealogies. We calculated log Bayes factors (LBF) and
model probability using the Bezier approximation of the marginal model
likelihood and the formula described in (Beerli \& Palczewski, 2010).

\textbf{Results}

\emph{DNA Sequencing}

We obtained mitochondrial genome sequences from 83 mountain treeshrew
individuals and 10 sequences from three congeners including the large
treeshrew (\emph{n} = 6), the small treeshrew (\emph{n} = 2) and the
ruddy treeshrew (\emph{n} = 2). Mitogenomes were sequenced to an average
depth of 50\selectlanguage{ngerman}×. Sequences were submitted to GenBank (Accession numbers
pending).

We sequenced UCEs from 80 mountain treeshrews (NCBI SRA number pending).
Each UCE-enriched library was sequenced with a mean of 2.3 million reads
(914,104--7,011,836), yielding a mean of 3,344 UCE loci (2,137--3,489)
per sample. The total number of UCE alignments that we used to generate
the pseudo-reference was 4,106, and the mean length was 495 base pairs
(149--2167 bps). After aligning reads to the pseudo-reference and
quality filtering, there were 7,861 SNPs including multiple SNPs per
locus. After removing loci with more than 10\% missing data across
individuals, 3,168 SNPs remained. The unlinked SNP dataset included
1,794 independent SNPs. Removing loci with missing data left 684
unlinked loci. In the phased pseudo-haplotype sequence alignment
dataset, after removing loci with less than one (114 loci) and more than
10 PIS (16 loci), 1,664 UCE alignments remained (Figure S1).

\emph{Kinship}

Pairwise kinship calculations revealed several groups of putatively
related individuals. After removing first-order relatives (\emph{n} = 22
with kinship \selectlanguage{english}[?] 0.18) from the dataset, 58 individuals remained,
including 33 from MT and 25 from MK.

\emph{Mitogenomes: genetic diversity, population structure, and
demographic inference}

There were 36 unique mitochondrial haplotypes in the dataset that
included close relatives (\emph{n} = 83), and 34 among the 58
``unrelated'' individuals. All subsequent analyses were performed with
the unrelated dataset. \emph{H} \textsubscript{d} was 0.977 (SD 0.008),
and \selectlanguage{greek}π \selectlanguage{english}was 0.00583 (SD 0.0006). Phylogenetic analyses recovered mountain
treeshrews as a monophyletic group with two deeply divergent but
sympatric lineages. The average number of nucleotide substitutions per
site between the two lineages is 0.013 (Figure S2). Outgroup
relationships were consistent with the topology in Roberts et al.
(2011). The BEAST dating analysis suggests that the lineages diverged
ca. 450,000 ybp (95\% Highest Posterior Density, HPD, 346,000--631,900
ybp) (Figure S3).

The median joining haplotype network (Figure 3) shows that the two
mountain treeshrew haplogroups are sympatric on both MT and MK. Three
haplotypes are found on both mountains (Table S3). Including related
individuals, haplogroups 1 and 2 are found in near equal proportion on
MK and MT (16 and 14 individuals, respectively), while haplogroup 1 is
more frequent on MT (46 out of 53 individuals) (Figure 4a). The AMOVA on
the unrelated dataset showed significant differentiation between the two
mountains (\emph{F} \textsubscript{ST} = 0.133, \emph{p} =.00812), with
13.3\% of variance accounted for by differences between mountains and
86.7\% within mountains. Dividing the population into high ([?] 2000 masl)
and low (\textless{} 2000 masl) elevation groups on each peak, 90.42\%
of the total variance is accounted for by within-group variation, and
9.58\% among (\emph{F} \textsubscript{ST} = 0.096, \emph{p =} 0.027).
Pairwise comparisons showed significant differences between high
elevation MK and low elevation MT (\emph{F} \textsubscript{ST} =
0.15,\emph{p} = 0.023) and high elevation MK and high elevation MT
(\emph{F} \textsubscript{ST} = 0.18, \emph{p} = 0.013).

Tajima's \emph{D} test was not significant, indicating a lack of
evidence for recent population contraction or expansion. In the Bayesian
skyline plot analysis, the 95\% HPD of the population change parameter
included zero, meaning that we cannot reject the hypothesis of zero
demographic changes in the last 60,000 years.

\emph{SNPs: Genetic diversity}

The nucleotide diversity of phased UCE pseudo-haplotypes (1,664
alignments concatenated end-to-end) for all 58 unrelated mountain
treeshrews was 0.001732 (SE 0.000022). Using the SNP dataset with no
missing data, average individual heterozygosity for all 80 individuals
was 0.2297 (SD 0.027), and the average inbreeding coefficient (\emph{F}
\textsubscript{IS}) was 0.0121 (SD 0.12). For the 58 unrelated
individuals, average heterozygosity was 0.2260 (SD 0.027), and\emph{F}
\textsubscript{IS} was 0.0187 (SD 0.12). Bartlett's test revealed that
the variances in observed and expected heterozygosity were not
significantly different (\emph{K\textsuperscript{2}} = 1.676,\emph{p} =
0.195). Average \emph{F} \textsubscript{IS} was higher on MK than MT,
but the difference was not significant (0.04 and 0.01 respectively,
Welch two-sample t-test \emph{p} = 0.2). Using the dataset with 10\%
missing data, the probability of two individuals having identical
genotypes assuming siblings are present (\emph{P\textsubscript{IDsib}} )
was 1.54 \selectlanguage{ngerman}× 10\textsuperscript{-199} .

\emph{Population Structure}

Both DAPC and STRUCTURE indicated that the most likely number of
population clusters was two and the second most likely was three, as
determined by BIC and the \selectlanguage{greek}Δ\selectlanguage{english}\emph{K} method, respectively. The
\selectlanguage{greek}Δ\selectlanguage{english}\emph{K}method is biased toward \emph{K} = 2 (Janes et al., 2017;
Campana, Hunt, Jones, \& White, 2011) and simulation studies have shown
that the mean probability (MeanLnP(K)) output from STRUCTURE performs
better in scenarios with high gene flow and low \emph{F}
\textsubscript{ST} (Latch, Dharmarajan, Glaubitz, \& Rhodes, 2006).
Because \emph{K} = 3 produced the highest MeanLnP(K) in the STRUCTURE
analysis (Table S4a), we consider this a relevant model and show the
proportion of individual membership in each cluster as defined by each
of the two analyses for both \emph{K} = 2 and \emph{K} = 3 (Figure 5).
We also ran STRUCTURE separately for individuals caught on MK and MT,
with settings described above except we ran simulations for \emph{K} =
1--5. We found no evidence of structure among MK individuals. MT
individuals were divided into two clusters - one with individuals [?] 2000
masl and one with individuals \textless{} 2000 masl, with individuals of
mixed ancestry at 2000 masl (Table S4b).

Cluster membership is mostly concordant between the DAPC and STRUCTURE
analyses, except STRUCTURE assigned mixed ancestry to many individuals
while DAPC did not. This is not unexpected as previous studies have
shown that DAPC may underestimate admixture (Frosch et al., 2014) while
STRUCTURE is more accurate at assigning mixed ancestry (Bohling, Adams
\& Waits, 2012). When \emph{K} = 2, individuals at 2000 and 2400 masl on
MT form a separate cluster from low elevation MT + MK (Figures 4 \& 5),
with mixed ancestry individuals at 2000 masl MT. For \emph{K} = 3, the
divisions are between high elevation MT, low elevation MT, and MK, with
individuals assigned mixed ancestry at low elevation MK (900 masl) and
2000 masl MT (Figure 5).

The PCA shows a similar pattern. Eigenvector 1 (7\% variation explained)
separates the two mountains, with overlap among individuals at 900 masl.
Eigenvector 2 (4\% variation explained) partially separates individuals
by elevation, with lower elevation individuals clustered together
(Figure 6). With \emph{K} = 2, after removing individuals that could not
be assigned to a STRUCTURE cluster (cutoff \emph{q}
\textsubscript{k}value \textless{} 0.6), \emph{F} \textsubscript{ST} is
0.05 (\emph{p} = 0.0001). The AMOVA showed 95\% of variation is
partitioned within clusters and 5\% between. With \emph{K} = 3, removing
individuals with\emph{q} \textsubscript{k} values \textless{}
0.6,\emph{F} \textsubscript{ST} between MK and low elevation MT was
0.035 (\emph{p} = 0.001), between MK and high elevation MT 0.092
(\emph{p} = 0.0001), and between low elevation MT and high elevation
MT\emph{F} \textsubscript{ST} = 0.065 (\emph{p} = 0.0005) (Tables 2a \&
2b). The AMOVA showed that 94\% of variation is distributed within
clusters, and 6\% among.

Including data for all 80 individuals, the Mantel test revealed a
significant, positive correlation between genotypic distance and
geographic distance (\emph{r} = 0.287, \emph{p} = 0.0001). Including the
58 unrelated individuals, the correlation was weaker but statistically
significant (\emph{r} = 0.05, \emph{p} \textless{} 0.0001). The
correlogram showed significant positive autocorrelation between
individuals at distances of 200 m and less (\emph{r} = 0.091, \emph{pr}
-rand [?] \emph{p} \emph{r} -data = 0.0001) and 200 m--1 km (\emph{r} =
0.036, \emph{p} \emph{r} -rand [?] \emph{p} \emph{r} -data = 0.0001);
autocorrelation was no longer significant at 2 km (\emph{r} = -0.001,
\emph{p} \emph{r} -rand [?] \emph{p} \emph{r} -data =0.598). At subsequent
distance classes (5, 10, 15, and 18 km), individuals have greater
genetic distance than expected at random; i.e., there is a significant
negative correlation between genetic and geographic distance (\emph{p}
\emph{r} -rand [?] \emph{r-} data = 0.009, 0.0001, 0.0001, 0.0001,
respectively) (Figure S4). The average geographic distance between pairs
of first-order relatives was 162.5 m, second order was 1.2 km, third
order was 4.8 km, and between distant or `unrelated' individuals was 12
km (Table 1). Differences between first and third, first and distant,
second and distant, and third and distant relatives were significant
(\emph{p} \textless{} 0.05).

\emph{Population and migration models}

Model 4 (Figure 2) was the best fit model as determined by Bayes Factors
in the MIGRATE-N analysis (Table S5). This result was consistent across
both runs of the program with alternative phasing. Model 4 divided the
population into three groups: high elevation MT, low elevation MT, and
MK, with high rates of bidirectional migration between all pairs. The
next best model was model 5, which divides the population into high
elevation MT and low elevation MT + MK with bidirectional migration
(Figure 2). Across models, the mean migration rate from high elevation
MT to MK was greater than from MK to MT (1.3--19\selectlanguage{ngerman}×, Table S6).

The results from NeEstimator suggest a larger effective population size
on MT (\textless{} 2000 masl + \selectlanguage{english}[?] 2000 masl) than MK despite the smaller
available habitat on MT (250 vs. 125, respectively) (Table 2 and
Supplemental Information).

\textbf{Discussion}

\emph{High gene flow across a heterogeneous landscape}

Seasonally stable temperatures and steep ecological gradients on
tropical mountains are predicted to limit dispersal across elevations
(Janzen, 1967; Ghalambor et al., 2006). This prediction is supported by
several studies that show limited gene flow and significant genetic
differentiation across elevations in taxa with continuous elevational
distributions (Bertrand et al., 2014; Mila, Warren, Heeb \& Thebaud,
2010; Gadek et al., 2018), sometimes with elevational parapatry between
sister species (Linck et al., 2019; Raxworthy et al. 2008a; DuBay \&
Witt, 2014). For example, Gadek et al. (2018) found cryptic population
genetic structure differentiating high from low elevation populations in
3 of 4 Peruvian songbird species with continuous elevational
distributions.

Our results are not consistent with a hypothesis of restricted gene flow
across elevational gradients. We report evidence of high gene flow
between MK and MT as well as between low and high elevations on both
peaks, indicating that neither the lowland habitat connecting the two
peaks, nor the steep elevational gradient across which mountain
treeshrews occur on each peak, has significantly limited effective
dispersal.

The mountain treeshrew population in KNP is best described as comprising
two or three clusters. The primary population division does not
correspond to the two peaks; rather, the summit region of MT (\selectlanguage{english}[?] 2000
masl) is distinct (Figure 5a). When dividing the population into three
clusters, individuals at low elevation MK show mixed ancestry with low
elevation MT (Figure 5b). Additionally, the MIGRATE-N analysis supports
the division of individuals into three population clusters (high
elevation MT, low elevation MT, and all of MK), with high migration
rates between all pairs (Figure 2).

The temperature lapse rate in KNP is -0.55degC per 100 meters of elevation
gain (Kitayama, 1992) which means that mountain treeshrews experience a
12.65degC average range in temperature between 900 and 3200 masl on MK.
This temperature range is higher than the thermal neutral zone for most
mammals that weigh \textless{} 1 kg (Khaliq, Hof, Prinzinger,
Bohning-Gaese, \& Pfenninger, 2014). On MT, mountain treeshrews
experience an 8.8degC range from 900 to 2500 masl (Camacho-Sanchez et al.,
2018). On MK, the species occupies four distinct vegetation zones, while
on MT it occupies three (Kitayama, 1992). If gene flow were restricted
due to limited elevational dispersal or selection against
cross-elevation migrants, we would expect to find greater genetic
differentiation on MK because of the broader elevational range and
associated diversity of environmental factors on this slope compared to
MT. However, the summit of MT was consistently recovered as the most
distinct population cluster while individuals caught along the entire
elevational gradient on MK form a single cluster (Figures 4 \& 5). This
suggests that elevation is not the primary variable influencing genetic
structure in KNP.

Tropical montane species tend to have narrower ecological tolerance and
elevational ranges than those in the temperate zone. McCain (2009)
reviewed 170 montane gradients across vertebrate taxa to test this
hypothesis and found that most groups showed the predicted pattern of
decreasing range size with decreasing latitude. However, in rodents
there was either no relationship or range size increased with declining
latitude. This finding could be explained by the presence of cryptic
species pairs or genetic differentiation separating low elevation
populations from high (den Tex, Thorington, Maldonado, \& Leonard, 2010;
Hinckley, Hawkins, Achmadi, Maldonado, \& Leonard, in review). However,
McCain (2009) hypothesized that rodents may be able to cope with the
lower temperatures associated with increasing altitude through
behavioral adaptations including burrowing. Further studies are
necessary to determine whether mountain treeshrews show behavioral or
phenotypic plasticity, adaptive phenotypic (e.g. Branch et al., 2017) or
genetic differentiation despite gene flow, or a generalist phenotype
(e.g. Feijo et al., 2019) that allows them to persist across such broad
environmental conditions which is rare among small mammals in this
landscape (Nor, 2001; Camacho-Sanchez et al., 2019).

\emph{Population genetic structure shaped by isolation-by-distance and
historical dynamics}

The population genetic pattern that we observe is partly consistent with
isolation-by-distance (IBD). The Mantel test showed a significant,
positive correlation between pairwise geographic distance and genetic
distance although the effect size is very small when individuals with
high kinship values are removed from the dataset. Individuals located
within 1 km have significantly lower pairwise genetic distances than
expected by chance, and this significance drops off by 2 km (Figure S4).
Additionally, pairwise \emph{F} \textsubscript{ST} between the
non-adjacent MT summit and MK clusters is greater than the value between
neighboring clusters (Table 2a). High rates of gene flow between
adjacent demes across the landscape, with relatively short dispersal
distances as suggested by the correlogram and ANOVA (Figure S4, Table
1), could have generated the clinal pattern we observe (Figure 5);
however, this does not explain the distinctiveness of the summit MT
cluster. The Euclidean distance between the lowest and highest sampled
points on MK (ca. 13.5 km) is greater than the distance between the
lowest and highest sampling points on MT (ca. 4.5 km), yet there is more
population genetic differentiation on MT. This indicates that the
structure we observe is not due to isolation-by-distance or
isolation-by-elevation alone, and that genetic similarity decays with
geographic distance at unequal rates in this landscape (Figure S5a).

Historical population dynamics, in addition to IBD, likely contributed
to the observed population genetic structure. Without data from other
locations in Borneo, it is difficult to determine what process(es)
generated the pattern. However, we suggest a plausible scenario given
known information about the relative ages of MT and MK and the degree of
divergence between the mountain treeshrew and its sister
species\emph{.}MT reached its current elevation earlier (ca. 11--7 Mya)
than MK (ca. 1 Mya) (Collenette, 1964; Hall, 1998; Liew et al. 2010).
This suggests that MT was available for colonization prior to the split
of mountain treeshrews from ruddy treeshrews ca. 4 Mya (Roberts et al.,
2011). If mountain treeshrews were resident on MT prior to the major
uplift of MK, and a second colonization event occurred later, this would
explain the signature of two population clusters. We find
higher-than-average genetic diversity among individuals at high
elevation MT despite its smaller habitat area (Figures 1a \& S5b), which
is consistent with our hypothesis that this region maintained a
relatively stable, or recently reduced, effective population size over
time relative to MK. Lower gene flow upslope to high elevation MT
relative to gene flow towards MK may have preserved the signature of
this cluster.

In the MIGRATE-N analysis, the estimated rate of migration upslope to
high elevation MT was lower than migration downslope and toward MK
(Table S6). We suggest that reduced gene flow to high elevation MT may
be related to shifts in the plant community that occur between 1450 masl
and the summit (van der Ent, Cardace, Tibbett, \& Echevarria, 2018).
Supporting this hypothesis, trapping success of mountain treeshrews and
other small mammals is low from 1500 to 1800 masl, and increases above
2000 masl (Camacho-Sanchez et al., 2019).

By contrast, lack of differentiation across MK could have been
influenced by an upslope shift at the mountain treeshrew's upper
elevational limit, enabled by climate warming and upslope shifts in
montane forest since the Last Glacial Maximum (LGM). During late
Quaternary glacial cycles, MK's summit alternated between ice-free
(during glacial minima) and ice-covered (during glacial maxima) periods
(Hall et al., 2008). After the LGM ca. 20,000 ybp, the ice began to melt
and by 9,200 ybp the summit was likely ice-free (Hall et al., 2008).
Upslope shifts in montane forest during this period of warming could
have enabled range expansion at high elevations by mountain
treeshrews\emph{.} Mountain treeshrews on MT likely did not experience a
concurrent upslope range shift since MT has a much lower summit than MK
which was never covered in ice (Hall et al., 2008). The lack of a
population expansion signature in the mountain treeshrew mitogenome data
could be explained by unrestricted gene flow between adjacent areas
during expansion (Pierce, Gutierrez, Rice, \& Pfennig, 2017). As
predicted for a recent expansion, we find lower-than-average genetic
diversity among high elevation MK individuals (\selectlanguage{english}[?]1600 masl) using
estimated effective migration surface modeling (Petkova, Novembre, \&
Stephens, 2016) to visualize genetic diversity on the landscape (Figure
S5b).

\emph{Mito-nuclear discordance}

The population genetic pattern inferred from our mitogenome data is
discordant with the nuclear SNP dataset, although it is not inconsistent
with a scenario of two colonization events to KNP. Phylogenetic analyses
revealed two divergent mitochondrial lineages within mountain
treeshrews; both lineages are found on both mountains, but haplogroups 1
and 2 are equally represented on MK while haplogroup 1 is more frequent
on MT (Figures 3 \& 4). As mentioned above, MT provided montane habitat
earlier than MK. If KNP were colonized a second time by mountain
treeshrews from another part of their range, such as the southern
portion of the Crocker Range, this would explain the presence of two
sympatric divergent lineages within Kinabalu Park. The greater frequency
of haplogroup 2 on MK could be explained by the closer geographic
proximity of MK to the Crocker Range (Figure 7) combined with
male-biased dispersal limiting the movement of haplogroup 2 from MK to
MT. Lack of recombination in the mitochondrial genome would have
retained the divergence between these the two lineages whereas
recombination in nuclear SNPs would result in genetic admixture between
the two groups. However, as noted above, the cluster on high elevation
MT has maintained a moderate level of differentiation.

As an alternative explanation, this pattern could have been caused by a
single colonization event of two sympatric lineages that diverged
elsewhere in Borneo, for example, due to isolation in interglacial
refugia and mixing during glacial maxima when montane forest was at its
maximum extent (Cannon, Morley, \& Bush, 2009; den Tex et al, 2010).
However, this scenario implies that the colonization of KNP by mountain
treeshrews would have occurred after the divergence between the two
lineages ca. 450,000 ybp, which is relatively recent compared to the age
of MT (at least 7 million years) and the age of the species (ca. 4
million years). Multiple colonization events to MK have been inferred in
other taxa, including plants in the genus \emph{Rhododendron} (Merckx et
al., 2015).

Gawin et al. (2014) documented a similar pattern in mountain blackeyes
(\emph{Chlorocharis emiliae} ) in Borneo; they found two divergent
mitochondrial haplogroups on MK, with one lineage sister to a lineage
found on Mt. Trus Madi, a mountain south of MK within the Crocker Range
(Figure 7). The pattern inferred from SNP data in a subsequent study was
not concordant, with a single lineage found on MK (Manthey et al.,
2017). This similar pattern may indicate a common colonization history
between mountain blackeyes and mountain treeeshrews. Future studies
should include broader geographic sampling of mountain treeshrews,
including individuals from across the Crocker Range, to test the
hypothesis of multiple colonization events and to determine the
phylogeographic history of this species in Borneo.

\emph{Utility of UCEs for fine-scale population genomics}

Since the original design of DNA hybridization probes to capture and
sequence UCEs for tetrapods (Faircloth et al., 2012), UCEs have proven
capable of resolving phylogenies at multiple timescales and across taxa
(e.g. Hawkins et al., 2016; McCormack et al., 2012; Smith et al., 2014;
Sun et al., 2014). UCEs can resolve phylogeographic patterns, population
structure, and demographic history of recently diverged species (e.g.
Harvey et al., 2016; Winker et al., 2018; Zarza et al., 2016). Here, we
show that sequence capture of ca. 5,000 UCEs yielded two highly
informative datasets suitable for landscape-scale population genomics.
These datasets were able to resolve patterns of population structure in
mountain treeshrews within KNP, an area of approximately 754
km\textsuperscript{2}. The SNP dataset provided sufficient statistical
power to identify individuals with high probability
(\emph{P\textsubscript{IDsib}} = 1.54 x 10\textsuperscript{-199} ) and
to identify putative family groups using pairwise kinship estimates.
This suggests the potential for UCEs to be used for population genomic
studies.

UCEs are valuable for studying species like the mountain treeshrew for
which few genomic resources are available. Although RAD-seq methods also
do not require reference genomes, UCE capture produces data with similar
information content and has several benefits over RAD-seq: 1) UCEs
enable direct comparison of inferences drawn from the same set of loci
across species, allowing conclusions to be drawn about the effects of
historical processes on diverse taxa (Lim et al., 2020), 2) UCE capture
offers repeatability such that studies can compare inferences for the
same species across time and geographic regions (Harvey et al., 2016),
and 3) UCE capture can be performed on low-quality DNA, including DNA
derived from historical museum specimens (Hawkins et al., 2016; Lim et
al., 2020; Lim \& Braun, 2016, Tsai et al., 2019).

\emph{Conservation Implications}

Anthropogenic GCC is one of the most significant threats to Earth's
biodiversity (Elsen \& Tingley, 2015; Bellard et al., 2014). Tropical
montane species are predicted to suffer the greatest effects because
they are likely to be endemic with restricted geographic ranges and
narrow climatic tolerance (Williams, Jackson, \& Kutzbach, 2007; La
Sorte \& Jetz, 2010). Several studies have documented upslope range
shifts in response to GCC-associated temperature and precipitation
changes (e.g. Freeman, Scholer, Ruiz-Gutierrez, \& Fitzparick, 2018;
Raxworthy et al., 2008b; Chen et al., 2009). These shifts can result in
reduced or fragmented habitat, decreasing genetic diversity and making
species vulnerable to extinction (Moritz et al., 2008; Williams et al.,
2007). In addition, in the case of locally adapted, highly structured
populations, increased dispersal and gene flow upslope could introduce
maladaptive genes (Weiss-Lehman \& Shaw, 2019). It is therefore
imperative to understand how tropical montane populations are structured
in order to predict, monitor, and mitigate the effects of GCC.

Here, we provide a foundational estimate of population structure,
genetic diversity, and historical gene flow for mountain treeshrews in
KNP. These data can be used to monitor GCC-induced population genetic
changes over time. Montane communities in KNP could experience an
upslope shift of ca. 490 m by 2100 CE (Still, 1999; Camacho et al.,
2018), assuming mild Intergovernmental Panel on Climate Change scenarios
(IPCC 2013, www.ipcc.ch/report/ar5/wg1/). Although the factors that
limit the mountain treeshrew at its lower elevational boundary are
unknown, assuming that the species tracks the predicted 490 m upslope
shift - whether because of climatic limitations or ecological
interactions with lowland species expanding upslope - we predict that it
will experience range contraction. The species already occupies the
upper elevational limits within KNP, so an upslope shift in the lower
bound of its distribution could not be countered with expansion at its
upper limit. The lack of strong population structure across elevations
means that upslope dispersal of lower elevation mountain treeshrews on
MK will likely not increase extinction risk by introducing maladaptive
genetic diversity. However, reduction in available habitat could make
the species vulnerable.

We also predict that mountain treeshrews would maintain connectivity
between MK and MT. However, the Crocker Range has few peaks above 1400
masl, and connectivity between KNP and the rest of the Crocker Range
could be severed (Figure 7). This highlights the importance of KNP as a
future refugium for montane species, as it contains the highest peak in
the region and the greatest high-elevation forested area. Conservation
efforts should focus on protecting forest habitat at 900 masl to
facilitate gene flow and preserve genetic diversity.

\textbf{Acknowledgements}

L. Olson, C. Thompson, and UMMZ provided \emph{Tupaia splendidula} DNA
extracts. S. Lindley and the FDA donated a MiSeq sequencing kit. T.
Giarla gracefully shared his SNP calling scripts. S. Castaneda Rico,
M-E. Ochirbat, M. Venkatraman, N. R. McInerney, A. M. Kearns, and A. W.
Kaganer provided helpful comments on the manuscript. K. Helgen and R.
Fleischer provided logistical support. Chien Lee provided the mountain
treeshrew photograph. This work was funded by the Spanish Government
(CGL2010-21524, CGL2017-86068-P), the Smithsonian Institution
(\emph{Building the Framework of Biodiversity Science: Next Generation
Phylogenomics} Smithsonian Grand Challenges Awards 2012 - 2014), the
National Science Foundation (1547168, 1717498), and the Department of
Biology at George Mason University. Logistical support was provided by
Laboratorio de Sistemas de Informacion Geografica y Teledeteccion de la
Estacion Biologica de Donana and Donana ICTS-RBD.

\textbf{Data Accessibility}

-DNA sequences: Genbank accessions XX-XX; NCBI SRA: XX

-UCE pseudo-haplotype multiple sequence alignments available on FigShare

\textbf{Author Contributions}

MTRH, MCS, and JAL organized field work; MTRH and MCS conducted field
work; MTRH, JAL, JEM, and LDP designed the study; JAL, JEM, and LLR
provided resources for the study; MTRH and LDP did laboratory work; LDP
wrote the manuscript and performed SNP analyses; MTRH performed
phylogenetic analyses; MCS contributed figures; TRW, HCL, and MGC
assisted in bioinformatics; MGC and JAW developed vcf2aln; all authors
participated in revisions and acceptance of the final version of the
manuscript.

\textbf{References}

Andermann, T., Fernandes, A. M., Olsson, U., Topel, M., Pfeil, B.,
Oxelman, B., Aleixo, A., Faircloth, B. C., \& Antonelli, A. (2018).
Allele phasing greatly improves the phylogenetic utility of
ultraconserved elements. \emph{Systematic Biology} , 68(1), 32--46.
https://doi.org/10.1093/sysbio/syy039

Beerli, P., \& Palczewski, M. (2010). Unified framework to evaluate
panmixia and migration direction among multiple sampling
locations.\emph{Genetics} , 185(1), 313--326.
https://doi.org/10.1534/genetics.109.112532

Beerli, P. (2005). Comparison of Bayesian and maximum-likelihood
inference of population genetic parameters. \emph{Bioinformatics} ,
22(3), 341--45. https://doi.org/10.1093/bioinformatics/bti803

Bellard, C., Leclerc, C., Leroy, B., Bakkenes, M., Veloz, S., Thuiller,
W., \& Courchamp, F. (2014). Vulnerability of biodiversity hotspots to
global change. \emph{Global Ecology and Biogeography} , 23(12),
1376--1386. https://doi.org/10.1111/geb.12228

Bertrand, J. A. M., Bourgeois, Y. X. C, Delahaie, B., Duval, T.,
Garcia-Jimenez, R., Cornuault, J., Heeb, P., Mila, B., Pujol, B., \&
Thebaud, C. (2014). Extremely reduced dispersal and gene flow in an
island bird. \emph{Heredity} , 112(2), 190--196.
https://doi.org/10.1038/hdy.2013.91

Bohling, J. H., Adams, J. R., \& Waits, L. P. (2013). Evaluating the
ability of Bayesian clustering methods to detect hybridization and
introgression using an empirical red wolf data set. \emph{Molecular
Ecology} , 22(1), 74--86. https://doi.org/10.1111/mec.12109

Bolger, A. M., Lohse, M., \& Usadel, B. (2014). Trimmomatic: A flexible
trimmer for Illumina sequence data. \emph{Bioinformatics} , 30(15),
2114--2120. https://doi.org/10.1093/bioinformatics/btu170

Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C., Xie, D.,
\ldots{}, Drummond, A. J. (2014). BEAST 2: A software platform for
Bayesian evolutionary analysis. \emph{PLoS Computational Biology} ,
10(4), e1003537. https://doi.org/10.1371/journal.pcbi.1003537

Branch, C. L., Jahner, J. P., Kozlovsky, D. Y., Parchman, T. L. \&
Pravosudov, V. V. (2017). Absence of population structure across
elevational gradients despite large phenotypic variation in mountain
chickadees (\emph{Poecile gambeli} ). \emph{Royal Society Open Science}
, 4(3), 170057. https://doi.org/10.1098/rsos.170057

Camacho-Sanchez, M., Quintanilla, I., Hawkins, M. T. R., Tuh Yit Yu, F.,
Wells, K., Maldonado, J. E. \& Leonard, J. A. (2018). Interglacial
refugia on tropical nountains: Novel insights from the summit rat
(\emph{Rattus baluensis} ), a Borneo mountain endemic. \emph{Diversity
and Distributions} , 24,1252-1266. https://doi.org/10.1111/ddi.12761

Camacho-Sanchez, M., \& Hawkins, M. T. R., Tuh Yit Yu, F., Maldonado, J.
E., \& Leonard, J. A. (2019). Endemism and diversity of small mammals
along two neighboring Bornean mountains. \emph{PeerJ} , 7,e7858.
https://doi.org/10.7717/peerj.7858

Campana, M. G., Hunt, H. V., Jones, H., \& White, J. (2011). CorrSieve:
Software for summarizing and evaluating Structure output.\emph{Molecular
Ecology Resources} , 11(2), 349--52.
https://doi.org/10.1111/j.1755-0998.2010.02917.x

Cannon, C. H., Morley, R. J., \& Bush, A. B. G. (2009). The current
refugial rainforests of Sundaland are unrepresentative of their
biogeographic past and highly vulnerable to
disturbance.\emph{Proceedings of the National Academy of Sciences} ,
106(27), 11188--11193. https://doi.org/10.1073/pnas.0809865106

Castillo Vardaro, J. A., Epps, C. W., Frable, B. W., \& Ray, C. (2018).
Identification of a contact zone and hybridization for two subspecies of
the American pika (\emph{Ochotona princeps} ) within a single protected
area. \emph{PLOS ONE} , 13(7), e0199032.
https://doi.org/10.1371/journal.pone.0199032

Castresana, J. (2000). Selection of conserved blocks from multiple
alignments for their use in phylogenetic analysis. \emph{Molecular
Biology and Evolution} , 17(4), 540--552.

Chen, I. C., Shiu, H. J., Benedick, S., Holloway, J. D., Chey, V. K.,
Barlow, H. S., Hill, J. K., \& Thomas, C. D. (2009). Elevation increases
in moth assemblages over 42 years on a tropical
mountain.\emph{Proceedings of the National Academy of Sciences} ,
106(5), 1479--1483. https://doi.org/10.1073/pnas.0809320106

Cheviron, Z. A., \& Brumfield, R. T. (2009). Migration-selection balance
and local adaptation of mitochondrial haplotypes in rufous-collared
sparrows (\emph{Zonotrichia capensis} ) along an elevational
gradient.\emph{Evolution} , 63(6), 1593--1605.
https://doi.org/10.1111/j.1558-5646.2009.00644.x

Chin, L., Moran, J. A., \& Clarke, C. (2010). Trap geometry in three
giant montane pitcher plant species from Borneo is a function of tree
shrew body size. \emph{New Phytologist} , 186(2), 461--470.
https://doi.org/10.1111/j.1469-8137.2009.03166.x

Collenette, P. (1964). A short account of the geology and geological
history of Mt. Kinabalu. \emph{Proceedings of the Royal Society of
London Series B, Biological Sciences} , 161, 56-63.

Combs, M., Puckett, E. E., Richardson, J., Mims, D., \& Munshi\selectlanguage{english}-South, J.
(2017). Spatial population genomics of the brown rat (\emph{Rattus
norvegicus} ) in New York City. \emph{Molecular Ecology} , 27(1),
83--98. https://doi.org/10.1111/mec.14437

Conomos, M. P., Reiner, A. P., Weir, B. S., \& Thornton, T. A. (2016).
Model-free estimation of recent genetic relatedness. \emph{The American
Journal of Human Genetics} , 98(1), 127--48.
https://doi.org/10.1016/j.ajhg.2015.11.022

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E.,
DePristo, M. A., \ldots{} 1000 Genomes Project Analysis Group. (2011).
The variant call format and VCFtools. \emph{Bioinformatics} , 27 (15),
2156--58. https://doi.org/10.1093/bioinformatics/btr330

den Tex, R. J., Thorington, R., Maldonado, J. E., \& Leonard, J. A.
(2010). Speciation dynamics in the SE Asian tropics: putting a time
perspective on the phylogeny and biogeography of Sundaland tree
squirrels, Sundasciurus.~\emph{Molecular Phylogenetics and Evolution}
,~\emph{55} (2), 711-720. doi: 10.1016/j.ympev.2009.12.023

Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., \&
Ovenden, J. R. (2014). NeEstimator v2: Re-implementation of software for
the estimation of contemporary effective population size (Ne) from
genetic data. \emph{Molecular Ecology Resources} , 14(1), 209--14.
https://doi.org/10.1111/1755-0998.12157

Dray, S., \& Dufour, A. B. (2007). The Ade4 package: Implementing the
duality diagram for ecologists. \emph{Journal of Statistical Software} ,
1(4), 1-20. https://doi.org/10.18637/jss.v022.i04

Drummond, A. J., \& Rambaut, A. (2007). BEAST: Bayesian Evolutionary
Analysis by Sampling Trees. \emph{BMC Evolutionary Biology} , 7(1), 214.
https://doi.org/10.1186/1471-2148-7-214

DuBay, S. G., \& Witt, C. C. (2014). Differential high-altitude
adaptation and restricted gene flow across a mid-elevation hybrid zone
in Andean tit-tyrant flycatchers. \emph{Molecular Ecology} , 23(14),
3551--3565. https://doi.org/10.1111/mec.12836

Earl, D. A., \& vonHoldt, B. M. (2012). STRUCTURE HARVESTER: A website
and program for visualizing STRUCTURE output and implementing the Evanno
method. \emph{Conservation Genetics Resources} , 4(2), 359--361.
https://doi.org/10.1007/s12686-011-9548-7

Elsen, P. R., \& Tingley, M. W. (2015). Global mountain topography and
the fate of montane species under climate change. \emph{Nature Climate
Change} , 5(8), 772--776. https://doi.org/10.1038/nclimate2656

Emmons, L. H. (2000). \emph{Tupai: A field study of Bornean treeshrews}
. Berkeley, CA: University of California Press.

Evanno, G., Regnaut, S., \& Goutdet, J. (2005). Detecting the number of
clusters of individuals using the software STRUCTURE: A simulation
study. \emph{Molecular Ecology} , 14(8), 2611--2620.
https://doi.org/10.1111/j.1365-294X.2005.02553.x

Excoffier, L., Foll, M., \& Petit, R. J. (2008). Genetic consequences of
range expansions. \emph{Annual Review of Ecology, Evolution, and
Systematics} , 40(1), 481--501.
https://doi.org/10.1146/annurev.ecolsys.39.110707.173414

Excoffier, L., \& Lischer, H. E. L. (2010). Arlequin suite ver. 3.5: A
new series of programs to perform population genetics analyses under
Linux and Windows. \emph{Molecular Ecology Resources} , 10(3), 564--567.
https://doi.org/10.1111/j.1755-0998.2010.02847.x

Falush, D., Stephens, M., \& Pritchard, J. K. (2003). Inference of
population structure using multilocus genotype data: Linked loci and
correlated allele frequencies. \emph{Genetics} , 164(4), 1567--1587.

Faircloth, B. C. (2016). PHYLUCE is a software package for the analysis
of conserved genomic loci. \emph{Bioinformatics} , 32(5), 786--88.
https://doi.org/10.1093/bioinformatics/btv646

Faircloth, B. C., McCormack, J. E., Crawford, N. G., Harvey, M. G.,
Brumfield, R. T., \& Glenn, T. C. (2012). Ultraconserved elements anchor
thousands of genetic markers spanning multiple evolutionary
timescales.\emph{Systematic Biology} , 61(5), 717--726.
https://doi.org/10.1093/sysbio/sys004

Feeley, K. J., Stroud, J. T., \& Perez, T. M. (2017). Most `global'
reviews of species' responses to climate change are not truly
global.\emph{Diversity and Distributions} , 23(3), 231--234.
https://doi.org/10.1111/ddi.12517

Feijo, A., Wen, Z., Cheng, J., Ge, D., Xia, L., \& Yang, Q. (2019).
Divergent selection along elevational gradients promotes genetic and
phenotypic disparities among small mammal populations. \emph{Ecology and
Evolution} , 9(12), 7080--7095. https://doi.org/10.1002/ece3.5273

Felsenstein, J. (2005). Accuracy of coalescent likelihood estimates: Do
we need more sites, more sequences, or more loci? \emph{Molecular
Biology and Evolution} , 23(3), 691--700.
https://doi.org/10.1093/molbev/msj079

Freeman, B. G., Scholer, M. N., Ruiz-Gutierrez, V., and Fitzpatrick, J.
W. (2018). Climate change causes upslope shifts and mountaintop
extirpations in a tropical bird community. \emph{Proceedings of the
National Academy of Sciences} , 115(47), 11982--11987.
https://doi.org/10.1073/pnas.1804224115

Frosch, C., Kraus, R. H. S., Angst, C., Allgower, R., Michaux, J.,
Teubner, J., \& Nowak, C. (2014). The genetic legacy of multiple beaver
reintroductions in central Europe. \emph{PLOS ONE} , 9(5), e97619.
https://doi.org/10.1371/journal.pone.0097619

Gadek, C. R., Newsome, S. D., Beckman, E. J., Chavez, A. N., Galen, S.
C., Bautista, E., \& Witt, C. C. (2018). Why are tropical mountain
passes `low' for some species? Genetic and stable-isotope tests for
differentiation, migration and expansion in elevational generalist
songbirds. \emph{Journal of Animal Ecology} , 87(3), 741--753.
https://doi.org/10.1111/1365-2656.12779

Gawin, D. F., Rahman, M. A., Ramji, M. F. S., Smith, B. T., Lim, H. C.,
Moyle, R. G., \& Sheldon, F.H. (2014). Patterns of avian diversification
in Borneo: The case of the endemic mountain black-eye
(\emph{Chlorocharis emiliae} ). \emph{The Auk} , 131(1), 86--99.
https://doi.org/10.1642/AUK-13-190.1

Ghalambor, C. K., Huey, R. B., Martin, P. R., Tewksbury, J. J., \& Wang,
G. (2006). Are mountain passes higher in the tropics? Janzen's
hypothesis revisited. \emph{Integrative and Comparative Biology} ,
46(1), 5--17. https://doi.org/10.1093/icb/icj003

Giarla, T. C., \& Esselstyn, J. A. (2015). The challenges of resolving a
rapid, recent radiation: Empirical and simulated phylogenomics of
Philippine shrews. \emph{Systematic Biology} , 64(5), 727-740.
https://doi.org/10.1093/sysbio/syv029

Gill, B. A., Kondratieff, B. C., Casner, K. L., Encalada, A. C.,
Flecker, A. S., Gannon, D. G., \ldots{}, Funk, W. C. (2016). Cryptic
species diversity reveals biogeographic support for the `mountain passes
are higher in the tropics' hypothesis. \emph{Proceedings of the Royal
Society B: Biological Sciences} , 283, 1832.
https://doi.org/10.1098/rspb.2016.0553

Goldberg, C. S., \& Waits, L. P. (2010). Quantification and reduction of
bias from sampling larvae to infer population and landscape genetic
structure. \emph{Molecular Ecology Resources} 10(2), 304--313.
https://doi.org/10.1111/j.1755-0998.2009.02755.x

Gueuning, M., Suchan, T., Rutschmann, S., Gattolliat, J., Jamsari, J.,
Kamil, A. I., \ldots{}, Alvarez, N. (2017). Elevation in tropical sky
islands as the common driver in structuring genes and communities of
freshwater organisms. \emph{Scientific Reports} , 7(1), 16089.
https://doi.org/10.1038/s41598-017-16069-y

Hall, R. (1998). The plate tectonics of Cenozoic SE Asia and the
distribution of land and sea. In R. Hall \& J. D. Holloway
(Eds.)\emph{Biogeography and the geological evolution of SE Asia} . (pp.
99-131). Backhuys Publishers: Leiden, The Netherlands.

Hall, R., Cottam, M., Suggate, S., Tongkul, F., Sperber, C., \& Batt, G.
E. (2009). \emph{The Geology of Mount Kinabalu} . Sabah Parks: Sabah,
Malaysia.

Han, K. H., Sheldon, F. H., \& Stuebing, R. B. (2000). Interspecific
relationships and biogeography of some Bornean tree shrews (Tupaiidae
Tupaia), based on DNA hybridization and morphometric
comparisons.\emph{Biological Journal of the Linnean Society} 70(1),
1--14. https://doi.org/10.1006/bijl.1999.0358

Harvey, M. G., Smith, B. T., Glenn, T. C., Faircloth, B. C., \&
Brumfield, R. T. (2016). Sequence capture versus restriction site
associated DNA sequencing for shallow systematics. \emph{Systematic
Biology} , 65(5), 910--924. https://doi.org/10.1093/sysbio/syw036

Hawkins, M. T. R. (2019). Family Tupaiidae. In R. A. Mittermeier \& D.
E. Wilson (Eds.) \emph{Handbook of the Mammals of the World Volume 8:
Insectivores, Sloths, and Colugos} (pp. 42--69). Lynx Edicions:
Barcelona, Spain.

Hawkins, M. T. R., Leonard, J. A., Helgen, K. M., McDonough, M. M.,
Rockwood, L. L., \& Maldonado, J. E. (2016). Evolutionary history of
endemic Sulawesi squirrels constructed from UCEs and mitogenomes
sequenced from museum specimens. \emph{BMC Evolutionary Biology} ,
16(1). https://doi.org/10.1186/s12862-016-0650-z

Hinckley, A., Hawkins, M. T. R., Achmadi, A., Maldonado, J. E., \&
Leonard, J. A. Ancient divergence driven by geographic isolation and
ecological adaptation in forest dependent Sundaland tree squirrels. (in
review) \emph{Frontiers in Ecology and Evolution} .

Jakobsson, M., \& Rosenberg, N. A. (2007). CLUMPP: A cluster matching
and permutation program for dealing with label switching and
multimodality in analysis of population structure.\emph{Bioinformatics}
, 23(14), 1801--1806. https://doi.org/10.1093/bioinformatics/btm233

Janes, J. K., Miller, J. M., Dupuis, J. R., Malenfant, R. M., Gorrell,
J. C., Cullingham, C. I., \& Andrew, R. L. (2017). The K = 2
conundrum.\emph{Molecular Ecology} , 26(14), 3594--3602.
https://doi.org/10.1111/mec.14187

Janzen, D. H. (1967). Why mountain passes are higher in the
tropics.\emph{The American Naturalist} , 101(919): 233--249.
https://doi.org/10.1086/282487

Jombart, T. (2008). Adegenet: An R package for the multivariate analysis
of genetic markers. \emph{Bioinformatics} , 24(11), 1403--1405.
https://doi.org/10.1093/bioinformatics/btn129

Katoh, K., Misawa, K., Kuma, K., \& Miyata, T. (2002). MAFFT: A novel
method for rapid multiple sequence alignment based on fast fourier
transform. \emph{Nucleic Acids Research} , 30(14), 3059--3066.
https://doi.org/10.1093/nar/gkf436

Khaliq, I., Hof, C., Prinzinger, R., Bohning-Gaese, K., \& Pfenninger,
M. (2014). Global variation in thermal tolerances and vulnerability of
endotherms to climate change. \emph{Proceedings of the Royal Society
Biological Sciences} , 281, 20141097.
https://doi.org/10.1098/rspb.2014.1097

Kitayama, K. (1992). An altitudinal transect study of the vegetation on
Mount Kinabalu, Borneo \emph{Vegetatio} , 102(2), 149--171.
www.jstor.org/stable/20046210

Kumar, S., Stecher, G., \& Tamura, K. (2016). MEGA7: Molecular
evolutionary genetics analysis version 7.0 for bigger
datasets.\emph{Molecular Biology and Evolution} , 33(7), 1870--1874.
https://doi.org/10.1093/molbev/msw054

La Sorte, F. A., \& Jetz, W. (2010). Projected range contractions of
montane biodiversity under global warming. \emph{Proceedings of the
Royal Society B: Biological Sciences} , 277(1699), 3401--10.
https://doi.org/10.1098/rspb.2010.0612

Lanfear, R., Calcott, B., Ho, S. Y. W., \& Guindon, S. (2012).
PartitionFinder: Combined selection of partitioning schemes and
substitution models for phylogenetic analyses. \emph{Molecular Biology
and Evolution} , 29(6), 1695--1701.
https://doi.org/10.1093/molbev/mss020

Latch, E. K., Dharmarajan, G., Glaubitz, J. C., \& Rhodes, O. E. (2006).
Relative performance of Bayesian clustering software for inferring
population substructure and individual assignment at low levels of
population differentiation. \emph{Conservation Genetics} , 7(2),
295--302. https://doi.org/10.1007/s10592-005-9098-1

Leigh, J. W., \& Bryant, D. (2015). Popart: Full-feature software for
haplotype network construction. \emph{Methods in Ecology and Evolution}
, 6(9), 1110--1116. https://doi.org/10.1111/2041-210X.12410

Lenoir, J., \& Svenning, J. C. (2015). Climate-related range shifts - a
global multidimensional synthesis and new research
directions.\emph{Ecography} , 38(1), 15--28.
https://doi.org/10.1111/ecog.00967

Li, H. (2013). Aligning sequence reads, clone sequences and assembly
contigs with BWA-MEM. \emph{ArXiv} 1303.3997v1.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N.,
\ldots{}, 1000 Genome Project Data Processing Subgroup. (2009). The
Sequence Alignment/Map format and SAMtools. \emph{Bioinformatics} ,
25(16), 2078--2079. https://doi.org/10.1093/bioinformatics/btp352

Librado, P., \& Rozas, J. (2009). DnaSP v5: A software for comprehensive
analysis of DNA polymorphism data. \emph{Bioinformatics} , 25(11),
1451--1452. https://doi.org/10.1093/bioinformatics/btp187

Liew, T. S., Schilthuizen, M., \& bin Lakim, M. (2010). The determinants
of land snail diversity along a tropical elevational gradient:
Insularity, geometry and niches. \emph{Journal of Biogeography} , 37(6),
1071--1078. https://doi.org/10.1111/j.1365-2699.2009.02243.x

Lim, H. C., and Braun, M. J. (2016). High-throughput SNP genotyping of
historical and modern samples of five bird species via sequence capture
of ultraconserved elements. \emph{Molecular Ecology Resources} , 16(5),
1204--1223. https://doi.org/10.1111/1755-0998.12568

Lim, H. C., Shakya, S. B., Harvey, M. G., Moyle, R. G., Fleischer, R.
C., Braun, M. J., \& Sheldon, F. H. (2020). Opening the door to greater
phylogeographic inference in Southeast Asia: Comparative genomic study
of five codistributed rainforest bird species using target capture and
historical DNA. \emph{Ecology and Evolution,} 10, 3222--3247.
https://doi.org/10.1002/ece3.5964

Linck, E., Freeman, B. G., \& Dumbacher, J. P. (2019). Speciation and
gene flow across an elevational gradient in New Guinea
kingfishers.\emph{BioRxiv} , 589044. https://doi.org/10.1101/589044

Manichaikul, A., Mychaleckyj, J. C., Rich, S. S., Daly, K., Sale, M., \&
Chen, W. M. (2010). Robust relationship inference in genome-wide
association studies. \emph{Bioinformatics} , 26 (22), 2867--2873.
https://doi.org/10.1093/bioinformatics/btq559

Manthey, J. D., Moyle, R. G., Gawin, D. F., Rahman, M. A., Ramji, M. F.
S., \& Sheldon, F. H. (2017). Genomic phylogeography of the endemic
mountain black-eye of Borneo (\emph{Chlorocharis emiliae} ): Montane and
lowland populations differ in patterns of Pleistocene
diversification.\emph{Journal of Biogeography} , 44(10), 2272--2283.
https://doi.org/10.1111/jbi.13028

Mason, N. A., Olvera-Vital, A., Lovette, I. J., \& Navarro-Siguenza, A.
G. (2018). Hidden endemism, deep polyphyly, and repeated dispersal
across the Isthmus of Tehuantepec: Diversification of the white-collared
seedeater complex (Thraupidae: \emph{Sporophila torqueola}
).\emph{Ecology and Evolution} , 8(3): 1867--1881.
https://doi.org/10.1002/ece3.3799

McCain, C. M. 2009. Vertebrate range sizes indicate that mountains may
be `higher' in the tropics. \emph{Ecology Letters} , 12(6), 550--560.
https://doi.org/10.1111/j.1461-0248.2009.01308.x

McCormack, J. E., Faircloth, B. C., Crawford, N. G., Gowaty, P. A.,
Brumfield, R. T., \& Glenn, T. C. (2012). Ultraconserved elements are
novel phylogenomic markers that resolve placental mammal phylogeny when
combined with species-tree analysis. \emph{Genome Research} , 22(4),
746--754. https://doi.org/10.1101/gr.125864.111

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K.,
Kernytsky, A., \ldots{}, DePristo, M. A. (2010). The Genome Analysis
Toolkit: A MapReduce framework for analyzing next-generation DNA
sequencing data. \emph{Genome Research,} 20(9), 1297--1303.
https://doi.org/10.1101/gr.107524.110

McKinney, G. J., Waples, R. K., Seeb, L. W., \& Seeb, J.E. 2017.
Paralogs are revealed by proportion of heterozygotes and deviations in
read ratios in genotyping-by-sequencing data from natural
populations.\emph{Molecular Ecology Resources} , 17(4): 656--669.
https://doi.org/10.1111/1755-0998.12613

Merckx, V. S. F. T., Hendriks, K. P., Beentjes, K. K., Mennes, C. B.,
Becking, L. E., Peijnenburg, K. T. C. A., \ldots{}, Schilthuizen, M.
(2015). Evolution of endemism on a young tropical mountain.\emph{Nature}
, 524(7565), 347--350. https://doi.org/10.1038/nature14949

Mila, B., Warren, B. H., Heeb, P., \& Thebaud, C. (2010). The geographic
scale of diversification on islands: Genetic and morphological
divergence at a very small spatial scale in the Mascarene grey white-eye
(Aves: \emph{Zosterops borbonicus} ). \emph{BMC Evolutionary Biology} ,
10(158). https://doi.org/10.1186/1471-2148-10-158.

Moritz, C., Patton, J., Conroy, C. J., Parra, J. L, White, G. C., \&
Beissinger, S.R. (2008). Impact of a century of climate change on small
mammal communities in Yosemite National Park, USA.
\emph{Science,}322(5899), 258--261.
https://doi.org/10.1126/science.1162547

Moritz, C. (1994). Defining `evolutionarily significant units' for
conservation. \emph{Trends in Ecology \& Evolution} , 9(10), 373--375.
https://doi.org/10.1016/0169-5347(94)90057-4

Muenchow, J., Dieker, P., Kluge, J., Kessler, M., \& von Wehrden, H.
(2018). A review of ecological gradient research in the tropics:
Identifying research gaps, future directions, and conservation
priorities. \emph{Biodiversity and Conservation} , 27(2), 273--285.
https://doi.org/10.1007/s10531-017-1465-y

Nor, S.M. (2001). Elevational diversity patterns of small mammals on
Mount Kinabalu, Sabah, Malaysia. \emph{Global Ecology and Biogeography}
, 10(1): 41--62. https://doi.org/10.1046/j.1466-822x.2001.00231.x

Oswald, J. A., Harvey, M. G., Remsen, R. C., Foxworth, D. U., Cardiff,
S. W., Dittmann, D. L., \ldots{}, Brumfield, R. T. (2016). Willet be one
species or two? A genomic view of the evolutionary history
of\emph{Tringa semipalmata} . \emph{The Auk} , 133(4), 593--614.
https://doi.org/10.1642/AUK-15-232.1

Payne, J., Francis, C. M., \& Phillipps, K. (2016). \emph{A Field Guide
to the Mammals of Borneo.} The Sabah Society: Kota Kinabalu, Sabah,
Malaysia.

Peakall, R., \& Smouse. P. E. (2012). GenAlEx 6.5: Genetic analysis in
Excel. Population genetic software for teaching and research-an
update.\emph{Bioinformatics} , 28(19), 2537--2539.
https://doi.org/10.1093/bioinformatics/bts460

Petkova, D., Novembre, J., \& Stephens, M. (2016). Visualizing spatial
population structure with estimated effective migration
surfaces.\emph{Nature Genetics} , 48(1): 94--100.
https://doi.org/10.1038/ng.3464

Phillipps, Q., \& Phillipps, K. (2016). \emph{Phillips' field guide to
the mammals of Borneo and their ecology.} Kota Kinabalu Natural History
Publications: Kota Kinabalu, Sabah, Malaysia.

Pierce, A. A., Gutierrez, R., Rice, A. M., \& Pfennig, K.S. (2017).
Genetic variation during range expansion: effects of habitat novelty and
hybridization. \emph{Proceedings of the Royal Society B: Biological
Sciences} , 284(1852). https://doi.org/10.1098/rspb.2017.0007

Polato, N. R., Gill, B. A., Shah, A. A., Gray, M. M., Casner, K. L.,
Barthelet, A., \ldots{}, Zamudio, K. R. (2018). Narrow thermal tolerance
and low dispersal drive higher speciation in tropical
mountains.\emph{Proceedings of the National Academy of Sciences} ,
115(49), 12471--12476. https://doi.org/10.1073/pnas.1809326115

Pritchard, J. K., Stephens, M., \& Donnelly, P. (2000). Inference of
population structure using multilocus genotype data. \emph{Genetics} ,
155(2), 945--959.

Raxworthy, C. J., Pearson, R. G., Zimkus, B. M., Reddy, S., Deo, A. J.,
Nussbaum, R. A., \& Ingram, C. M. (2008a). Continental speciation in the
tropics: contrasting biogeographic patterns of divergence in
the\emph{Uroplatus} leaf-tailed gecko radiation of
Madagascar.\emph{Journal of Zoology} , 275(4), 423--440.
https://doi.org/10.1111/j.1469-7998.2008.00460.x

Raxworthy, C. J., Pearson, R. G., Rabibisoa, N., Rakotondrazafy, A. M.,
Ramanamanjato, J. B., Raselimanana, A. P., \ldots{}, Stone, D. A.
(2008b). Extinction vulnerability of tropical montane endemism from
warming and upslope displacement: A preliminary appraisal for the
highest massif in Madagascar. \emph{Global Change Biology} , 14(8),
1703--1720. https://doi.org/10.1111/j.1365-2486.2008.01596.x

Roberts, T. E., Lanier, H. C., Sargis, E. J., \& Olson, L. E. (2011).
Molecular phylogeny of treeshrews (Mammalia: Scandentia) and the
timescale of diversification in Southeast Asia. \emph{Molecular
Phylogenetics and Evolution} , 60(3), 358--372.
https://doi.org/10.1016/j.ympev.2011.04.021

Ronquist, F., \& Huelsenbeck, J. P. (2003). MrBayes 3: Bayesian
phylogenetic inference under mixed models. \emph{Bioinformatics} ,
19(12): 1572--1574. https://doi.org/10.1093/bioinformatics/btg180

Schmitz, J., Ohme, M., \& Zischler, H. (2000). The complete
mitochondrial genome of \emph{Tupaia belangeri} and the phylogenetic
affiliation of Scandentia to other Eutherian orders. \emph{Molecular
Biology and Evolution} , 17(9), 1334--43.
https://doi.org/10.1093/oxfordjournals.molbev.a026417

Smith, B. T., Harvey, M. G., Faircloth, B. C., Glenn, T. C., \&
Brumfield, R. T. (2014). Target capture and massively parallel
sequencing of ultraconserved elements for comparative studies at shallow
evolutionary time scales. \emph{Systematic Biology} , 63(1), 83--95.
https://doi.org/10.1093/sysbio/syt061

Still, C. J., Foster, P. N., \& Schneider, S. H. (1999). Simulating the
effects of climate change on tropical montane cloud
forests.\emph{Nature} , 398(6728), 608--610.
https://doi.org/10.1038/19293

Sturges, H. (1926). The choice of a class-interval. \emph{Journal of the
American Statistical Association,} 21, 65--66.

Sun, K., Meiklejohn, K. A., Faircloth, B. C., Glenn, T. C., Braun, E.
L., \& Kimball, R. T. (2014). The evolution of peafowl and other taxa
with ocelli (eyespots): a phylogenomic approach. \emph{Proceedings of
the Royal Society B: Biological Sciences} , 281(1790).
https://doi.org/10.1098/rspb.2014.0823

Tsai, W. L. E., Mota-Vargas, C., Rojas-Soto, O., Bhowmik, R., Liang, E.
Y., Maley, J. M., \ldots{}, McCormack, J. E. (2019). Museum genomics
reveals the speciation history of Dendrortyx wood-partridges in the
mesoamerican highlands. \emph{Molecular Phylogenetis and Evolution} ,
136, 29--34. https://doi.org/10.1016/j.ympev.2019.03.017

van der Ent, A., Cardace, D., Tibbett, M., \& Echevarria, G. (2018).
Ecological implications of pedogenesis and geochemistry of ultramafic
soils in Kinabalu Park (Malaysia). \emph{Catena} 160, 154--169.
https://doi.org/10.1016/j.catena.2017.08.015

Waits, L. P., Luikart, G., \& Taberlet, P. (2001). Estimating the
probability of identity among genotypes in natural populations: Cautions
and guidelines. \emph{Molecular Ecology} , 10(1), 249--256.
https://doi.org/10.1046/j.1365-294X.2001.01185.x

Waples, R. S., \& Do, C. (2008). Ldne: A program for estimating
effective population size from data on linkage
disequilibrium.\emph{Molecular Ecology Resources} , 8(4), 753--756.
https://doi.org/10.1111/j.1755-0998.2007.02061.x

Weiss-Lehman, C., \& Shaw, A. K. (2019). Spatial population structure
determines extinction risk in climate-induced range shifts. \emph{The
American Naturalist} , 706259. https://doi.org/10.1086/706259

Williams, J. W., Jackson, S. T., \& Kutzbach, J. E. (2007). Projected
distributions of novel and disappearing climates by 2100
AD.\emph{Proceedings of the National Academy of Sciences} , 104(14),
5738--5742. https://doi.org/10.1073/pnas.0606292104

Winker, K., Glenn, T. C., \& Faircloth, B.C. (2018). Ultraconserved
elements (UCEs) illuminate the population genomics of a recent,
high-latitude avian speciation event. \emph{PeerJ} , 6, e5735.
https://doi.org/10.7717/peerj.5735

Yu, H. T. (1995). Patterns of diversification and genetic population
structure of small mammals in Taiwan. \emph{Biological Journal of the
Linnean Society} , 55(1): 69--89.
https://doi.org/10.1111/j.1095-8312.1995.tb01050.x

Zarza, E., Faircloth, B. C., Tsai, W. L. E., Bryson R. W., Klicka, J.,
\& McCormack, J. E. (2016). Hidden histories of gene flow in highland
birds revealed with genomic markers. \emph{Molecular Ecology} , 25(20):
5144--5157. https://doi.org/10.1111/mec.13813

\textbf{Table 1.} Geographic distances between pairs of individuals with
different levels of estimated relatedness based on analysis with the
KING software. ANOVA and Tukey Honest Significant Differences test
showed significant differences in distances between all pairs except
first and second-order relatives and second and third-order relatives
(\emph{p} \textless{} 0.05).\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
\textbf{Kinship} & \textbf{Relatedness} & \textbf{n} & \textbf{Average
distance (m)} & \textbf{Median distance (m)} & \textbf{Max distance
(m)}\tabularnewline
\midrule
\endhead
\textgreater{} 0.18 & First order & 42 & 162.5 & 100.8 &
570.1\tabularnewline
0.177 - 0.0884 & Second order & 56 & 1247 & 322.8 & 25850\tabularnewline
0.0883 - 0.0442 & Third order & 112 & 4819 & 793.0 &
26490\tabularnewline
\textless{} 0.044 & Distant & 2932 & 12310 & 15960 &
29430\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2.} Effective population sizes and pairwise\emph{F}
\textsubscript{ST} of population clusters with \textbf{a)}\emph{K} = 3
and \textbf{b)} \emph{K} = 2. MK = Mount Kinabalu, MT = Mount
Tambuyukon. \emph{N} \textsubscript{e} estimates are on the diagonal
with 95\% CI in parenthesis; \emph{F} \textsubscript{ST}estimates are
below the diagonal, with associated \emph{p} -values above the diagonal.

\textbf{a} .\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
\emph{K} = 3 & MK & MT \textless{} 2000masl & MT \textgreater{}
2000masl\tabularnewline
\midrule
\endhead
MK (n=22) & 125 (105--152) & 0.00120 & 0.00010\tabularnewline
MT \textless{} 2000masl (n=19) & 0.035 & 202 (157--282) &
0.00050\tabularnewline
MT \textgreater{} 2000masl (n=11) & 0.092 & 0.065 & 48
(40--59)\tabularnewline
\bottomrule
\end{longtable}

\textbf{b.}\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\begin{minipage}[b]{0.30\columnwidth}\raggedright\strut
\emph{K} = 2\strut
\end{minipage} & \begin{minipage}[b]{0.30\columnwidth}\raggedright\strut
MK + MT \textless{} 2000masl\strut
\end{minipage} & \begin{minipage}[b]{0.30\columnwidth}\raggedright\strut
MT \textgreater{} 2000masl\strut
\end{minipage}\tabularnewline
\midrule
\endhead
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
MK + MT \textless{} 2000masl (n=36)\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
180 (160--205)\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0001\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
MT \textgreater{} 2000masl (n=18)\strut
\end{minipage} & \begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
0.050\strut
\end{minipage} & \begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
57 (52--63)\strut
\end{minipage}\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figure 1. a (left)} Map of mountain treeshrew distribution
(inset modified from IUCN 2019), and map of Kinabalu Park, Sabah,
Borneo. Park boundaries are demarcated by dashed lines, transects by
black lines, and sampling locations by white circles. Shading indicates
lower and upper portion of mountain treeshrew habitat (900--2000 masl
and \textgreater{}2000 masl, respectively) \textbf{b (right).} Image of
a mountain treeshrew and a pitcher plant (\emph{Nepenthes lowii} ),
Kinabalu National Park (Photo credit: Chien C. Lee). The two species
have a mutualistic relationship in which mountain treeshrews feed on the
sugary secretions provided by the plant and in turn provide the plant
nitrogen through feces (Chin, Moran, \& Clarke, 2010).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2.} Population structure and migration models evaluated
using MIGRATE-N. Model ranks are based on Bayes factors; log marginal
likelihood values are listed in Table S4. MT, Mt. Tambuyukon; MK, Mt.
Kinabalu; High [?] 2000 masl; Low \textless{} 2000 masl.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3.} Median joining network of 34 mitogenome sequences.
Dashed lines represent the number of base pair differences between
haplotypes except in cases where the number of differences exceeds 40.
Colors correspond to the two mountains (MT, orange; MK, blue). Note that
the two haplogroups are not shown to scale and are separated by 186 bp
substitutions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure} \textbf{4.} Elevational profiles of Mt. Kinabalu and Mt.
Tambuyukon with the distribution of mitochondrial haplogroups per
elevation shown below each transect line and SNP clusters above.
Mitogenome pie charts indicate, for each elevation, the number of
treeshrews sampled with a haplotype from mitochondrial haplogroup 1
(light grey) and 2 (dark grey). SNP pie charts indicate for each
elevation the proportion of ancestry assigned to cluster 1 (light grey)
and cluster 2 (dark grey).

\textbf{Hosted file}

\verb`image5.emf`  available at \url{https://authorea.com/users/313978/articles/444406-high-gene-flow-across-heterogeneous-tropical-montane-environments-in-a-bornean-endemic-small-mammal}

\textbf{Figure 5.} Cluster membership based on DAPC and STRUCTURE
analyses for \textbf{a)} \emph{K} = 2 and \textbf{b)} \emph{K} = 3. Each
horizontal line represents a single individual with shading showing how
much of each's ancestry can be attributed to each cluster. Individuals
are arranged from high elevation Mt. Tambuyukon to low elevation Mt.
Tambuyukon followed by low elevation Mt. Kinabalu to high elevation Mt.
Kinabalu.

\textbf{Hosted file}

\verb`image6.emf`  available at \url{https://authorea.com/users/313978/articles/444406-high-gene-flow-across-heterogeneous-tropical-montane-environments-in-a-bornean-endemic-small-mammal}

\textbf{Figure 6.} PCA Plot with individuals caught on Mt. Kinabalu
shaded in blue and Mt. Tambuyukon in orange. PC1 explains 7\% of
variance; PC2 explains 4\%.

\textbf{Hosted file}

\verb`image7.emf`  available at \url{https://authorea.com/users/313978/articles/444406-high-gene-flow-across-heterogeneous-tropical-montane-environments-in-a-bornean-endemic-small-mammal}

\textbf{Figure 7.} Estimated mountain treeshrew habitat in the year 2100
CE assuming mild IPCC scenarios. Light grey shows [?] 900 masl, which is
current mountain treeshrew habitat; dark grey indicates [?] 1400 masl,
which is potential mountain treeshrew habitat in 2100 assuming mild
climate change as projected by IPCC. Protected areas are demarcated by
dashed lines. Transects sampled in this study are shown in black, with
sampling locations indicated by white circles.

\textbf{Hosted file}

\verb`image8.emf`  available at \url{https://authorea.com/users/313978/articles/444406-high-gene-flow-across-heterogeneous-tropical-montane-environments-in-a-bornean-endemic-small-mammal}

\selectlanguage{english}
\FloatBarrier
\end{document}

