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ö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ński, Sękiewicz, Walas, Boratyński, & Dering, 2019; Pan, Hü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 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 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 alpineNebria , members of the N. ingens complex has atrophied
hind wings that impede flight, restricting gene flow and distributional
range size (Kavanaugh 1985). In fact, the 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 N. ingens and 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 Nebria ingenscomplex 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 (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; Figure 1 ). Finally, persistence in drainage basins is
proposed, as we expect ancestral populations of the N. ingenscomplex were dependent on these geographical features to obtain suitable
microhabitats.