Introduction
Cataloguing and describing species is a crucial first step towards
understanding Earth’s biodiversity. Species are the foundation of
biological questions including evolutionary processes, ecological
systems, physiological mechanisms, and are key to formulating
conservation efforts. However, determining species limits remains a
contentious topic (Freudenstein et al., 2016; Hey, 2001; Wake, 2006),
particularly for species delimitation in taxa prone to cryptic
diversity. (i.e., morphologically indistinguishable species yet exhibit
extensive molecular divergence; Bickford et al., 2007). Morphological
homogeneity may be the result of recent divergence (i.e., not enough
time has elapsed for morphological traits to evolve) or niche
conservatism (i.e., two species inhabit ecologically similar niches but
geographic separation led to genetic divergence due to lack of gene
flow; Wiens & Graham, 2005). Cryptic species cannot be distinguished by
morphology alone; thus, additional data types, such as molecular and
ecological data, are necessary for robust species delimitation (Bickford
et al., 2007; Stockman & Bond, 2007). For example, recent studies using
molecular data in addition to other types of data have unveiled cryptic
species across many animal taxa, including birds (e.g., Battey &
Klicka, 2017; Garg et al., 2016), arthropods (e.g., Bond et al., 2001;
Bond & Sierwald, 2002; Bond & Stockman, 2008; Derkarabetian et al.,
2011; Derkarabetian & Hedin, 2014; Garrick et al., 2018; Starrett et
al., 2018; Starrett & Hedin, 2007; Zhang & Li, 2014), amphibians
(e.g., Chan et al., 2017; Moritz et al., 2016; Ortega-Andrade et al.,
2015; Ramírez-Reyes et al., 2017; Reynolds et al., 2012; Weisrock &
Larson, 2006), cnidarians (e.g., Holland et al., 2004), and annelids
(e.g., Barroso et al., 2010). These studies highlight that species
delimitation utilizing only morphology potentially underestimates
species diversity in taxa with relative morphological homogeneity.
Spiders in the infraorder Mygalomorphae (Bond et al. 2012; Opatova et
al. 2019) present a perplexing situation for robust species delimitation
when compared to other more diverse araneomorph spiders (as well as
other taxa). Araneomorph spider species delimitation relies
predominately on morphological criteria (e.g., distinctive differences
in male genitalia, body size, etc.; Dupérré and Tapia 2015; Richardson
2016; also see Bond et al. in review for summary) whereas mygalomorph
spiders are relatively morphologically homogeneous (Bond & Stockman,
2008; Hamilton et al., 2014; Harvey et al., 2018; Hendrixson & Bond,
2005a, 2005b; Huey et al., 2019; Opatova & Arnedo, 2014b) and subject
to significant population genetic structuring at microgeographical
scales (Bond et al., 2001; Hedin et al., 2015; Starrett et al., 2018).
Mygalomorph life history traits that include limited dispersal abilities
(with few exceptions; see Coyle 1983; 1985), relatively long life spans
(10-30 years), habitat specialization, and site fidelity altogether make
them ideal organisms for studying population divergence and
phylogeography (Hamilton et al., 2014; Hedin et al., 2013; Hendrixson et
al., 2013; Hendrixson & Bond, 2005a, 2005b; Opatova & Arnedo, 2014a;
Satler et al., 2011; Starrett & Hedin, 2007; Stockman & Bond, 2007).
Although morphology tends to under-split mygalomorph species, DNA
barcoding and related approaches such as GMYC tend to over-split species
owing to their population structuring at very fine geographical scales;
it has been long recognized (Bond et al. 2001) that mygalomorph
populations are highly structured over relatively short distances.
Hamilton et al. (2014) posited that GMYC greatly overestimated species
diversity in the tarantula genus, Aphonopelma (recognizing 114
species), whereas more integrative type approaches recovered fewer
species (34). In addition to single locus approaches, multispecies
coalescent methods using many loci are more inclined to detect
population structure rather than speciation events (Sukumaran &
Knowles, 2017). As a result, such attributes of mygalomorph
morphological similarity and population structure typify a system
requiring more integrative approaches when delimiting species; that is,
it is important to evaluate molecular, geographical, ecological, and
morphological lines of evidence as opposed to relying on only one data
type, analysis, or conceptual approach (e.g., Bond and Stockman 2008;
Hamilton et al. 2014; Hedin, Carlson, and Coyle 2015).
The focus of our study is the Antrodiaetus unicolor mygalomorph
species complex. Like many related mygalomorph spider species, they have
relatively long life spans (> 10 years), high site and
microhabitat fidelity (e.g. mesic forests and stream banks), and limited
dispersal capabilities (Coyle, 1971; Hendrixson & Bond, 2005a, 2005b).
Sympatry within the species complex occurs where Antrodiaetus
unicolor and A. microunicolor are co-distributed across part of
their ranges in the eastern United States (Figure 1; Hendrixson and Bond
2005a; 2005b). Antrodiaetus microunicolor is found only in the
southwestern region of the Appalachian Mountains, whereas A.
unicolor is centered in the central and southern regions of the
Appalachian Mountains with peripheral populations extending as far west
as the Ozarks in Arkansas and east near the Atlantic Coast, with
populations as far south as the Gulf Coast, and as far north to the
central-eastern corner of Pennsylvania. In the mesic forests of the
Appalachians these spiders are primarily found along creek banks or
steep ravines, whereas in the peripheral populations they are isolated
in small ravines in humid hardwood forests (Coyle, 1971; Hendrixson &
Bond, 2005a). These spiders build subterranean burrows covered by a
unique collapsible collar door where they sit and wait for prey; as a
result, they rarely depart from their burrows unless disturbed, or to
seek a mate in the case of mature males (Coyle, 1971; Hedin et al.,
2019).
Antrodiaetus spiders show both morphological stasis and variation
across their distribution (Coyle, 1971). Individuals from localities
separated by hundreds of kilometers may be morphologically
indistinguishable, yet spiders at the same location may exhibit
significant disparities in size and coloration (Coyle, 1971; Hendrixson
& Bond, 2005a). Due to difficulty interpreting this intraspecific
variation, Coyle (1971) was conservative when revising the genus by
maintaining all populations of A. unicolor as one species. Later,
Hendrixson and Bond (2005a) used morphological and behavioral data to
distinguish two forms different in size, coloration, and setal patterns
from the Coweeta Long Term Ecological Research site in southwestern
North Carolina and described A. microunicolor as a new species. A
subsequent molecular analysis using two genetic markers (mtDNA gene
cytochrome oxidase I and nuclear ribosomal RNA gene 28S) to evaluate
species boundaries between these two forms showed that A.
unicolor is paraphyletic with respect to A. microunicolor(Hendrixson & Bond, 2005b), suggesting that these spiders possibly
comprise a complex of multiple cryptic species.
With the advent of next-generation sequencing methods, it is now
feasible to generate genomic-scale data for non-model organisms (Baird
et al., 2008; Faircloth et al., 2012; Lemmon et al., 2012). These much
larger, more comprehensive data sets provide a framework for more
rigorous tests of species boundaries in systems where previously only a
few targeted loci were available. RADseq (restriction-site associated
DNA sequencing) approaches are one of the most widely used techniques
for generating a reduced representation of the nuclear genome with
extensively sampled homologous loci (Andrews et al., 2016; Baird et al.,
2008). Various RADseq techniques have been utilized for addressing
population genetic and phylogeographic studies (Andrews et al., 2016;
Baird et al., 2008; Emerson et al., 2010), reconstructing phylogenetic
relationships among both closely (e.g., Eaton et al. 2016) and distantly
(e.g., Leaché et al. 2015) related species, and evaluating species
boundaries and speciation processes (e.g., Battey & Klicka, 2017; De
Jesús-Bonilla et al., 2019; Delgado‐Machuca et al., 2019).
We generated RADseq data using 3RAD, a three-enzyme protocol that
reduces DNA chimaera and adapter-dimer formation (Bayona-Vásquez et al.,
2019), to investigate species boundaries and phylogenetic relationships
within the A. unicolor species complex. Specifically, we employed
the cohesion species concept (Templeton, 1989), which defines a species
as a set of populations that derives from a single evolutionary lineage
and meet the criteria of being genetically exchangeable and/or
ecologically interchangeable. This integrative approach is particularly
insightful when evaluating species limits in morphologically cryptic
taxa with high population genetic structuring (e.g., Bond & Stockman,
2008; Hendrixson et al., 2013; Stockman & Bond, 2007). To evaluate
cohesion species identity within the A. unicolor species complex,
we assessed the amount of genetic population structure using clustering
analyses to identify the number of evolutionary lineages (i.e., genetic
exchangeability), and constructed niche-based distribution or ecological
niche models (ENMs) for each lineage. ENMs were then compared to
evaluate overlap/similarity to assess potential ecological
interchangeability.