\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[T2A]{fontenc}
\usepackage[ngerman,polish,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Drainage basins serve as multiple glacial refugia for alpine habitats in
the Sierra Nevada Mountains, California}



\author[1]{Yi-Ming Weng}%
\author[2]{David Kavanaugh}%
\author[3]{Sean Schoville}%
\affil[1]{University of Wisconsin Madison Graduate School}%
\affil[2]{California Academy of Sciences}%
\affil[3]{University of Wisconsin Madison}%


\vspace{-1em}



  
  \date{May 23, 2024}


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





\selectlanguage{english}
\begin{abstract}
The evolutionary histories of alpine species are often directly
associated with responses to glaciation. Deep divergence among
populations and complex patterns of genetic variation have been inferred
as consequences of persistence within glacier boundaries (i.e. on
nunataks), while shallow divergence and limited genetic variation is
assumed to result from expansion from large refugia at the edge of ice
shields (i.e. massifs de refuge). However, for some species, dependence
on specific microhabitats could profoundly influence their spatial and
demographic response to glaciation, and such a simple dichotomy may
obscure the localization of actual refugia. In this study, we use the
Nebria ingens complex (Coleoptera: Carabidae), a water-affiliated ground
beetle lineage, to test how drainage basins are linked to their observed
population structure. By analyzing mitochondrial COI gene sequences and
genome-wide single nucleotide polymorphisms, we find that the major
drainage systems of the Sierra Nevada Mountains in California best
explain the population structure of the N. ingens complex. In addition,
we find that an intermediate morphotype within the N. ingens complex is
the product of historical hybridization of N. riversi and N. ingens in
the San Joaquin basin during glaciation. This study highlights the
importance of considering ecological preferences in how species respond
to climate fluctuations and provides an explanation for discordances
that are often observed in comparative phylogeographic studies.%
\end{abstract}%



\sloppy


\textbf{Introduction}

The geographical distribution and diversification patterns of alpine
species are commonly associated with historical glaciation cycles
(Wallis, Waters, Upton, \& Craw, 2016). During periodic glacial climate
fluctuations in the Pliocene and Pleistocene epochs, alpine species
diversified as a consequence of population isolation, punctuated
demographic change, and distributional range shifts (Knowles, 2001;
Jansson \& Dynesis, 2002; Schoville and Roderick, 2009).
Lineage-specific patterns of genetic structure and genetic diversity are
often linked to the location of glacial refugia and subsequent processes
of recolonization after the retreat of glaciers (Sch\selectlanguage{ngerman}önswetter, Stehlik,
Holderegger, \& Tribsch, 2005). For example, European alpine species
that retain substantial genetic diversity in a heterogeneous pattern
across the alpine zone have been associated with persistence in
high-elevation ice-free refugia (i.e. nunataks) during glaciation
(Stehlik, Blattner, Holderegger, \& Bachmann, 2002; Schneeweiss, \&
Schoenswetter, 2011). More generally, species that possess deep and
extensive lineage diversity are often inferred to have persisted in
multiple refugia, where isolated populations undergo strong genetic
drift due to the smaller population size (Schoville \& Roderick 2010;
Homburg et al., 2013; Rovito \& Schoville 2017). In contrast,
populations that show shallower genetic structure across alpine
distributions are inferred to have persisted in large, interconnected
populations outside of the ice shields during glacial periods (i.e.
massifs de refuge, or periglacial refugia). These two scenarios, at
their extremes referred to as 'nunatak survival' and 'massifs de
refuge', have been tested as competing hypotheses in many different
alpine species (Wachter et al., 2016; Rovito \& Schoville 2017;
Kosi\selectlanguage{polish}ń\selectlanguage{english}ski, S\selectlanguage{polish}ę\selectlanguage{english}kiewicz, Walas, Boraty\selectlanguage{polish}ń\selectlanguage{english}ski, \& Dering, 2019; Pan, H\selectlanguage{ngerman}ülber,
Willner, \& Schneeweiss, 2020). Both hypotheses are supported in
different study systems, with emerging genomic datasets often favoring
more complex population histories (Lohse, Nicholls, \& Stone, 2011).

Although this dichotomy of glacial survival considers the refugia being
either outside or within the boundary of glaciation, which seems
plausible for many alpine species, it might be insufficient for
reconstructing the location of glacial refugia for species that live in
specific habitats (Theissinger et al., 2013). Many alpine species
require specific habitat requirements that are unlikely to be commonly
found within or at the edge of large ice sheets. For example, the
discordance in population structure of two co-distributed sedge species
from the Rocky Mountains has been associated with their different
microhabitat preferences (Massatti \& Knowles 2014). Interestingly,
microhabitat partitioning is not rare for alpine species due to strong
environmental gradients and interspecific interactions (Gereben, 1995;
Kleckova, Konvicka, \& Klecka, 2014; Kikvidze et al., 2015). Subtle
differences in microhabitat preference could be responsible for the
discordant phylogeographic histories (Alvarez et al., 2009), even in
co-distributed, interacting species (DeChaine \& Martin, 2006). In
addition, as specific microhabitat requirements further constrain the
spatial extent of suitable habitat in glacial refugia, the population
structure of microhabitat specialists tends to show deeper
diversification (Schoville \& Roderick, 2010 ; Rovito \& Schoville,
2017). Therefore, species' microhabitat could be a principal driver of
population structure and evolutionary history, and in this study, we
test whether microhabitat usage of a cold-specialized alpine species
shapes lineage diversity patterns.

To test this hypothesis, we examine the population history of alpine
ground beetles in the \emph{Nebria ingens} species complex (Coleoptera:
Carabidae), because they possess the characteristics of microhabitat
specialization, host independence, limited dispersal ability, and strong
genetic differentiation (Schoville, Roderick, \& Kavanaugh, 2012).
Members of the \emph{N. ingens} complex are restricted to high-elevation
riparian, waterfall and 'seep' habitat (formed by melting snow and
glaciers), where they forage as generalist predators and scavengers (of
mostly invertebrates). While their habitat utilization and distribution
is not host dependent, they are patchily distributed due to their narrow
habitat preferences and limited dispersal ability. Like other
alpine\emph{Nebria} , members of the \emph{N. ingens} complex has
atrophied hind wings that impede flight, restricting gene flow and
distributional range size (Kavanaugh 1985). In fact, the \emph{N.
ingens} complex is limited to the Sierra Nevada Mountains of California
and provides an interesting case study of recent alpine speciation
(Schoville et al., 2012). The deep divergence between \emph{N. ingens}
and \emph{N. riversi} , estimated roughly around a million years ago,
suggests divergence in allopatry, yet several populations with
intermediate phenotypes, and discordant mitochondrial and nuclear data,
imply a more complex evolutionary history within this species complex.

Here, we reexamine the evolutionary history of the \emph{Nebria
ingens}complex using more extensive population sampling and genome-wide
genetic markers, to identify how glacial refugia and post-glacial
recolonization affect contemporary genetic structure. As population
structure is the product of both contemporary gene flow and historical
divergence (Paun, Schönswetter, Winkler, Consortium, \& Tribsch, 2008),
evolutionary history can be complicated to interpret. Conventionally, we
can apply coalescent-based approaches on multilocus data to
simultaneously estimate the demographic parameters such as effective
population size, divergence time, and migration rate (Carling, Lovette,
\& Brumfield, 2010; Qu et al., 2012; Rougemont \& Bernatchez, 2018).
However, alpine species typically possess complex population structure
that is difficult to model (due to the large number of free parameters)
without strong simplifying assumptions (e.g. Rovito and Schoville,
2018). Here, we employ a hypothesis testing framework based on
population structure models (\textbf{Table 1} ) to test competing
scenarios of 1) a single or a few core glacial refugia, 2)
drainage-associated refugia, or 3) survival within the extent of the
glacial ice sheets ('nunatak survival'). For example, one or a few
ancestral populations are expected if the entire species complex shifted
into periglacial refugia and remained interconnected during glaciation.
The glacial survival hypothesis (nunatak) would be supported if
ancestral populations and observed genetic diversity patterns are
directly associated with known ice free areas within the glacier
boundary (Rood, Burbank, \& Finkel, 2011e; \textbf{Figure 1} ). Finally,
persistence in drainage basins is proposed, as we expect ancestral
populations of the \emph{N. ingens}complex were dependent on these
geographical features to obtain suitable microhabitats.

\textbf{Materials and Methods}

\emph{Study System and testing for premating barriers to gene flow}

In previous genetic work, Schoville et al. (2012) revealed that there
are two major genetic lineages within the \emph{N. ingens} complex
associated with the allopatric species \emph{N. ingens} and \emph{N.
riversi} , but morphologically and genetically intermediate populations
suggested some gene flow. The morphologically intermediate populations
were readily distinguished by color pattern, as they were dull green,
whereas \emph{N. riversi} was bright metallic green and \emph{N.
ingens}was matte black. Here we use the same convention to group
individuals into three morphotypes. In order to understand whether these
morphotypes are capable of genetic exchange, we used a simple mating
experiment to assess premating isolation among the three morphotypes
(\emph{N. riversi} , intermediate, and \emph{N. ingens} ). Reduced
copulation across the three morphotypes would indicate limited
interspecific gene flow. In this experiment, samples from three
populations (sites 1, 8, and 27;\textbf{Figure 1} , \textbf{Table 2} )
were chosen to represent \emph{N. riversi} , the intermediate
morphotype, and \emph{N. ingens} , respectively. However, due to
conservation concerns regarding the small local population sizes of
beetles in the field, our experiment was constrained by sample size to
10-20 beetles per population. All beetles were collected in July 2018,
maintained at \textasciitilde{}5\textsuperscript{o}C in a refrigerator
with a mealworm diet prior to experiments, and all the trials were
conducted within a week following field collection. For each trial, one
male and one female were placed within a large yogurt cup (10 cm in
diameter and 4 cm in depth) for 30 minutes. All experimental
combinations of morphotypes were conducted (for both males and females),
including experiments within each morphotype as controls, to determine
the rate of successful copulation. A successful copulation was defined
as the male mounting the female, with direct observation of his
genitalia penetrating the female. Unfortunately, subsequent attempts to
rear the progeny from the crosses failed in the lab.

\emph{Genetic Sampling, DNA extraction and Sequencing}

In this study, we included 566 beetles from 27 populations of
the\emph{N. ingens} complex in the Sierra Nevada, California
(\textbf{Figure 1} ). The sample beetles were collected with permits
granted by Sequoia and Kings National Parks (Study\# SEKI-0091),
Yosemite National Park (Study\# YOSE-00093), and the California Fish and
Game Department (\#SC-006997). Genomic DNA from these samples was
extracted using the QIAamp 96 DNA Blood Kit (QIAGEN) or DNeasy Blood \&
Tissue Kits (QIAGEN) for both genotyping by sequencing (GBS) and
mitochondrial COI gene sequencing. In addition to 219 samples from the
study of Schoville et al. (2012), 104 beetles were selected to sequence
an approximately 645 base pair region (after trimming the primer
sequences) of the mitochondrial COI gene using the universal arthropod
primers Jerry (5'-CAACATTTATTTTGATTTTTTGG-3') and Pat (5'-TCCAATGCA
CTAATCTGCCATATTA-3') (Simon et al. 1994).). The resulting dataset of 323
COI gene sequences were aligned for downstream analyses. For nuclear
genomic analyses, 381 samples were selected for genotyping by sequencing
(GBS) (Elshire et al, 2011), with 3 individuals replicated. The
restriction enzyme ApeKI was used to reduce the genome to an appropriate
size distribution and 150 base pair single-end Illumina sequencing
libraries were prepared by the Biotechnology Center of University of
Wisconsin-Madison. All genetic data are deposited at NCBI (see Data
Accessibility statement).

\emph{Processing of Genomic Data}

Genome-wide single nucleotide polymorphisms (SNPs) were called from the
GBS data using Stacks v1.48 (Catchen, Hohenlohe, Bassham, Amores, \&
Cresko, 2013). To generate a large set of reference loci (ustacks), a
minimum depth of three reads was used to construct stacks of loci and
two mismatches were allowed between two stacks (alleles). To build the
catalog of loci among individuals (cstacks), two mismatches were allowed
between loci. The loci were then filtered using the following criteria:
loci were removed if (1) present in fewer than five populations, or (2)
present in fewer than 30\% of the individuals in a population. Before
calling genotypes, further criteria were used to remove loci. First, the
loci where more than 75\% of individuals had a confounded match
(multiple loci match to single catalog locus in an individual), or with
likelihood less than -8.0, were removed. Second, SNPs with an error rate
higher than 0.1 were dropped, with the error rate determined by the
three replicated individuals. This resulted in a \emph{semi-filtered
dataset} comprised of 311 individuals and 21,238 SNPs. In order to
reduce any potential confounding effects from linkage disequilibrium or
selection among SNPs, we further filtered this dataset such that only
one SNP with the least missing data was retained for each GBS locus, and
SNPs were removed that violated Hardy-Weinberg equilibrium at the
population level based on a p-value of 0.001. This
final\emph{fully-filtered} \emph{dataset} comprised 311 individuals and
12,498 SNPs, and is used for downstream analyses that require
independent loci.

\emph{Phylogenetic tree reconstruction}

For the mitochondrial COI dataset, the alignment of 323 sequences was
used to construct a Bayesian phylogenetic tree using the software Beast2
(Bouckaert et al, 2014), with out-group sequences of \emph{N. purpurata}
, \emph{N. meanyi} , \emph{N. vandykei} , \emph{N. praedicta} , and
\emph{N. metallica} downloaded from GenBank (accession numbers:
MN346054, KU641246, KU641248, JQ711382, KU875540). The model HKY+G was
selected as the best substitution model based on model tests conducted
with the software MEGA v7 (Kumar, Stecher, \& Tamura, 2016). A clock
rate for the COI gene was set to 0.0113 per site per million years based
on the best available COI mutation rate estimate for beetles, in this
case related carabid beetles (Andújar, Serrano, \& Gómez-Zurita, 2012).
Ten million steps of a Markov chain Monte Carlo (MCMC) were run to
estimate model parameters and trees, with samples drawn every 1000
steps. Two separate runs were conducted. The posterior samples were
assessed for convergence using TRACER v1.7 (Rambaut, Drummond, Xie,
Baele, \& Suchard, 2018) and each parameter value was assessed to ensure
a high effective sample size (\textgreater{}200). Additionally, a
statistical parsimony cladistic analysis, or TCS network, was
reconstructed using Popart-1.7 (Clement, Snell, Walker, Posada, \&
Crandall, 2002; Leigh \& Bryant, 2015). Based on the divisions in the
TCS network, haplotypes were assigned to northern or southern
haplogroups.

A population tree based on genome-wide SNP data (12,498
SNPs,\emph{fully-filtered} \emph{dataset} ) was constructed using
maximum likelihood in the program IQtree version 1.6.12 (Nguyen,
Schmidt, Von Haeseler, \& Minh, 2015). The substitution model GTR+ASC
was selected and the branch support was estimated based on 1,000
bootstrap replicates. The output tree was then visualized and colored
using FigTree v1.4.1 (Rambaut, 2012). Additionally, we inferred the
species tree using the program SNAPP (Bryant, Bouckaert, Felsenstein,
Rosenberg, \& RoyChoudhury, 2012), implemented in BEAST2 version 2.6.2.
The\emph{fully-filtered} \emph{dataset} was used, but due to
computational constraints, one random individual was selected from each
of eight populations. Individuals were chosen from populations with
relatively pure ancestry coefficients to avoid violating assumptions in
the species tree approach, namely the effect of gene flow from other
populations (see the result of sNMF analysis; Bryant et al., 2012,
Stoltz et al., 2019). SNAPP was run for 30 million Markov chain Monte
Carlo (MCMC) generations, with the custom program SNAPPER employed to
run SNAPP at faster speed (Stoltz et al., 2019). Two independent runs
were conducted and the output trees were visualized after 10\% burnin
using DensiTree version 2.01 (Bouckaert 2010).

\emph{Estimate of Ancestry Coefficients and Population Structure}

Sparse non-negative matrix factorization (sNMF) from the R package LEA
v1.4.0 (Frichot, Mathieu, Trouillon, Bouchard, \& François, 2014) was
used to reconstruct ancestry coefficients among sampled individuals in
the \emph{fully-filtered dataset} . To evaluate the most likely number
of ancestral populations (K), we computed the cross-entropy criterion
for K=1 to K=20 with ten replicates each. Since the lowest cross-entropy
was not clear, we present multiple levels of K to interpret the ancestry
coefficients (see results). The summed ancestry coefficients for each
population were plotted in geographical pie charts to visualize
geographical population structure (\textbf{Figure S1} ).

\emph{Testing spatial clustering by genetic cline}

To test whether the genetic clusters from sNMF are the consequence of
natural genetic clusters or isolation by distance (IBD) with
discontinuous sampling, we used the R package conStruct version 1.0.4
(Bradburd, Coop, \& Ralph, 2018). We compare the non-spatial (cluster)
and spatial (IBD) models for K=1 to K=10 ancestral populations. The
input allele frequencies of SNPs from \emph{fully-filtered dataset} were
calculated using VCFtools version 0.1.17 (Danecek et al., 2011), and the
geographical distances among populations were calculated using R package
geodist version 0.0.4 (Karney, 2013). In the conStruct package, we ran
cross-validation analysis through K=3 to K=10 with ten replicates and
100,000 iterations of the MCMC chains for each K. The predictive
accuracies were estimated with 90\% training and 10\% test data.
Two-sample t-tests were used to test differences between spatial and
non-spatial models for each K. Moreover, we also consider the preferred
models by examining the layer contributions from K=1 to K=10 using the
construct function, with two MCMC chains of 100,000 iterations
(Bradburd, Coop, \& Ralph, 2018).

\emph{Estimate of demographic history}

In order to better understand the demographic history of the \emph{N.
ingens} complex, we employed Stairway plot version 2 (Liu and Fu 2015,
Liu 2020) to reconstruct the historical demography of focal populations
using the \emph{semi-filtered dataset} . Similar to SNAPP, the eight
populations with pure ancestry coefficients were used to minimize any
bias of recent gene flow. The input folded site allele frequency
spectrum was estimated for each population using the custom program
easySFS (\emph{https://github.com/isaacovercast/easySFS}). We used a
point mutation rate of 2.8×\(10^{-9}\) per site per generation,
which was estimated based on \emph{Drosophila} and is close to the
nuclear mutation rate of other non-social insects (Keightley, Ness,
Halligan, \& Haddrill, 2014, Keightley et al., 2015, Oppold and
Pfenninger 2017). Other parameters were set to default.

\emph{Testing refugia models with TMRCA of paired populations}

The inferred ancestral clusters from sNMF and conStruct were further
tested as refugial populations by comparing the divergence times of
paired populations. Under a drainage refugia hypothesis, we expect to
find older divergence times between populations from different drainages
systems than within drainage systems, as the two populations retreated
to separate glacial refugia and were isolated for a longer period of
time. For the preferred ancestral clustering model (K=6), we assigned
populations to major drainage basins in the Sierra Nevada based on the
locations of their collection sites (\textbf{Table 2} ). If populations
had nearly equal distance from two different drainages, they were
assigned to both drainages. As several alternative ancestral clustering
models (K=4 and K=5) from sNMF were difficult to distinguish from K=6,
we examined the effect of this analysis under these alternative
clustering models (\textbf{Figure S1} , \textbf{Table 1} ). We grouped
populations from Tuolumne and Merced drainages into the same refugium
(K=5), and further grouped populations from the Kern and Kaweah
drainages into single refugium (K=4). Finally, we considered a three
refugia model, where populations were grouped within the same
morphotype.

Paired populations were chosen by grouping collecting sites with their
nearest four to seven geographical neighbors, including sites from the
same and/or different drainages or morphotypes. The time of most recent
common ancestor (TMRCA) of the paired populations was estimated using
the joint site frequency spectrum in \selectlanguage{greek}δ\selectlanguage{english}a\selectlanguage{greek}δ\selectlanguage{english}i version 2.0.3 (Gutenkunst,
Hernandez, Williamson, \& Bustamante, 2009). Since there is a potential
bias in using filtered genotype data (Warmuth and Ellegren 2019), we
estimated the two-dimensional allele frequency spectrums (2d-sfs)
directly from the GBS reads. The 2d-sfs estimation from GBS reads
requires a reference genome, so we generated a de novo assembly
of\emph{N. ingens} from a single male individual using 130x coverage
Illumina NovaSeq pair-end reads. We used SPAdes version 3.14.1
(Bankevich et al., 2012) to assemble the genome (more details regarding
genome assembly are described in the Supplementary Methods). The GBS
reads were mapped to the reference genome using STAR version 2.7.3a and
used to calculate the optimized 2d-sfs using angsd version 0.931, after
filtering out the low quality reads using default quality control
settings (Dobin et al., 2013, Korneliussen, Albrechtsen, \& Nielsen,
2014). The TMRCA, as well as the migration rate and effective population
size of each population, were approximated with the 2d-sfs using
dadi\_pipeline version 3.1.4 (Gutenkunst et al., 2009, Portik et al.,
2017). In the dadi\_pipeline program, we selected the divergence with
continuous symmetric migration model and ran this 10 times in a 4-run
(10, 20, 30, 40 replications) or 5-run (10, 20, 30, 40, 50 replications)
optimization, with grid size starting from the number that was rounded
up from the larger allele size to the nearest ten, until the likelihoods
and optimized parameters reached a stable state. The parameters from the
optimization with highest likelihood score were selected to represent
the best optimization. The standard deviations of each set of parameters
was estimated using 100 nonparametric bootstraps of the 2d-sfs generated
with \selectlanguage{greek}δ\selectlanguage{english}a\selectlanguage{greek}δ\selectlanguage{english}i, and then the Godambe Information Matrix (GIM) was used to
analyze the uncertainty of the parameters based on the 100 bootstrap
2d-sfs (Coffman, Hsieh, Gravel, \& Gutenkunst, 2016).

The TMRCA estimations of paired populations were used as dependent
variable, with grouping status (from the same or different drainage,
cluster, or morphotype) used as a categorical independent variable and
geographical distance as a covariate in an analysis of covariance
(ANCOVA). For each ANCOVA regression model, AIC and BIC were calculated,
and Vuong's test was used to compare those models using the R package
nonnest2 version 0.5-5 (Vuong, 1989; Merkle \& You, 2018).

\emph{Testing effects of post glacial population structure}

Compared to the estimation of divergence time, the fixation
index\emph{F} \textsubscript{ST} is more sensitive to gene flow among
populations, where paired populations experiencing gene flow are
expected to have lower \emph{F} \textsubscript{ST} values (in a 0 to 1
scale). Therefore, if gene flow occurred across drainage basins after
populations recolonized high elevation areas (i.e. across mountain
ridges), we expect to see little to no discrepancy in comparisons of
populations from the same or different watersheds. For the same set of
population pairs of the best model in the TMRCA analysis, we calculated
the means of Weir \& Cockerham's \emph{F} \textsubscript{ST} in VCFtools
version 0.1.17 (Weir \& Cockerham 1984, Danecek et al., 2011) using
the\emph{semi-filtered dataset} . Similar to the TMRCA, the ANCOVA was
used to test whether population pairs from different drainages show
higher\emph{F} \textsubscript{ST}, with geographic distance between
paired populations as a covariate.

\textbf{Results}

\emph{Mating trials}

Mating trials demonstrated that each morphotype within the \emph{N.
ingens} complex could successfully copulate with its own morphotype and
with other morphotypes (\textbf{Table 3} ). While we observed a
reduction in successful mating in experiments involving \emph{N. ingens}
, this also occurred when \emph{N. ingens} males were paired with
\emph{N. ingens} females, suggesting that it may have been linked to
reproductive status (\emph{e.g} . reproductive diapause). Successful
copulation occurred across all combinations of morphotypes, although not
always for both male and female directed crosses. We note, however, that
sample sizes are small.

\emph{Mitochondrial phylogeny}

The COI phylogenetic tree has two major lineages, a northern and
southern clade with divergence time of 0.52 million years ago, and both
clades contain individuals from the intermediate morphotype
(\textbf{Figure 2} ). There is clear geographical structure, with
populations from the same drainage belonging to the same clades and even
the same haplotype. For example, there is a distinct subclade in the
northern clade that includes samples from the Tuolumne and Merced
drainages, and excludes other haplotypes found in the northern clade.
Similarly, several clades in the southern lineage consist of samples
only from the Kern or Kern/Kaweah drainage. The TCS network shows that
the two distinct haplogroups are separated by eight substitutions
(\textbf{Figure S2} ). Haplotypes assigned to the northern clade are
widely distributed, although many are unique and occur locally. These
haplotypes are found in the Tuolumne (\emph{N. riversi} ), San Joaquin
(intermediate morphotype), Kings, and Kern (\emph{N. ingens} )
drainages, including 13 out of the 21 populations of intermediate and
\emph{N. ingens} sample sites, as far as the southernmost population
(site27). On the other hand, the southern clade consists only of samples
from San Joaquin (intermediate morphotype), Kings, Kern, and Kaweah
(\emph{N. ingens} ), and no southern clade haplotypes are observed in
northern populations of \emph{N. riversi} . (\textbf{Table S1} ).

\emph{Nuclear Phylogeny}

The species tree inferred using SNAPP, based on individuals from eight
sample sites, also shows two major clades. However, a northern clade
includes only samples from the Tuolumne, Merced, and San Joaquin
drainages, whereas a southern clade only includes samples from Kings,
Kaweah, and Kern drainages (\textbf{Figure 3} ). The species tree
topology has strong posterior branch support for all nodes, and suggests
that morphologically intermediate populations from the San Joaquin
watershed are genetically closer to the northern \emph{N.
riversi}individuals. In contrast, the concatenated Maximum likelihood
analysis using IQtree does not show two major clades, but instead
represents the evolutionary divergence as hierarchically nested branches
of populations that reflect geographical proximity (\textbf{Figure 4} ).
Bootstrap support is high for populations, though weak for a few
branches joining different populations. It is also worthwhile to note
that the tree topology largely supports a clear geographical subdivision
by drainage system (\textbf{Figure 4} ).

\emph{Population Structure and Ancestry Coefficients}

The clustering of individuals into ancestral populations was estimated
with sNMF using cluster numbers from 1 to 20. The optimal number of
clusters was unclear as the cross-entropy continuously decreases up to
K=11, although the relative change in cross-entropy values is smaller
after K=6 (\textbf{Figure S3} ). In the K=3 model, clusters are
consistent with morphotypes (\emph{N. riversi} , intermediate,
and\emph{N. ingens} ), but with some populations showing some admixture
among these major lineages. These admixed populations are located in a
geographical transition zone where separate lineages could have
experienced historical gene flow (\textbf{Figure S1} ). In the K=4
model, the southern cluster splits into two clusters, revealing a
distinct ancestral subpopulation located in the Kings drainage. In the
K=5 model, the southernmost cluster further subdivides into west
(Kaweah) and east (Kern) ancestral populations. In the K=6 model, the
northern sample sites (sites 3 and 4) split from other northern
populations and thus subdivide \emph{N. riversi} (sites 1 to 6) into two
clusters. From K=7 to K=11, new clusters mostly emerge around the
contact zones of major drainages (\textbf{Figure S1} ).

\emph{Testing spatial clustering by genetic cline}

Comparisons between spatial and non-spatial models using conStruct
reveal that the spatial model is always favored (\textbf{Figure} 5), as
the accuracies of the spatial model are significantly higher than that
of non-spatial model in K=1 to K=8. In addition, the analysis of layer
contributions for K=1 to K=10 shows that additional layers contribute
very little at K=5 and become negligible at K=7 (\textbf{FigureS4} ).
Viewed in combination with the results from sNMF, K=4 to K=6 are
preferred models.

\emph{Estimate of demographic history}

In the analyses of demographic change, the stairway plots of the eight
populations show similar patterns, with population size declining
dramatically from around 100 thousand to a few thousands individuals
within the last 30 thousand years (\textbf{Figure 6} ). According to the
stairway plots, the sample sites 1 and 3 (\emph{N. riversi} ) have
similar patterns of population decline starting around 25 thousand years
ago, but the ancestral population was larger in site 3, a site located
in the contact zone of the Merced and San Joaquin drainage systems
(\textbf{Figure 6A \& B} ). The sites 9, 10, and 12 represent
intermediate morphotypes located within San Joaquin drainage, and they
share similar patterns of decline starting around 10 thousand years ago
(\textbf{Figure 6C, D, \&E} ). The three sites of \emph{N. ingens} ,
from the Kings (site 17), Kaweah (site 24), and Kern (site 27), shows
declines ranging from 20 thousand to two thousand years ago
(\textbf{Figure 6F, G, \&H} ).

\emph{Testing glacial refugia and post glacial recolonization}

The TMRCA estimation of 96 pairs of populations were used to conduct
ANCOVA for models of K=3 to K=6. In all tested models, the interaction
terms are not significant and thus the type II sum of squares was used
in ANCOVA tests. The K=6 drainage refugia model includes 30 pairs of
populations from the same drainage system and 66 pairs of populations
from different drainage systems. This model has the lowest AIC and BIC,
and explains the data significantly better than either K=3 or K=4
population structure models based on Vuong's test (\textbf{Table 1} ).
However, we cannot distinguish the K=5 and K=6 population structure
model, as they possess the same set of paired populations due to the
lack of multiple population samples form the Merced drainage.

In the best fit K=6 model, the result of ANCOVA shows that TMRCA
estimates are significantly older (\textgreater{} 20,000 years ago) for
pairs of populations from different drainages than from the same
drainage (\textasciitilde{} 10,000 years ago; p-value \textless{}
0.001), regardless of whether geographical distance is included as a
covariant or not (\textbf{Figure 7A} ). However, there are a few of
population pairs from the same drainage that do have older TMRCA,
notably sites in the Kings or Kern drainages. Similarly, there are a few
outliers with recent TMRCA when pairs of different drainages are
compared, including sites from the San Joaquin and Kings drainages. In
contrast, there is no significant difference between the two types of
population pairs when measuring Weir and Cockerham's\emph{F}
\textsubscript{ST} (p-value \textgreater{} 0.05). Geographical distance
has a high \emph{R\textsuperscript{2}} value (0.35 and 0.55) when
included as a covariate, implying genetic variation is consistent with
isolation by distance (\textbf{Figure 7} B).

\textbf{Discussion}

The main goal of this study was to test competing hypotheses for the
location of glacial refugia in an alpine beetle species, including a
hypothesis focusing on the role of drainage basins as refugia. By
focusing on analyses of population structure and TMRCA estimates from
paired populations of the cold-specialized, riparian beetles in the
alpine \emph{N. ingens}complex, we found strong evidence that population
structure and patterns of lineage divergence are correlated with
drainage basins in the Sierra Nevada Mountains, California. This
supports the hypothesis that drainage basins, and by extension
microhabitat preferences of a species, allow for persistence in multiple
refugia across the mountain range. We also provide evidence that
post-glacial gene flow can readily obscure the signature of glacial
refugia, if not carefully dissected, and that episodes of glaciation
enhanced lineage diversity in this alpine species complex through both
isolation and hybridization.

\emph{Disentangling sources of population structure}

The population structure of a species reflects both historical
divergence processes and contemporary gene flow among subpopulations,
and it requires careful analysis to disentangle these evolutionary
processes (Nielsen, \& Slatkin, 2007; Sch\selectlanguage{ngerman}önswetter \& Schneeweiss 2019).
Summary statistics, such as Wright's \emph{F} \textsubscript{ST}, are
useful for describing genetic drift and gene flow under equilibrium
models, but provide little information about deeper evolutionary events
that shape patterns of genetic variation. Clustering based algorithms,
which infer individual ancestry components according to a model, are
also challenging to interpret when population divergence is continuous
in space (François \& Durand, 2010; but see Bradburd et al. 2018).
Coalescent approaches leverage the distribution of gene genealogies to
separately estimate the population parameters, including migration rate
and divergence time, and can be used to distinguish models of migration
from divergence (Nielsen \& Wakeley, 2001). However, coalescent
approaches are often difficult to implement when evolutionary history is
spatially complex, as the number of free parameters grows rapidly
(Beerli, 2009).

In this study, we leveraged aspect of each approach to disentangle the
evolutionary history of the \emph{N. ingens} complex. First, we used the
clustering algorithm sNMF to identify population structure, though we
found that the log likelihood increased with an increasing number of
clusters. A large shift in K occurred as different morphotypes were
distinguished genetically (K=3), but models K=4 to K=6 had improved log
likelihood scores and revealed further population subdivision reflecting
major drainage systems in the Sierra Nevada. More subtle patterns of
admixture among drainages continued to emerge as the number of clusters
increases (\textbf{Figure S3} ), but this appears to be the product of
gene flow and uneven sampling (continuous population structure). Indeed,
by explicitly comparing spatial and non-spatial cluster models using
conStruct, we found that the spatial models are generally more accurate
than the non-spatial models (\textbf{Figure 5} ), indicating the
presence of isolation by distance among sample sites. Analysis of
pairwise\emph{F} \textsubscript{ST} values among populations also
strongly supports gradual genetic divergence with geographical distances
(\emph{R\textsuperscript{2}} = 0.35; \textbf{Figure 7} ).

We further assessed patterns of population divergence using a model
fitting framework based on coalescence theory (Rougemont \& Bernatchez,
2018). We compared ANCOVA regressions on the paired populations' TMRCA
for different levels of population structure. For TMRCA, K=6 and K=5
show the best AIC, BIC, and results from the Vuong's test. Although K=6
and K=5 are indistinguishable owing to the lack of multiple population
samples from the Merced drainage, both support the model of drainage
associated glacial refugia (\textbf{Table 1} ).

\emph{Differentiating drainage based multiple refugia from nunatak
refugia}

Although our results suggest multiple ancestral populations and high
levels of genetic diversity in some high elevation populations, which
could fit the predictions of glacial survival within the ice sheet
(nunatak survival hypothesis), the pattern of several isolated glacial
refugia of \emph{N. ingens} complex is fundamentally different from
nunataks. Conventionally, nunatak refugia refer to the small, ice-free
geographical areas within the glacial ice sheets. In such habitats,
populations would be restricted to small size and exposed to extremely
challenging environmental conditions (cold and windy, and thus xeric),
before expanding in postglacial periods (Dahl 1987; Holderegger \&
Thiel\selectlanguage{english}-Egenter 2009). In the case of the \emph{N. ingens} complex, high
elevation sites with elevated genetic diversity reflect recent admixture
among ancestral population clusters (conStruct results, see\textbf{Table
1} ). Ancestral populations are not associated with ice-free areas in
the north and central part of the Sierra Nevada (\textbf{Figure 1} ).
Similar contact zones are found in a large set of alpine plants from the
European Alps (Thiel\selectlanguage{english}-Egenter et al., 2011), and a recent synthesis of
alpine studies suggests this might be a common phenomenon across long
linear mountain ranges (Wallis et al., 2016). Many alpine or sub-alpine
species that have patterns consistent with multiple glacial refugia are
certainly not related to persistence in nunataks (Kubow, Robinson,
Shama, \& Jokela, 2010; Homburg et al., 2013), suggesting that the deep
divergence of population structure or multiple ancestral populations
might not be sufficient to be directly linked to nunatak refugia.

The other common character of populations recolonized from the nunataks
is the vestige of a bottleneck during glacial periods (Segarra\selectlanguage{english}-Moragues,
Palop\selectlanguage{english}-Esteban, Gonzalez\selectlanguage{english}-Candelas, \& Catalan, 2007; Kosi\selectlanguage{polish}ń\selectlanguage{english}ski et al.,
2019). Our data suggest the opposite has happened, as the stairway plots
show drastic population decline initiating at the last glacial maximum
and accelerating into the present (\textbf{Figure 7} ). Similarly,
estimates of ancestral population size (reference population,
N\textsubscript{ref}) in \selectlanguage{greek}δ\selectlanguage{english}a\selectlanguage{greek}δ\selectlanguage{english}i were several fold larger than the
descendant populations (\textbf{Table S2} ). Although other studies
report glacial bottlenecks or postglacial population increases
(Bidegaray-Batista et al., 2016; Sim, Hall, Jex, Hegel, \& Coltman,
2016; Starcov\selectlanguage{ngerman}á, Vohralík, Kry\selectlanguage{polish}š\selectlanguage{english}tufek, Bolf\selectlanguage{ngerman}íková, \& Hulva, 2016; Weng,
Yang, \& Yeh, 2016; Huang et al., 2017), it seems logical that alpine
species would undergo founder events as they recolonize high elevation
areas following glacier retreat (Schoville and Roderick, 2009).

We acknowledge that the continuous spatial structure of population could
possibly lead to a false signature of recent population decline for
stairway plot analysis and possibly other sfs-based analyses including
\selectlanguage{greek}δ\selectlanguage{english}a\selectlanguage{greek}δ\selectlanguage{english}i (Battey, Ralph, \& Kern, 2020). However, there are several reasons
we interpret the larger ancestral population sizes estimated by the
stairway plot and \selectlanguage{greek}δ\selectlanguage{english}a\selectlanguage{greek}δ\selectlanguage{english}i to represent a true demographic decline. First,
the census population size of populations within the \emph{N.
ingens}complex is mostly very small and populations are isolated from
each other by complex topographic terrain (\emph{e.g} . suitable habitat
occurs within the first few hundred meters from the headwaters of
streams or small seep areas within glacier cirques). Thus, estimates of
a contemporary effective population size numbering in the hundreds or
several thousand is realistic in our opinion. Second, all eight sample
sites show a consistent trend of population decline, even though they
differ in their neighborhood size and geographical position in the
geographical cline, which are major factors that bias the stairway plot
estimates in continuously structured populations (Battey et al., 2020).
Lastly, the cool and wet glacial period preceding the glacial maximum
would likely be optimal for a cold-specialized, riparian beetle, and
thus explains the larger population size at this time. A few other
studies have shown population expansion during glaciation followed by
postglacial population declines in alpine and subalpine species
(Galbreath, Hafner, \& Zamudio, 2009; Huang et al., 2016; Zaman, Hubert,
\& Schoville, 2019).

\emph{Divergence, admixture and speciation}

The role of climate fluctuations in causing secondary contact of
recently diverged species is well known (Hewitt 1999), with the outcomes
varying from free admixture, to directional gene flow, to complete
isolation or even hybrid speciation. Several cases demonstrating the
range of these outcomes are known from alpine taxa within the Sierra
Nevada Mountains (Gompert, Fordyce, Forister, Shapiro, \& Nice, 2006;
Rovito 2010; Rubidge, Patton, \& Moritz, 2014). The degree and
geographical scale of admixture likely varies across taxa due to a
number of factors, including the dispersal ability of species, habitat
suitability and connectivity, as well as the reproductive compatibility
of lineages and fitness of hybrids (Barton \& Hewitt 1985; Dufresnes,
Berroneau, Dubey, Litvinchuk, \& Perrin, 2020).

In this study, we found that an intermediate morphotype in the \emph{N.
ingens} complex (Schoville et al. 2012) extends to populations
throughout the Central Sierra Nevada, encompassing the entire San
Joaquin drainage and extending south to contact zones in the Kings River
watershed (sites 14 and 15; \textbf{Figure 1} ). The distribution of the
morphological phenotype is different from the observed genetic patterns,
where mitochondrial haplotype variation suggests introgression
from\emph{N. riversi} into nearly the entire southern range of \emph{N.
ingens} (\textbf{Figure 2 \& S2} ), and genome-wide nuclear markers
suggest that the intermediate morphotype is an independent lineage (site
9, 10 and 12; \textbf{Figure 3} ). For example, in the sNMF analysis
(from K=3 to K=11), most individuals representing the intermediate
morphotype do not show extensive ancestry from either \emph{N.
riversi}or \emph{N. ingens} (\textbf{Figure S1} ). Additionally, TMRCA
estimates between populations of the intermediate morphotype and either
species mostly predate the last glacial maximum. Despite this genetic
divergence, we found no evidence of premating isolation among the three
lineages (\textbf{Table 3} ), although our sample sizes were low.

Based on these results, we infer that 1) an early glacial cycle caused
deep divergence to form between northern and southern populations, which
formed \emph{N. ingens} and \emph{N. riversi} . In a subsequent glacial
cycle, 2) a contact zone formed near the San Joaquin drainage and led to
the formation of the intermediate morphotype. This lineage diversity was
maintained during the last glacial maximim and 3) allowed for the
formation of more complex population structure as each major lineage
recolonized high elevation habitat. It is clear that limited gene flow
occurs in parts of the alpine range, leading to the formation of a
genetic cline in genetic divergence measures such as\emph{F}
\textsubscript{ST}. The general pattern we observe is consistent with
models of diversification linked to the orbitally forced range dynamic
(ORD) model (Dynesius \& Jansson, 2000), and the more recently coined
'mix-isolation-mix' model proposed by He et al. (2019).

\emph{Implications for alpine species conservation}

Our finding of a drastic population decline since the last glaciation
retreat, coupled with field observations of extremely limited habitat
and small population sizes, suggests that this species complex may be
vulnerable to localized extinction. For many \emph{Nebria} species in
North America, lower elevational distributions have shifted upward in
the past 40 years (Kavanaugh 2015). If population range shifts are
associated with warmer temperatures in the Holocene, it is very likely
that populations will continue to shrink with ongoing global warming.
This could be exacerbated by the loss of Sierra Nevada glaciers and a
reduction in winter snowpack (Basagic \& Fountain, 2011; Belmecheri,
Babst, Wahl, Stahle, \& Trouet, 2016). Although \emph{N. ingens} complex
is not an economically significant animal, it is difficult to predict
the consequence of a general predator disappearing from a comparatively
depauperate alpine ecosystem.

\textbf{Acknowledgements}

We are grateful to the National Parks of the United States including
Yosemite, Sequoia, and King Canyon National Parks for providing
permission to conduct research in the parks. Thanks also to Ministry of
Education Republic of China (Taiwan) for offering a scholarship to YMW.
This research is also supported by the Valentine Eastern Sierra Reserve
(UC Santa Barbara Natural Reserve System) through VESR Graduate Student
Funding to YMW, and a UW Research Competition Award to SDS. We thank
Benton Veire and Matthew Medeiros for helping us to conduct the field
work, and Ph.D. committees of YMW for providing suggestions to this
study. The authors have no conflicts of interests to declare.

\textbf{References}

Alvarez, N., Thiel-Egenter, C., Tribsch, A., Holderegger, R., Manel, S.,
Sch\selectlanguage{ngerman}önswetter, P., Taberlet, P., Brodbeck, S., Gaudeul, M., Gielly, L.,
Küpfer, P., Mansion, G., Negrini, R., Paun, O., Pellecchia, M., Rioux,
D., Schüpfer, F., Van Loo, M., Winkler, M., Gugerli, F., \& IntraBioDiv
Consortium (2009). History or ecology? Substrate type as a major driver
of spatial genetic structure in Alpine plants.~\emph{Ecology Letters}
,~\textbf{12(7)} , 632- 640.

Andújar, C., Serrano, J. \& Gómez-Zurita, J. (2012). Winding up the
molecular clock in the genus \emph{Carabus} (Coleoptera: Carabidae):
assessment of methodological decisions on rate and node age
estimation.\emph{BMC Evolutionary Biology} , \textbf{12} , 40.

Bankevich, A., S. Nurk, Antipov, D., Gurevich, A. A., Dvorkin, M.,
Kulikov, A. S., Lesin, V. M., Nikolenko, S. I., Pham, S., Prjibelski, A.
D., Pyshkin, A. V., Sirotkin, A. V., Vyahhi, N., Tesler, G., Alekseyev,
M. A., \& Pevzner, P. A. (2012). SPAdes: a new genome assembly algorithm
and its applications to single-cell sequencing. \emph{Journal of
Computational Biology} , \textbf{19} , 455- 477.

Barton, N. H., \& Hewitt, G. M. (1985). Analysis of hybrid
zones.\emph{Annual Review of Ecology and Systematics} , \textbf{16} ,
113- 148.

Basagic, H. J., \& Fountain, A. G. (2011). Quantifying 20th century
glacier change in the Sierra Nevada, California. \emph{Arctic,
Antarctic, and Alpine Research} , \textbf{43(3)} , 317- 330.

Battey, C. J., Ralph, P. L., \& Kern, A. D. (2020). Space is the place:
effects of continuous spatial structure on analysis of population
genetic data. \emph{Genetics} , \textbf{215} , 193- 214.

Beerli, P. (2015). How to use MIGRATE or why are Markov chain Monte
Carlo programs difficult to use. \emph{Population Genetics for Animal
Conservation} , \textbf{17} , 42- 79.

Belmecheri, S., Babst, F., Wahl, E. R., Stahle, D. W., \& Trouet, V.
(2016). Multi-century evaluation of Sierra Nevada snowpack. \emph{Nature
Climate Change} , \textbf{6(1)} , 2- 3.

Bidegaray-Batista, L., Sánchez-Gracia, A., Santulli, G., Maiorano, L.,
Guisan, A., Vogler, A. P., \& Arnedo, M. A. (2016). Imprints of multiple
glacial refugia in the Pyrenees revealed by phylogeography and
palaeodistribution modelling of an endemic spider. \emph{Molecular
Ecology} , \textbf{25} , 2046- 2064.

Bohonak, A. J. (1999). Dispersal, gene flow, and population
structure.\emph{The Quarterly Review of Biology} , \textbf{74} , 21- 45.

Bouckaert, R. R. (2010). DensiTree: making sense of sets of phylogenetic
trees. \emph{Bioinformatics} , \textbf{26} , 1372- 1373.

Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C-H., Xie, D.,
Suchard, MA., Rambaut, A., \& Drummond, A. J. (2014). BEAST 2: A
Software Platform for Bayesian Evolutionary Analysis. \emph{PLoS
Computational Biology} , \textbf{10(4)} , e1003537.
doi:10.1371/journal.pcbi.1003537

Bradburd, G., Coop, G., \& Ralph, P. (2018). Inferring continuous and
discrete population genetic structure across space. \emph{Genetics}
,\textbf{210} , 33- 52.

Bryant, D., Bouckaert, R., Felsenstein, J., Rosenberg, N. A., \&
RoyChoudhury, A. (2012). Inferring species trees directly from biallelic
genetic markers: bypassing gene trees in a full coalescent
analysis.\emph{Molecular Biology and Evolution} , \textbf{29} , 1917-
1932.

Carling, M. D., Lovette, I. J., \& Brumfield, R. T. (2010). Historical
divergence and gene flow: coalescent analyses of mitochondrial,
autosomal and sex-linked loci in \emph{Passerina buntings}
.\emph{Evolution} , \textbf{64} : 1762- 1772.

Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., \& Cresko, W. A.
(2013). Stacks: an analysis tool set for population
genomics.\emph{Molecular Ecology} , \textbf{22} , 3124- 3140.

Clement, M., Snell, Q., Walker, P., Posada, D., \& Crandall, K. (2002).
TCS: estimating gene genealogies. Brigham Young University, Provo, UT.

Coffman, A. J., Hsieh, P. H., Gravel, S., \& Gutenkunst, R. N. (2016).
Computationally Efficient Composite Likelihood Statistics for
Demographic Inference. \emph{Molecular Biology and Evolution}
,\textbf{33(2)} , 591- 593. https://doi.org/10.1093/molbev/msv255

Dahl, E. (1987). The nunatak theory reconsidered. \emph{Ecological
Bulletins} , 77-94.

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E.,
DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S.
T., McVean, G., Durbin, R., \& 1000 Genomes Project Analysis Group.
(2011). The variant call format and VCFtools. \emph{Bioinformatics}
,\textbf{27} , 2156- 2158.

DeChaine, E. G., \& Martin, A. P. (2006). Using coalescent simulations
to test the impact of Quaternary climate cycles on divergence in an
alpine plant\selectlanguage{english}-insect association. \emph{Evolution} , \textbf{60(5)} ,
1004- 1013.

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha,
S., Batut, P., Chaisson, M., \& Gingeras, T. R. (2013). STAR: ultrafast
universal RNA-seq aligner. \emph{Bioinformatics} , \textbf{29} , 15- 21.

Dufresnes, C., Berroneau, M., Dubey, S., Litvinchuk, S. N., \& Perrin,
N. (2020). The effect of phylogeographic history on species boundaries:
a comparative framework in \emph{Hyla} tree frogs. \emph{Scientific
Reports} , \textbf{10} , 1- 12.

Dynesius, M., \& Jansson, R. (2000). Evolutionary consequences of
changes in species' geographical distributions driven by Milankovitch
climate oscillations. \emph{Proceedings of the National Academy of
Sciences} , \textbf{97(16)} , 9115- 9120.

Elshire, R.J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K.,
Buckler, E. S., \& Mitchell, S. E. (2011). A robust, simple
genotyping-by-sequencing (GBS) approach for high diversity
species.\emph{PLoS ONE} , \textbf{6} , e19379.

Francois, O. \& Durand, E. (2010). Spatially explicit Bayesian
clustering models in population genetics. \emph{Molecular Ecology
Resources} , \textbf{10(5)} , 773- 784

Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., \& Francois, O.
(2014). Fast and efficient estimation of individual ancestry
coefficients. \emph{Genetics} , 196, 973-983.

Galbreath, K. E., Hafner, D. J., \& Zamudio, K. R. (2009). When cold is
better: climate-driven elevation shifts yield complex patterns of
diversification and demography in an alpine specialist (American
pika,\emph{Ochotona princeps} ). \emph{Evolution} , \textbf{63} , 2848-
2863.

Gereben, B. A. (1995). Co-occurrence and microhabitat distribution of
six \emph{Nebria} species (Coleoptera: Carabidae) in an alpine glacier
retreat zone in the Alps, Austria. \emph{Arctic and Alpine Research}
,\textbf{27(4)} , 371- 379.

Gompert, Z., Fordyce, J. A., Forister, M. L., Shapiro, A. M., \& Nice,
C. C. (2006). Homoploid hybrid speciation in an extreme
habitat.\emph{Science} , \textbf{314(5807)} , 1923- 1925.

Gutenkunst, R. N., Hernandez, R. D., Williamson, S. H., \& Bustamante,
C. D. (2009). Inferring the joint demographic history of multiple
populations from multidimensional SNP frequency data. \emph{PLoS
Genetics} , 5, e1000695.

He, Z., Li, X., Yang, M., Wang, X., Zhong, C., Duke, N. C., Wu, C. I.,
\& Shi, S. (2019). Speciation with gene flow via cycles of isolation and
migration: insights from multiple mangrove taxa. \emph{National Science
Review} , \textbf{6(2)} , 275- 288.

Hewitt, G. M. (1999). Post-glacial re-colonization of European
biota.\emph{Biological Journal of the Linnean Society} , \textbf{68} ,
87- 112.

Holderegger, R., \& Thiel\selectlanguage{english}-Egenter, C. (2009). A discussion of different
types of glacial refugia used in mountain biogeography and
phylogeography\emph{. Journal of Biogeography} , \textbf{36} , 476- 480.

Homburg, K., Drees, C., Gossner, M. M., Rakosy, L., Vrezec, A., and
Assmann, T. (2013). Multiple glacial refugia of the low-dispersal ground
beetle \emph{Carabus irregularis} : molecular data support predictions
of species distribution models. \emph{PLoS One} , \textbf{8} , e61185.

Huang, C. C., Hsu, T. W., Wang, H. V., Liu, Z. H., Chen, Y. Y., Chiu, C.
T., Huang, C. L., Hung, K. H., \& Chiang, T. Y.. (2016). Multilocus
analyses reveal postglacial demographic shrinkage of \emph{Juniperus
morrisonicola} (Cupressaceae), a dominant alpine species in
Taiwan.\emph{PLoS One} , \textbf{11} , e0161713.

Huang, Z. S.,. Yu, F. L, Gong, H. S., Song, Y. L., Zeng, Z. G., \&
Zhang, Q. (2017). Phylogeographical structure and demographic expansion
in the endemic alpine stream salamander (Hynobiidae:\emph{Batrachuperus}
) of the Qinling Mountains. \emph{Scientific Reports} , \textbf{7} ,
1871.

Jansson, R., \& Dynesius, M. (2002). The fate of clades in a world of
recurrent climatic change: Milankovitch oscillations and
evolution.\emph{Annual Review of Ecology and Systematics} ,
\textbf{33(1)} , 741-777.

Kavanaugh, D. H. (1985). On wing atrophy in carabid beetles (Coleoptera:
Carabidae), with special reference to Nearctic \emph{Nebria} . In G. E.
Ball (Ed.), \emph{Taxonomy, Phylogeny and zoogeography of beetles and
ants.} \emph{A Volume Dedicated to the Memory of Philip Jackson
Darlington, Jr.(1904-1983). Series Entomologica} , \textbf{33} ,
408-431.

Kavanaugh, D. H. (2015). A new species of \emph{Nebria} Latreille
(Insecta: Coleoptera: Carabidae: Nebriini) from the Spring Mountains of
Southern Nevada. \emph{Proceedings of the California Academy of
Sciences} , \textbf{62} , 441--450.

Keightley, P. D., Ness, R. W., Halligan, D. L., \& Haddrill, P. R.
(2014). Estimation of the spontaneous mutation rate per nucleotide site
in a \emph{Drosophila melanogaster} full-sib family. \emph{Genetics}
,\textbf{196} , 313-320.

Keightley, P. D., Pinharanda, A., Ness, R. W., Simpson, F.,
Dasmahapatra, K. K., Mallet, J., Davey, J. W., \& Jiggins, C. D. (2015).
Estimation of the spontaneous mutation rate in \emph{Heliconius
melpomene} . \emph{Molecular Biology and Evolution} , \textbf{32} ,
239-243.

Kikvidze Z, Brooker RW, Butterfield BJ, Callaway RM, Cavieres LA, Cook
BJ, Lortie CJ, Michalet R, Pugnaire FI, Xiao S, Anthelme F, Bjork RG,
Cranston BH, Gavilan RG, Kanka R, Lingua E, Maalouf JP, Noroozi J,
Parajuli R, Phoenix GK, Reid A, Ridenour WM, Rixen C, Schob C. 2015. The
effects of foundation species on community assembly: A global study on
alpine cushion plant communities. \emph{Ecology} , \textbf{96} , 2064--
2069.

Kleckova, I., Konvicka, M., \& Klecka, J. (2014). Thermoregulation and
microhabitat use in mountain butterflies of the genus \emph{Erebia} :
importance of fine-scale habitat heterogeneity. \emph{Journal of Thermal
Biology} , \textbf{41} , 50- 58.

Knowles, L. L. (2001). Did the Pleistocene glaciations promote
divergence? Tests of explicit refugial models in montane
grasshopprers.\emph{Molecular Ecology} , \textbf{10(3)} , 691- 701.

Korneliussen, T. S., Albrechtsen, A., \& Nielsen, R. (2014). ANGSD:
Analysis of Next Generation Sequencing Data. \emph{BMC Bioinformatics}
,\textbf{15} , 356.

Kosi\selectlanguage{polish}ń\selectlanguage{english}ski, P., S\selectlanguage{polish}ę\selectlanguage{english}kiewicz, K., Walas, \selectlanguage{polish}Ł., \selectlanguage{english}Boraty\selectlanguage{polish}ń\selectlanguage{english}ski, A., \& Dering, M.
(2019). Spatial genetic structure of the endemic alpine plant\emph{Salix
serpillifolia} : genetic swamping on nunataks due to secondary
colonization? \emph{Alpine Botany} , \textbf{129} , 107- 121.

Karney, C. F. F. (2013). Algorithms for geodesics. \emph{Journal of
Geodesy} , \textbf{87} , 43- 55.

Kubow, K. B., Robinson, C. T., Shama, L. N., \& Jokela, J. (2010).
Spatial scaling in the phylogeography of an alpine caddisfly, Allogamus
uncatus, within the central European Alps. \emph{Journal of the North
American Benthological Society} , \textbf{29} , 1089- 1099.

Kumar, S., Stecher, G., \& Tamura, K. (2016). MEGA7: Molecular
Evolutionary Genetics Analysis Version 7.0 for Bigger
Datasets.\emph{Molecular Biology and Evolution} , \textbf{33} , 1870-
1874.

Leigh, J. W., \& Bryant, D. (2015). PopArt: full-feature software for
haplotype network construction. \emph{Methods in Ecology and Evolution}
,\textbf{6 (9)} , 1110- 1116.

Liu, X. (2020). Human prehistoric demography revealed by the polymorphic
pattern of CpG transitions. \emph{Molecular Biology and Evolution} ,
msaa112

Liu, X., \& Fu, Y. X. (2015). Corrigendum: Exploring population size
changes using SNP frequency spectra. \emph{Nature Genetics} ,\textbf{47}
, 1099.

Lohse, K., Nicholls, J. A., \& Stone, G. N. (2011). Inferring the
colonization of a mountain range---refugia vs. nunatak survival in high
alpine ground beetles. \emph{Molecular Ecology} , \textbf{20(2)} , 394-
408.

Merkle, E., \& You, D. (2018). nonnest2: tests of non-nested models.
Retrieved 04 August, 2020, from
https://cran.r-project.org/web/packages/nonnest2

Massatti, R., \& Knowles, L. L. (2014). Microhabitat differences impact
phylogeographic concordance of codistributed species: Genomic evidence
in montane sedges (\emph{Carex} L.) from the Rocky
Mountains.\emph{Evolution} , \textbf{68} , 2833- 2846.

Nguyen, L.T., Schmidt, H. A., Von Haeseler, A., \& Minh, B. Q. (2015).
IQ-TREE: a fast and effective stochastic algorithm for estimating
maximum-likelihood phylogenies. \emph{Molecular Biology and Evolution}
,\textbf{32} , 268- 274.

Nielsen, R., \& Wakeley, J. (2001). Distinguishing migration from
isolation: a Markov chain Monte Carlo approach. \emph{Genetics}
,\textbf{158(2)} , 885- 896.

Nielsen, R., \& Slatkin, M. (2007). Likelihood analysis of ongoing gene
flow and historical association. \emph{Evolution} , \textbf{54(1)} , 44-
50.

Oppold, A. M., \& Pfenninger, M. (2017). Direct estimation of the
spontaneous mutation rate by short-term mutation accumulation lines
in.\emph{Evolution Letters} , \textbf{1} , 86- 92.

Pan, D., H\selectlanguage{ngerman}ülber, K., Willner, W., \& Schneeweiss, G. M. (2020). An
explicit test of Pleistocene survival in peripheral versus nunatak
refugia in two high mountain plant species. \emph{Molecular Ecology}
,\textbf{29} , 172- 183.

Paun, O., Schönswetter, P., Winkler, M., Consortium, I. \& Tribsch, A.
(2008), Historical divergence vs. contemporary gene flow: evolutionary
history of the calcicole \emph{Ranunculus alpestris} group
(Ranunculaceae) in the European Alps and the Carpathians.\emph{Molecular
Ecology} , \textbf{17} : 4263- 4275.

Portik, D. M., Leaché, A. D., Rivera, D., Barej, M. F., Burger, M.,
Hirschfeld, M., Rödel, M. O., Blackburn, D. C., \& Fujita, M. K. (2017).
Evaluating mechanisms of diversification in a Guineo\selectlanguage{english}-Congolian tropical
forest frog using demographic model selection. \emph{Molecular Ecology}
,\textbf{26} , 5245- 5263.

Qu, Y., Zhang, R., Quan, Q., Song, G., Li, S.H. \& Lei, F. (2012).
Incomplete lineage sorting or secondary admixture: disentangling
historical divergence from recent gene flow in the Vinous\selectlanguage{english}-throated
parrotbill (\emph{Paradoxornis webbianus} ). \emph{Molecular Ecology}
,\textbf{21} : 6117- 6133.

Rambaut, A. (2012). FigTree v1. 4. Available from
https://beast.community/figtree

Rambaut, A., Drummond, A. J., Xie, D., Baele, G., \& Suchard, M. A.
(2018). Posterior summarization in Bayesian phylogenetics using Tracer
1.7. \emph{Systematic Biology} , \textbf{67} , 901.

Rood, D. H.,~Burbank, D. W., \&~Finkel, R. C.~(2011).~Chronology of
glaciations in the Sierra Nevada, California,
from~\textsuperscript{10}Be surface exposure dating.~\emph{Quaternary
Science Reviews} ,~\textbf{30} :~646- 661.

Rougemont, Q., \& Bernatchez, L. (2018). The demographic history of
Atlantic salmon (\emph{Salmo salar} ) across its distribution range
reconstructed from approximate Bayesian computations. \emph{Evolution}
,\textbf{72} : 1261- 1277

Rovito, S. M. (2010). Lineage divergence and speciation in the Web\selectlanguage{english}-toed
Salamanders (Plethodontidae: \emph{Hydromantes} ) of the Sierra Nevada,
California. \emph{Molecular Ecology} , \textbf{19(20)} , 4554- 4571.

Rovito, S. M., \& Schoville, S. D. (2017). Testing models of refugial
isolation, colonization and population connectivity in two species of
montane salamanders. \emph{Heredity} (Edinb), \textbf{119} , 265- 274.

Rubidge, E. M., Patton, J. L., \& Moritz, C. (2014). Diversification of
the Alpine Chipmunk, \emph{Tamias alpinus} , an alpine endemic of the
Sierra Nevada, California. \emph{BMC evolutionary biology}
,\textbf{14(1)} , 34.

Schonswetter, P., \& Schneeweiss, G. M.. (2019). Is the incidence of
survival in interior Pleistocene refugia (nunataks) underestimated?
Phylogeography of the high mountain plant. \emph{Ecology Evolution}
,\textbf{9} , 4078- 4086.

Schoville, S. D., \& Roderick, G. K. (2009). Alpine biogeography
of\emph{Parnassian} butterflies during Quaternary climate cycles in
North America. \emph{Molecular Ecology} , \textbf{18} , 3471- 3485.

Schoville, S. D., \& Roderick, G. K. (2010). Evolutionary
diversification of cryophilic \emph{Grylloblatta} species
(Grylloblattodea: Grylloblattidae) in alpine habitats of
California.\emph{BMC Evolutionary Biology} , \textbf{10} , 163.

Schoville, S. D., Roderick, G. K., \& Kavanaugh, D. H. (2012). Testing
the `Pleistocene species pump'in alpine habitats: lineage
diversification of flightless ground beetles (Coleoptera: Carabidae:
Nebria) in relation to altitudinal zonation. \emph{Biological Journal of
the Linnean Society} , \textbf{107} , 95- 111.

Schneeweiss, G. M., \& Schoenswetter, P. (2011). A re\selectlanguage{english}-appraisal of
nunatak survival in arctic\selectlanguage{english}-alpine phylogeography. \emph{Molecular
Ecology} , \emph{20} (2), 190-192.

Segarra\selectlanguage{english}-Moragues, J. G., Palop\selectlanguage{english}-Esteban, M., Gonzalez\selectlanguage{english}-Candelas, F., \&
Catalan, P. (2007). Nunatak survival vs. tabula rasa in the Central
Pyrenees: a study on the endemic plant species \emph{Borderea
pyrenaica}(Dioscoreaceae). \emph{Journal of Biogeography} , \textbf{34}
, 1893- 1906.

Schonswetter, P., Stehlik, I., Holderegger, R., \& Tribsch, A. (2005).
Molecular evidence for glacial refugia of mountain plants in the
European Alps. \emph{Molecular Ecology} , \textbf{14(11)} , 3547- 3555.

Sim, Z., Hall, J. C., Jex, B., Hegel, T. M., \& Coltman, D. W. (2016).
Genome-wide set of SNPs reveals evidence for two glacial refugia and
admixture from postglacial recolonization in an alpine
ungulate.\emph{Molecular Ecology} , \textbf{25} , 3696- 3705.

Simon, C., Frati, F., Beckenbach, A., Crespi, B., Liu, H., \& Flook, P.
(1994). Evolution, weighing, and phylogenetic utility of mitochondrial
gene sequences and a compilation of conserved polymerase chain reaction
primers. \emph{Annals of the Entomological Society of America}
,\textbf{87} , 651- 701.

Starcova, M., Vohralik, V., Kry\selectlanguage{polish}š\selectlanguage{english}tufek, B., Bolf\selectlanguage{ngerman}íková, B. \selectlanguage{polish}Č., \selectlanguage{english}\& Hulva,
P. (2016). Phylogeography of the Alpine shrew, \emph{Sorex
alpinus}(Soricidae, Mammalia). \emph{Journal of Vertebrate Biology}
,\textbf{65} , 107- 116.

Stehlik, I., Blattner, F. R., Holderegger, R., \& Bachmann, K. (2002).
Nunatak survival of the high Alpine plant \emph{Eritrichium nanum} (L.)
Gaudin in the central Alps during the ice ages. \emph{Molecular Ecology}
, \textbf{11(10)} , 2027- 2036.

Stoltz, M., Bauemer, B., Bouckart, R., Fox, C., Hiscott, G., \& Bryant,
D. (2019). Bayesian inference of species trees using diffusion models.
arXiv preprint arXiv:1909.07276.

Theissinger, K., B\selectlanguage{ngerman}álint, M., Feldheim, K. A., Haase, P., Johannesen, J.,
Laube, I., \& Pauls, S. U. (2013). Glacial survival and post\selectlanguage{english}-glacial
recolonization of an arctic--alpine freshwater insect
(\emph{Arcynopteryx dichroa} , Plecoptera, Perlodidae) in
Europe.\emph{Journal of Biogeography} , \textbf{40} , 236- 248.

Thiel\selectlanguage{english}-Egenter, C., Alvarez, N., Holderegger, R., Tribsch, A., Englisch,
T., Wohlgemuth, T., Colli, L., Gaudeul, M., Gielly, L., Jogan, N.,
Linder, H.P., Negrini, R., Niklfeld, H., Pellecchia, M., Rioux, D.,
Schonswetter, P., Taberlet, P., van Loo, M., Winkler, M., and Gugerli,
F. (2011), Break zones in the distributions of alleles and species in
alpine plants. \emph{Journal of Biogeography} , \textbf{38} : 772- 782.

Vuong, Q. H. (1989). Likelihood ratio tests for model selection and
nonnested hypotheses. \emph{Econometrica} , \textbf{57(2)} , 307- 333.

Wachter, G. A., Papadopoulou, A., Muster, C., Arthofer, W., Knowles, L.
L., Steiner, F. M., \& Schlick-Steiner, B. C. (2016). Glacial refugia,
recolonization patterns and diversification forces in
Alpine-endemic\emph{Megabunus} harvestmen. \emph{Molecular Ecology} ,
\textbf{25} , 2904- 2919.

Wallis, G. P., Waters, J. M., Upton, P., \& Craw, D. (2016). Transverse
Alpine Speciation Driven by Glaciation. \emph{Trends in Ecology \&
Evolution} , 31, 916-926.

Warmuth, V. M., \& Ellegren, H. (2019). Genotype\selectlanguage{english}-free estimation of
allele frequencies reduces bias and improves demographic inference from
RADSeq data. \emph{Molecular Ecology Resources} , \textbf{19} , 586-
596.

Weir, B. S., \& Cockerham, C. C. 1984. Estimating F\selectlanguage{english}-statistics for the
analysis of population structure. \emph{Evolution} , \textbf{38} , 1358-
1370.

Weng, Y. M., Yang, M. M., \& Yeh, W. B. (2016). A comparative
phylogeographic study reveals discordant evolutionary histories of
alpine ground beetles (Coleoptera, Carabidae). \emph{Ecology and
Evolution} , \textbf{6} , 2061- 2073.

Zaman, K., Hubert, M. K., \& Schoville, S. D. (2019). Testing the role
of ecological selection on colour pattern variation in the butterfly
Parnassius clodius. \emph{Molecular Ecology} , \textbf{28} , 5086- 5102.

\textbf{Data Accessibility}

All the genetic data sequenced in this study have been deposited in
National Center for Biotechnology Information (NCBI) including COI gene
in GenBank (accession number: MT676423-MT676429) and next generation
sequences reads under BioProject (PRJNA645878).

\textbf{Author Contributions}

This study was conceptually designed by YMW and SDS. The data collection
and analyses were conducted by YMW. All authors contributed to the
writing of the manuscript, with the first draft completed by YMW as part
of his dissertation.

\textbf{Supplementary Material}

The supplementary materials are deposited in the following file:
MEC\_Supplemental\_Information\_Weng etal.docx

\textbf{Figure Legends}

\textbf{Figure 1} . Sample sites of the \emph{Nebria ingens} complex in
the Sierra Nevada, California. Collecting sites of \emph{N. riversi} ,
the intermediate morphotype, and \emph{N. ingens} are depicted with
solid triangles, squares, and circles, respectively. The survey sites
where no populations were located (despite survey efforts) are denoted
by x marks. The extent of last glacial maximum is shown as a grey line
(Rood, Burbank, \& Finkel, 2011). Geographical names for each site can
be found in \textbf{Table 2} .

\textbf{Figure 2} . The Bayesian mitochondrial COI gene tree of
the\emph{Nebria ingens} complex with posterior probabilities listed
above the nodes and divergence time listed beneath the nodes. Divergence
time is estimated using a strict molecular clock (0.0113 substitutions
per lineage per million years) and is denoted by the mean and 95\%
confidence interval. The two major clades correspond to a northern clade
consisting of \emph{N. riversi} , intermediate morphotypes, and
some\emph{N. ingens} samples, and a southern clade consisting solely of
the intermediate morphotype and \emph{N. ingens} .

\textbf{Figure 3.} The SNAPP species tree of individuals from eight
sites using the \emph{fully-filtered} SNP dataset. The thick black line
shows the consensus tree, and the numbers at each node denote the
posterior probability.

\textbf{Figure 4} . Population structure of the \emph{N. ingens} complex
based on sNMF clustering analysis for K=6 (left) and K=10 (right). The
concatenated ML tree is shown on the K=6 map with bootstrap support for
major nodes. The colors of population clusters and clades reflect the
major drainage basins of the Sierra Nevada, and the drainage is roughly
circled by the dashed lines.

\textbf{Figure 5.} The model accuracies estimated using conSturct. For
the K=3 to K=8, spatial models are significantly more accurate than
non-spatial models (\emph{p} \textless{}0.05).

\textbf{Figure 6.} Stairway plots of eight selected sites including site
1 (A), site 3 (B), site 9 (C), site 10 (D), site 12 (E), site 17 (F),
site 24 (G), and site 27 (H). The solid lines, dashed lines, and dotted
lines represent the median, 12.5\% and 87.5\%, and 2.5\% and 97.5\%
confidence intervals, respectively. All sites show dramatic population
declines in past 30 thousand years.

\textbf{Figure 7.} Measurements of TMRCA (A) and\emph{F}
\textsubscript{ST} (B) for paired populations from the same (solid
circle) and different (empty circle) drainages, plotted against
geographical distance. For TMRCA (A), there is a significant difference
between same and different drainage population pairs, with no
significant correlation between TMRCA and geographical distance. The
TMRCA across different drainages mostly predates the last glacial
maximum (LGM), which is denoted by the grey box. For\emph{F}
\textsubscript{ST} (B), there is no significant difference between
population pairs within and among drainages, but the correlation between
\emph{F} \textsubscript{ST} and geographic distance is significant.

\textbf{Table 1. Hypothesis testing framework for competing glacial
refugia hypotheses.}\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllll@{}}
\toprule
Model name & Periglacial refugia & Periglacial refugia & Periglacial
refugia & Periglacial refugia & Drainage basin refugia & Drainage basin
refugia & Drainage basin refugia & Drainage basin refugia & Drainage
basin refugia & Nunatak refugia\tabularnewline
\midrule
\endhead
Submodel & Submodel & single refugium & two refugia & three refugia &
three refugia & refugia in 4 drainge basins & refugia in 5 drainge
basins & refugia in 6 drainge basins & Multiple refugia in small ice
free areas & Multiple refugia in small ice free areas\tabularnewline
Hypothesis & Hypothesis & The entire species complex shares one glacial
refugium, in the periglacial region & Two species retreated to different
periglacial refugia & \emph{Nebria riversi}, \emph{N. ingens}, and the
intermediate morphotype retreated to three different periglacial refugia
& \emph{Nebria riversi}, \emph{N. ingens}, and the intermediate
morphotype retreated to three different periglacial refugia & The
glacial refugia were within drainage basins: Tuolumne + Merced, San
Joaquin, King, and Kern + Kaweah & The glacial refugia were within
drainage basins: Tuloumne + Merced, San Joaquin, King, Kern, and Kaweah
& The glacial refugia were within drainage basins: Tuolumne, Merced, San
Joaquin, King, Kern, and Kaweah & Multiple glacial refugia persisted
within the ice free areas inside the glacial extent & Multiple glacial
refugia persisted within the ice free areas inside the glacial
extent\tabularnewline
Expected \# of clusters & Expected \# of clusters & K=1 & K=2 & K=3 &
K=3 & K=4 & K=5 & K=6 & K\textgreater{}6 &
K\textgreater{}6\tabularnewline
sNMF & sNMF & unsupported & unsupported & supported & supported &
supported & supported & supported & supported & supported\tabularnewline
spatial-model of conStruct & spatial-model of conStruct & unsupported &
unsupported & partially supported & partially supported & supported &
supported & partially supported & unsupported &
unsupported\tabularnewline
AIC (TMRCA) & AIC (TMRCA) & \emph{NA} & \emph{NA} & 2081.652 & 2081.652
& 2083.925 & 2062.997+ & 2062.997 & \emph{NA} & \emph{NA}\tabularnewline
BIC (TMRCA) & BIC (TMRCA) & \emph{NA} & \emph{NA} & 2091.909 & 2091.909
& 2094.183 & 2073.255+ & 2073.255 & \emph{NA} & \emph{NA}\tabularnewline
Model comparison & Model comparison & \emph{NA} & \emph{NA} & \emph{NA}
& \emph{NA} & =K3 & \textgreater{}K4 & =K5 & \emph{NA} &
\emph{NA}\tabularnewline
\bottomrule
\end{longtable}

+ K=5 has same set of paired population due to the lack of samples from
Merced drainage.

\textbf{Table2. Collecting information for each sample site.}\selectlanguage{english}
\begin{longtable}[]{@{}lllllll@{}}
\toprule
Code & Site name & Drainage system & Sample size (nuclear) & Sample size
(mtDNA) & Latitude & Longitude\tabularnewline
1 & Conness Lake & Tuolumne & 32 & 22 & 37.94329 &
-119.28616\tabularnewline
2 & Kuna Lake & Tuolumne & 0 & 1 & 37.84874 & -119.25929\tabularnewline
3 & Donohue Pass & Merced/San Joaquin & 16 & 16 & 37.75289 &
-119.24847\tabularnewline
4 & Lyell Peak & Merced/San Joaquin & 16 & 15 & 37.74956 &
-119.25991\tabularnewline
5 & Ritter Range & San Joaquin & 0 & 0 & 37.71442 &
-119.20355\tabularnewline
6 & Ottoway Lake & Merced & 0 & 3 & 37.63967 & -119.39713\tabularnewline
7 & Recess Lakes & San Joaquin & 3 & 10 & 37.43124 &
-118.78462\tabularnewline
8 & Ruby Lake & San Joaquin & 16 & 12 & 37.41155 &
-118.77399\tabularnewline
9 & Italy Lake & San Joaquin & 10 & 10 & 37.34084 &
-118.76735\tabularnewline
10 & Selden Pass & San Joaquin & 14 & 17 & 37.28632 &
-118.88537\tabularnewline
11 & Piute Pass & San Joaquin & 22 & 18 & 37.23683 &
-118.679944\tabularnewline
12 & Lamarck Lakes & San Joaquin & 17 & 20 & 37.20853 &
-118.65117\tabularnewline
13 & Hungry Packer Lake & San Joaquin/Kings & 11 & 11 & 37.16094 &
-118.64077\tabularnewline
14 & Treasure Lake & Kings & 9 & 11 & 37.14502 &
-118.58003\tabularnewline
15 & Sam Mack Lake & Kings & 14 & 17 & 37.11192 &
-118.51174\tabularnewline
16 & Dusy Pass & Kings & 2 & 6 & 37.09534 & -118.53421\tabularnewline
17 & Taboose Pass & Kings & 9 & 13 & 36.97993 &
-118.39877\tabularnewline
18 & Sixty Lakes & Kings & 4 & 5 & 36.82479 & -118.43382\tabularnewline
19 & Sphinx Lakes & Kings/Kern & 2 & 2 & 36.71427 &
-118.51487\tabularnewline
20 & Milly's Footpass & Kings/Kern & 12 & 12 & 36.69039 &
-118.43453\tabularnewline
21 & North Forester & Kings/Kern & 4 & 4 & 36.68882 &
-118.37648\tabularnewline
22 & South Forester & Kern & 16 & 25 & 36.63116 &
-118.34489\tabularnewline
23 & Wright Lakes & Kern & 12 & 16 & 36.59707 & -118.6658\tabularnewline
24 & Pear Lake & Kaweah & 12 & 14 & 36.59707 & -118.6658\tabularnewline
25 & Crabtree Lakes & Kern & 14 & 14 & 36.54445 &
-118.32363\tabularnewline
26 & Monarch Lake & Kaweah/Kern & 6 & 6 & 36.50159 &
-118.24711\tabularnewline
27 & Army Pass & Kern & 38 & 19 & 36.41429 & -118.54783\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table3. Successful copulation rate of \emph{N. riversi},
intermediate morphotypes, and \emph{N. ingens}}\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
& \emph{N. riversi} male & intermediate male & \emph{N. ingens}
male\tabularnewline
\emph{N. riversi} female & 3/3 & 7/7 & 1/5\tabularnewline
intermediate female & 4/4 & 6/7 & 0/5\tabularnewline
\emph{N. ingens} female & 1/3 & 6/7 & 1/5\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figures-new/Screen-Shot-2020-09-09-at-5.00.00-PM}

\end{center}
\end{figure}

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

