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.