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.