Introduction
Species delimitation has increasingly greater biological, social, economic, and political consequences amidst the rapid habitat and biodiversity losses of our planet (Coates et al. 2018; Stantonet al. 2019). While species comprise the operational units of modern conservation initiatives (Mace 2004; Magurran 2013; MacDonald et al. 2017; IPBES 2019), the procedures employed in species delimitation and identification suffer from a number of limitations, in part because the processes of speciation are affected by diverse historical, genetic, ecological, and stochastic factors (Sites & Marshall 2003; Wiens 2004; de Queiroz 2007; Wilkins 2009; Loera et al. 2012). This has led to a plethora of species concepts that differ in their emphases on morphological, genetic, or ecological information for resolving divergences and distinguishing species, each of which may produce inconsistent or contradictory species assessments (Mayr 1942, 1957, 1963; Mayden 1997; Wiens 2004; de Queiroz 1998; Freudenstein et al. 2016).
Assessing species boundaries can be particularly challenging for taxa that are not easily morphologically distinguished or have superficially similar habitat preferences. In this context, quantitative comparisons of genomic and ecological divergences can aid in the resolution of historical and ongoing speciation events, particularly for allopatric or parapatric sister species (Carstens et al. 2013). Furthermore, the integration of these two sources of information permits insight into whether ecologically-based divergent selection has contributed to speciation and reinforcement of a species’ genomic integrity (Sperling 2003; Graham et al. 2004). For example, significant niche divergence between sister species would suggest that ecologically-based divergent selection is a principal process underlying speciation (i.e., ecological speciation), while niche conservatism suggests that other processes are implicated (e.g., speciation without selection or mutation-order speciation, sensu Nosil 2012). Incorporation of ecological information into genetic species delimitation frameworks can thereby provide valuable inferences on modes of speciation and strengthen phylogenetic inferences when traditional approaches have not convincingly done so.
The development of high-throughput DNA sequencing techniques such as RADseq (Baird et al. 2008) and related methods have enabled genotyping of thousands of single nucleotide polymorphisms (SNPs), greatly extending our capacity to detect recent, fine-scale genomic divergences in non-model organisms (Andrews et al. 2016). Such genomic data have proven valuable for clarifying population dynamics and re-assessing species limits in taxa that have been historically difficult to characterize using other approaches, either due to recent diversification, morphological ambiguity, historical introgression, or some combination of these factors (Hohenlohe et al. 2013; Wagneret al. 2013; Escudero et al. 2014; Vargas et al.2017; Abdelkrim et al. 2018; Hinojosa et al. 2019; Hundsdoerfer et al. 2019). Alongside genomic advances, continued development of ecological niche models (ENMs) has facilitated an integrative approach to inferring processes that contribute to ecological diversification and reinforcement of species (Schluter 2001, 2009; Manel et al. 2003; Balkenhol et al. 2016). ENMs, generally parameterized by relating species’ occurrences to geographic and environmental factors, may be used to quantify ecological niches, habitat associations, and potential geographic distributions of single species (Austin 1985; Peterson 2001; Guisan & Zimmermann 2000; Elithet al. 2006, 2011; Phillips et al . 2006; Zimmermann et al. 2010). However, ENMs may also be used in a comparative framework to assess niche divergence and conservatism among recently diverged evolutionary lineages, often delineated on the basis of genetic data (Sites & Marshall 2003; Graham et al. 2004, Kozak & Wiens 2006, Bond & Stockman 2008, Jezkova et al. 2009, Loera et al.2012; Newton et al. 2020).
The butterfly genus Speyeria Scudder, 1872 (Lepidoptera: Nymphalidae) is well known for its phenotypic variability and ambiguous evolutionary relationships among component lineages (dos Passos & Grey 1947; Moeck 1975; Dunford 2009). A total of 16 species are currently recognized in Speyeria , as well as over 110 morphologically variable subspecies and several species complexes with poorly understood evolutionary relationships (dos Passos & Grey 1947; Scott et al.1998; Dunford 2009; Pelham 2019). Among these, S. hesperis(Edwards, 1864) and S. atlantis (Edwards, 1862) form a large complex containing 26 subspecies (Pelham 2019). These include five subspecies in S. atlantis that are broadly distributed in conifer woodlands across North America from the Rocky Mountains to Newfoundland. More variation is taxonomically recognized in S. hesperis , which has 21 subspecies occurring in drier meadows and open forests throughout western North America and east to South Dakota and southeastern Manitoba (Pelham 2019). The two species contact each other in mixed forest areas from Manitoba to British Columbia and south along the Rocky Mountains to Colorado, exhibiting substantial morphological similarity between species in some areas as well as variation within species (dos Passos & Grey 1947; Moeck 1975; Dunford 2009).
Some taxonomic treatments have considered S. hesperis to be a subspecies of S. atlantis based on overall morphological similarity (Grey 1951; Miller and Brown 1981; Hammond et al.2013). Recent work recognizes these taxa as distinct species, based on assessments of morphological and genetic divergence, as well as an apparent lack of hybridization between sympatric populations (Campbellet al. 2017; de Moya et al. 2017; Campbell et al.2019; Riva et al. 2019; Thompson et al. 2019). However, the taxa remain difficult to reliably identify using morphology alone (Scott et al. 1998; Opler & Warren 2005). Additionally, Campbellet al. (2019) have recently shown, based on a limited number of specimens, substantial genomic divergence in SNPs between S. hesperis populations occurring north and east of the Rocky Mountains (throughout British Columbia, Alberta, and Saskatchewan in Canada and Montana and South Dakota in the United States; hereafter referred to as the “northern lineage”), and those occurring throughout the southwestern and Great Basin regions of the United States (“southern lineage”). Further phylogenetic complexity is provided by interactions with species that have not historically been considered part of this complex, including intermediates between S. hesperis and S. zerene (Boisduval, 1852) (Campbell et al. 2019). This may have important implications for conservation initiatives for S. zerene , which has multiple subspecies experiencing significant population declines in western regions where S. hesperis andS. zerene co-occur (McHugh et al. 2013; Sims 2017). While these genetic assessments have helped clarify some taxonomic ambiguities, there have been no attempts to assess whether genomic divergences correspond to variation in ecological niches and habitat associations.
Our objective is thus to provide a proof-of-concept for the use of ecological modelling to strengthen genomic assessments of species boundaries and to clarify some of the extrinsic factors involved in the diversification of the S. atlantis-hesperis species complex. We use de novo SNPs to recover distinct genetic clusters of populations that maintain their genomic integrity in regions of contact (Sperling 2003), testing alternate hypotheses on species delimitation and phylogeographic factors that may have contributed to recovered genetic patterns. We additionally use ENMs to compare ecological niches of evolutionary lineages identified on a genomic basis to infer whether ecologically-based divergent selection is a likely contributor to their speciation and reinforcement of genomic integrity. Our integration of these methods demonstrates broad utility for the reconciliation of species concepts, such as the genomic integrity (Sperling 2003) and ecological species concepts (Van Valen 1976; Andersson 1990; Nosil 2012), that should contribute to stable species delimitations.