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.