# A Framework for Resolving Cryptic Species: A Case Study from the Lizards of the Australian Wet Tropics

A Framework for Resolving Cryptic Species: A Case Study from the Lizards of the Australian Wet... Abstract As we collect range-wide genetic data for morphologically-defined species, we increasingly unearth evidence for cryptic diversity. Delimiting this cryptic diversity is challenging, both because the divergences span a continuum and because the lack of overt morphological differentiation suggests divergence has proceeded heterogeneously. Herein, we address these challenges as we diagnose and describe species in three co-occurring species groups of Australian lizards. By integrating genomic and morphological data with data on hybridization and introgression from contact zones, we explore several approaches—and their relative benefits and weaknesses—for testing the validity of cryptic lineages. More generally, we advocate that genetic delimitations of cryptic diversity must consider whether these lineages are likely to be durable and persistent through evolutionary time. Exome capture, cryptic species, phylogeography, species delimitation, squamates, taxonomy Cryptic species, or taxa that are morphologically similar but genetically divergent, exemplify the two major challenges of species delimitation. First, species form on a continuum (Darwin 1859; Mayr 1942; Mallet 1995; De Queiroz 2007). As populations differentiate across space and time, they gradually become more divergent. As reflected in the debate over defining operational taxonomic units via DNA barcoding (Moritz and Cicero 2004), deciding how much divergence is sufficient to name lineages as species can be arbitrary. As with morphologically distinct lineages, diagnosing cryptic lineages presents this challenge because they fall through the full range of the divergence continuum (Hedgecock and Ayala 1974; Gómez et al. 2002; McDaniel and Shaw 2003). Second, speciation proceeds heterogeneously across many axes. We typically recover correlations across axes—e.g., rates of trait evolution such as song and mitochondrial divergence (Winger and Bates 2015). However, when axes of differentiation are discordant—for example, when phenotypic disparity is high and genetic divergence is low – the status of lineages becomes ambiguous. By definition, cryptic lineages have diverged heterogeneously (Bickford et al. 2007)—they are genetically-distinct groups that exhibit little or no morphological divergence. The taxonomic process of naming a species is a binary exercise—either a lineage is recognized as a species or not—and accounting for heterogeneity in this binary framework remains a challenge. Biodiversity researchers increasingly face these challenges because we are increasingly discovering new cryptic lineages (Bickford et al. 2007). A confluence in genetic advances and broader geographic sampling has led to rapid increases in the number of cryptic species, such that a quarter of articles published in Zoological Record Plus mention cryptic species (Bickford et al. 2007; see also de León and Poulin 2016). Cryptic species comprise a significant proportion of the diversity in some regions (e.g., tropics; Smith et al. 2008) and taxonomic groups (e.g., reptiles; Oliver et al. 2010), and recently, putative cryptic species have been identified in high-profile threatened species like orangutans (Nater et al. 2017). These findings, and their implications for evolutionary biology and conservation (Frankham et al. 2012), emphasize the need for a more rigorous framework to assess the taxonomic status of cryptic lineages (Adams et al. 2014; Struck et al. in press). In this work, we propose a framework that diagnoses those cryptic lineages that are expected to be sufficiently durable to contribute to build-up of diversity over time and space. Because speciation is a continuum, we expect that many nascent species are lost to hybridization and extinction as part of the protracted speciation process (Rosindell et al. 2010; Dynesius and Jansson 2014). As such, we adopt the biological species concept (BSC, (Mayr 1942)), which defines species as units that exhibit barriers to reproduction and are thus more likely to persist through time. Although some might find this definition overly restrictive, we apply it here in hopes of avoiding “taxonomic over-inflation” (Isaac et al. 2004). However, how do we diagnose populations that are likely to have substantial reproductive isolation (RI)? In allopatry, the degree of morphological difference is expected to correlate with the extent of RI (Mayr 1942; Bolnick et al. 2006; Funk et al. 2006). Thus, when genetic and phenotypic divergence concur, species delimitation is typically uncontroversial. For cryptic species, where we cannot use phenotypic divergence as a proxy for RI, we must instead use multiple lines of evidence to assess likelihood of strong RI. A popular approach to assess cryptic lineages is to apply statistical species delimitation to multilocus genetic data (Fujita et al. 2012; Carstens et al. 2013), which some argue makes species delimitation more objective (Rannala 2015). An illuminating line of research has explored the parameter space under which these coalescent-based methods are expected to return statistically robust results (Zhang et al. 2011; Olave et al. 2014). What remains to be seen is if these statistically robust lineages are also biologically robust (Sukumaran and Knowles 2017). That is, will newly delimited lineages remain distinct through changing geographies and environments, or will they be mere evolutionary ephemera lost to hybridization and/or extinction (Seehausen et al. 2008; Rosenblum et al. 2012; Dynesius and Jansson 2014)? A stronger and more direct approach to delimitation is to test for strongly restricted gene flow in sympatry or parapatry (Richardson et al. 1986; Adams et al. 2014). This can be done either by sampling a modest number of individuals in sympatry or by assessing the extent of genetic introgression through intensive analysis of contact zones between parapatric taxa. However, when candidate taxa are allopatric, assessing the likelihood of strong RI is even more challenging, as long recognized under operational versions of the BSC. To the extent that RI is time-dependent (Coyne and Orr 1989; Sasa et al. 1998; Fitzpatrick 2002; Roux et al. 2016), another approach is to extrapolate from closely-related taxa where the relationship between divergence and extent of RI has already been determined. In contrast to a “bar-coding gap” (Hebert et al. 2004b), this approach uses genome-scale evidence to assess the likelihood of strong RI rather than patterns of genetic divergence across nominal species vs. populations. We explore these multiple approaches to species delimitation through the study of morphologically cryptic, phylogeographic lineages in the lizards of the Australian wet tropics (AWT). The AWT is a narrow region of rainforest in northeast Queensland, Australia (Fig. 1). During repeated glacial cycles in the Quaternary (Graham et al. 2006), the rainforest and, accordingly, the species endemic to these rainforests, were split across two major refugia. Populations of these rainforest species diverged across these refugia; comparative data have recovered deep phylogeographic splits within species across more than twenty taxa (Moritz et al. 2009). Morphological analyses show limited phenotypic divergence among phylogeographic lineages (Schneider and Moritz 1999; Hoskin et al. 2005; Hoskin et al. 2011). Subsequent contact zone studies showed that some of these phylogeographic lineages are reproductively isolated (Phillips et al. 2004; Hoskin et al. 2005; Singhal and Moritz 2013). Figure 1. View largeDownload slide Geographic distribution of lineages in the three groups; boundaries were inferred after sequencing an average of 27 individuals for mtDNA and 12 for nDNA (Supplementary Table S1, Fig. S2 available on Dryad). Crosses represent localities for samples included in this study. A) The five lineages in the ‘Carlia rubrigularis’ group: C. wundalthini, C. rubrigularis N, C. rubrigularis S, C. rhomboidalis N, and C. rhomboidalis S, B) the four lineages in the ‘Lampropholis coggeri’ group: L. coggeri N, L. coggeri C, L. coggeri S, and L.coggeri EU, and C) the four lineages in the ‘L. robertsi’ group: L. robertsi CU, L. robertsi BFAU, L.robertsi TU, L. robertsi BK. In the inset map, the distribution of the rainforest is shown in light green. Like many phylogeographic lineages, these lineages are geographically circumscribed, and their ranges either are geographically proximate or narrowly overlapping. Pictures courtesy of B. Phillips, C. Peng, and S. Zozaya (L-R). Figure 1. View largeDownload slide Geographic distribution of lineages in the three groups; boundaries were inferred after sequencing an average of 27 individuals for mtDNA and 12 for nDNA (Supplementary Table S1, Fig. S2 available on Dryad). Crosses represent localities for samples included in this study. A) The five lineages in the ‘Carlia rubrigularis’ group: C. wundalthini, C. rubrigularis N, C. rubrigularis S, C. rhomboidalis N, and C. rhomboidalis S, B) the four lineages in the ‘Lampropholis coggeri’ group: L. coggeri N, L. coggeri C, L. coggeri S, and L.coggeri EU, and C) the four lineages in the ‘L. robertsi’ group: L. robertsi CU, L. robertsi BFAU, L.robertsi TU, L. robertsi BK. In the inset map, the distribution of the rainforest is shown in light green. Like many phylogeographic lineages, these lineages are geographically circumscribed, and their ranges either are geographically proximate or narrowly overlapping. Pictures courtesy of B. Phillips, C. Peng, and S. Zozaya (L-R). Herein, we focus on three species groups—the ‘Carlia rubrigularis’, ‘Lampropholis coggeri’, and ‘Lampropholis robertsi’ groups—which are part of the broader radiation of Eugongylus lizards (family: Scincidae) (Skinner et al. 2011). These groups are ecologically-similar; all are small (30–55 mm) skinks of the leaf-litter in the rainforests of the AWT. Previous multilocus analyses revealed several phylogeographic lineages within each of these nominal taxa (Dolman and Moritz 2006; Bell et al. 2010). Like many phylogeographic units, these lineages are mostly morphologically cryptic and geographically circumscribed, and their ranges are either geographically proximate or narrowly overlapping. With an eye to integrative taxonomy (Padial et al. 2010), we synthesize data on genetics, morphology, and reproductive isolation assessed in contact zones to resolve the species status of lineages within these three groups and to formally revise their taxonomy. More generally, we use these lizards as a data-rich case study to explore the challenges of delimiting species among cryptic lineages that are parapatric or allopatric. Methods Sampling, Data Collection, and Data Processing In this study, we analyze genetic data for individuals across three species groups across five nominal species and 13 putative lineages. The ‘Carlia rubrigularis’ group consists of five lineages: C. rubrigularis, northern Wet Tropics (N); C. rubrigularis, southern Wet Tropics (S); C. rhomboidalis, northern mid-east Queensland (N); C. rhomboidalis, southern mid-east Queensland (S); and C. wundalthini at Cape Melville (Fig. 1). The ‘Lampropholis coggeri’ group consists of four lineages: L. coggeri, northern Wet Tropics (N); L. coggeri, central Wet Tropics (C); L. coggeri, southern Wet Tropics (S); and L. coggeri in the Mt Elliot uplands (EU). The montane ‘Lampropholis robertsi’ group consists of four allopatric lineages: L. robertsi, Carbine Tableland uplands (CU); L. robertsi, Thornton Peak uplands (TU); L. robertsi, Mt Bellenden Ker (BK); and L. robertsi, Mt Bartle Frere and southern Atherton Tablelands (BFAU). These lineages had been previously described through the sequencing of an average of 27 geographically dispersed individuals per lineage for mtDNA and 12 for multilocus nDNA (Dolman and Moritz 2006; Bell et al. 2010). Because these lineages were sampled extensively in previous genetic analyses, and because they are circumscribed geographically (Fig. 1), we limited our sampling to a few individuals per lineage (Fig. 1; Supplementary Tables S1 and S2 available on Dryad at http://dx.doi.org/10.5061/dryad.g7v1b). Further, across these lineages, three lineage-pairs meet in hybrid zones, two of which are very narrow ($$<$$1 km) and the other of which is wider ($$<$$10 km) (Singhal and Moritz 2013; Singhal and Bi 2017). Given that introgression is geographically limited relative to lineage ranges, we selected individuals well removed from contemporary contact zones to estimate species trees, test for lineage-wide introgression, and to apply statistical delimitation methods. We further included two outgroups: C. storri for the ‘C. rubrigularis’ group and L. amicula for the ‘L. coggeri’ and ‘L. robertsi’ groups. For L. coggeri N and C, we supplemented our sampling by including previously-published transcriptome data (Singhal 2013), which we analyzed using the same approach outlined below. In total, we sampled 25 individuals (Supplementary Tables S1 and S2 available on Dryad). We sequenced a homologous set of 3320 loci across all individuals using an exome capture approach previously described (Bragg et al. 2016). Briefly, we identified homologous exons across the Eugongylus skink clade from transcriptome data of three species. After filtering the exons for GC-content and length, we included probes specific to each of the three species in the capture array. For each sample, we extracted DNA (Sunnucks and Hales 1996) and generated uniquely-barcoded libraries (Meyer and Kircher 2010). Individuals were then pooled in equimolar amounts along with other Eugongylus group taxa (Bragg et al. 2018), and we captured our target exons using a SeqCap EZ Developer Probes following manufacturer”s instructions. The captured libraries were subsequently sequenced along with samples for other projects on the Illumina HiSeq2000 and 2500 for 100 bp paired-end reads. Our bioinformatics pipeline is designed for both population genetic and phylogenetic analyses and thus generates both variant data for each lineage and locus alignments across all haplotypes in a given group. The basic approach follows, with further details available in Singhal et al. (2017). Following de-multiplexing, we trimmed low quality sequence and adaptor sequence from reads using Trimmomatic v0.36 and merged overlapping reads using pear v0.9.10 (Zhang et al. 2013; Bolger et al. 2014). We generated an assembly for each individual with Trinity v2.3.2 and identified assembled contigs homologous to our original targets with blat v36x1 (Kent 2002; Grabherr et al. 2011). For each exon, we picked the best matching contig across all individuals in the lineage to generate a pseudo-reference genome. We then annotated the coding sequence for these exons using exonerate v2.2.0 (Slater and Birney 2005). To generate variant data, we mapped reads from each individual to the pseudo-reference genome using bwa v0.7.12 (Li and Durbin 2009), called variant and invariant sites using gatk v3.6 Unifiedgenotyper, filtered out sites with quality $$<$$20, quality depth (QD) $$<$$5, and coverage $${<}20{\times}$$, and then determined haplotypes using gatkReadBackedPhasing (McKenna et al. 2010). Eight percent of sites were unphased; for these sites, we randomly phased them with respect to other phased blocks. We then generated multi-lineage alignments across all haplotypes with mafft v7.294 (Katoh et al. 2002). Like many phylogeographic studies, we first identified these lineages by sequencing mitochondrial loci. For these taxa, mtDNA and nDNA are highly congruent except within narrow contact zones, and mtDNA provides our most complete understanding of geographic limits (Dolman and Moritz 2006; Bell et al. 2010). Accordingly, we downloaded all geo-referenced NADH dehydrogenase subunit 4 (ND4) data for these groups from GenBank (Supplementary Table S3 available on Dryad). We aligned the data using mafft and identified the coding sequence using exonerate. Phylogenetic Analyses We reconstructed the evolutionary history of 13 lineages and two outgroups by using starbeast2 v0.13.5, a coalescent multilocus method (Ogilvie et al. 2017). Individuals were assigned to lineages following their ‘putative lineage’ designations (Supplementary Table S2 available on Dryad). We filtered loci to only include those that were $$\geqslant$$75% complete across samples and to remove loci $$\geqslant$$1500 bp because longer loci are more likely to capture recombination events. Because of the computational demands of running starbeast2, we generated three random subsamples of 200 loci each from the remaining loci. We ran starbeast2 on these random subsets for 500e6 generations sampling every 1e5 generations. Each locus was assigned to its own strict clock and a GTR model of molecular evolution. Because we lack robust age constraints for nodes in this species tree, we instead inferred branch lengths in units of substitutions per site. For the mitochondrial phylogenetic analysis, we determined the best-fitting partitioning strategy using PartitionFinder2 (Lanfear et al. 2016). We then inferred the mtDNA gene trees using MrBayes v3.2.6, running two runs of four chains each for 50e6 steps (Ronquist et al. 2012). We set the branch length prior to exponential (100); the default prior overestimates branch lengths when the majority of bifurcations occur within-species (Brown et al. 2009). Population Genetic Analyses Our population genetic analyses were aimed at describing basic patterns of diversity, divergence, and current and historical introgression among the lineages in these groups. For each lineage, we calculated within-lineage genetic diversity ($$\pi$$; Nei and Li 1979) for both the exome and mtDNA data, across silent sites only. For each pairwise-comparison between lineages in a species group, we calculated raw and net divergence ($$d_{xy}$$ and $$d_{a})$$ for the exome and mtDNA data (Nei and Li 1979), again across silent sites only. The patterns of divergence and diversity among individuals within a lineage and across lineages confirm our lineage assignments based on mtDNA and previous multilocus nuclear data (Dolman and Moritz 2006; Bell et al. 2010). To test for historical introgression, we used the D-statistic (Durand et al. 2011). For the topology [(P1, P2), P3), outgroup], the D-statistic distinguishes if P1 and P3 exhibit incomplete lineage sorting or introgression by comparing site patterns across these four tips. We conducted this test across all possible comparisons within species groups, except for sister taxa, which cannot be accommodated in the D-statistic framework. Based on the overall species tree topology (Fig. 2), we determined the appropriate species to use as P2 (Supplementary Table S5 available on Dryad). For all Lampropholis comparisons, we used L. amicula as the outgroup, and for Carlia, we used C. storri. Because of our limited within-lineage sampling, we implemented the version of the test designed for fixed variants. Thus, we removed any site that was missing or polymorphic for any lineage in a species group. For the remaining sites, we calculated the D-statistic across all possible species comparisons. To assess significance, we calculated the standard deviation across 200 bootstraps and used a one-tailed $$z$$-test (Eaton and Ree 2013). For L. robertsi, we further tested for introgression using the D$$_{\mathrm{FOIL}}$$ approach designed for five-taxon symmetric topologies (Pease and Hahn 2015); we could not apply this method to other groups because their topologies are asymmetric. This method confirmed our D-statistic results, so we do not discuss them further. Figure 2. View largeDownload slide Species tree for the included taxa as inferred using STARBEAST2 with 200 randomly selected exons. This topology and branching times are robust across multiple random samples (Supplementary Fig. S2 available on Dryad). Nodes with $$\geqslant$$95% local posterior probability are indicated with filled circles. Nodes for which BPP inferred $$\geqslant$$95% probability of a speciation event are shown by open circles. STACEY results match BPP results and are not shown. Arrows indicate pairwise relationships for which there is significant evidence for historical introgression (Supplementary Table S6 available on Dryad). The phylogeographic lineages included in this study all diverged over a relatively narrow span of time. Figure 2. View largeDownload slide Species tree for the included taxa as inferred using STARBEAST2 with 200 randomly selected exons. This topology and branching times are robust across multiple random samples (Supplementary Fig. S2 available on Dryad). Nodes with $$\geqslant$$95% local posterior probability are indicated with filled circles. Nodes for which BPP inferred $$\geqslant$$95% probability of a speciation event are shown by open circles. STACEY results match BPP results and are not shown. Arrows indicate pairwise relationships for which there is significant evidence for historical introgression (Supplementary Table S6 available on Dryad). The phylogeographic lineages included in this study all diverged over a relatively narrow span of time. Finally, we collated previously-published data on reproductive isolation at three contact zones: L. coggeri N and C, L. coggeri C and S, and C. rubrigularis N and S (Phillips et al. 2004; Dolman 2008; Singhal and Moritz 2012, 2013; Singhal and Bi 2017). These studies sampled densely through each contact zone to infer current rates of hybridization and to determine patterns of introgression across the genome and geography. Statistical Species Delimitation One of the most common ways to validate cryptic lineages is through multilocus coalescent-based (MSC) approaches (Fujita et al. 2012). Accordingly, and as recommended by Rannala (2015), we used two MSC approaches (bpp v3.3a and stacey v1.2.4) to test species boundaries across these groups (Yang and Rannala 2014; Jones 2017). Using the same filtering as in our starbeast2 analyses, we selected three random samples of 100 loci per species group and generated input files for bpp and stacey. We used the species tree inferred from the starbeast2 analyses (see Phylogenetic Analyses) as the guide tree for BPP. For the ‘L. coggeri’ group, we used two topologies of the species tree that reflected uncertainty in the placement of L. coggeri EU (Supplementary Fig. S1 available on Dryad). We then ran bpp for 500,000 generations across three sets of priors to ensure our results were robust to prior specification. These priors were: 1) $$\rm{\theta} \sim$$ (2, 2000), $$\tau \sim$$ (2, 2000), 2) $$\theta \sim$$ (1, 10), $$\tau \sim$$ (1, 10), and 3) $$\theta \sim$$ (1, 10), $$\tau \sim$$(2, 2000). We ensured that we had 20–80% acceptance rate; having too high or too low of acceptance rates can affect results. We ran stacey for 1e7 generations. Each locus was set to have its own clock rate and own substitution model under a HKY model. Priors were: collapse height $$=$$ 0.001, growth rate $$\sim$$ lognormal(mean $$=$$ 5, sd $$=$$ 2); collapse weight $$\sim$$ uniform (0,1); population prior scale $$\sim$$ lognormal (mean $$= -7$$, sd $$=$$ 2), and relative death rate $$\sim \beta (\alpha = 1, \beta = 8$$). Species delimitations were determined using speciesDA with a burn-in of 10% and a collapse height of 0.0001 (Jones 2017). Analyses showed that results were robust across collapse heights from 0.0001 to 0.001. Morphological Data and Analysis We assessed morphological differences in scale and body measurement traits between the major genetic lineages in each of the three species groups in the Wet Tropics region. For the ‘C. rubrigularis’ species group, we excluded C. wundalthini and C. rhomboidalis because they differ in breeding colors and to some degree body size and shape (Dolman 2008; Hoskin 2014). Further, we combined individuals from L. coggeri N and C, from L. robertsi CU and TU, and L. robertsi BK and BFAU because our analyses of the genomic evidence suggested each of these lineage pairs are not distinct and should be collapsed (see the Discussion). A single investigator took scale counts and body measurements on an average of 26 specimens per lineage (Supplementary Appendix 1, Table S4 available on Dryad). The traits measured summarize morphological characteristics that affect how lizards function in their environment (e.g., relative limb length affects locomotion; Losos 2011) and are standardly used to delimit skink species (Ingram and Covacevich 1989; Ingram 1991; Greer 1997). The scale traits counted were the number of supraciliaries, infralabials, supralabials, midbody scale rows, paravertebrals, and lamellae under the fourth toe (Hoskin 2014). We found no to little variation for these scale traits within species groups and thus performed no further analyses (Supplementary Table A4 available on Dryad). We measured five aspects of body size and shape: total head and body length (snout to vent length, SVL), distance between the front and hindlimbs (axilla to groin length, AG), length of the hindlimb (L2), head length (HL), and head width (HW) (Hoskin 2014). Only adults were measured, defined as individuals with SVL greater than 48 mm, 44 mm, and 32 mm in C. rubrigularis, L. robertsi, and L. coggeri, respectively. For each species group, we used multivariate analyses to determine if lineages differed significantly in size or shape. We first used a principal component analysis (PCA) to summarize across all five body measurements (Jolliffe 2002). We then tested if body size (PC1) and body shape (PCs 2–5) varied significantly across lineages. All analyses were nested by region within lineage to account for among-region variation within lineages. “Region” represents different mountain ranges, tablelands, and lowland areas, and roughly matched the subregions defined for the Wet Tropics (Bell et al. 2010). For body size, we used a nested univariate analysis of variance (ANOVAs) of PC1 between lineages within each species group. For body shape, we used a nested multivariate analysis of variance (MANOVA) of PCs 2–5. Significance was assessed using Roy”s Greatest Root. For those species groups where we included more than two lineages, we performed nested univariate (size) or multivariate (shape) planned contrasts to test which lineages differed significantly. All analyses were conducted in sas v9.2. Results Analysis of Genetic Data Our exome capture approach worked well across all lineages. On average, we recovered an average of 2.29 Mb per individual across 2668 loci at an average coverage of 112$$\times$$ (Supplementary Table S2 available on Dryad). The inferred topology is well-supported and is consistent with previous phylogenetic hypotheses based on mtDNA data (Fig. 2, Supplementary Fig. S2 available on Dryad; Bell et al. 2010) and other phylogenomic analyses (Bragg et al. 2018). Branching times and topologies were quantitatively and qualitatively similar across replicate analyses of starbeast2 (Supplementary Fig. S1 available on Dryad). The branching times between these lineages all occur within a narrow range and highlight that two lineages of C. rubrigularis as currently recognized are polyphyletic. Further, the splitting patterns generally agree with the biogeographic relationships between species. For example, L. robertsi BFAU is sister to L. robertsi BK, and the two lineages occur as adjacent montane isolates (Figs. 1 and 2). The exception to this congruence is a leapfrog distribution of L. coggeri EU, which is sister to L. coggeri N and C rather than the geographically adjacent L. coggeri S (Figs. 1 and 2). We inferred a 4$$\times$$ to 9$$\times$$ range of genetic divergences between lineages within nominal species, with nuclear $$d_{xy}$$ at silent sites ranging from 0.5% to 1.95% and nuclear $$d_{a}$$ ranging from 0.18% to 1.68% (Fig. 3, Supplementary Table S5 available on Dryad). Lineages associated with small ranges such as C. wundalthini and the montane lineages of L. robertsi showed lower levels of within-population diversity. Measures of mtDNA and nuclear $$d_{xy}$$ and $$d_{a}$$ were correlated ($$d_{xy}$$; $$r = 0.61$$, $$P$$-value $$=$$ 0.016; $$d_{a}$$; $$r = 0.57$$, $$P$$-value $$=$$ 0.025; Fig. 3A, B). As for inferred branching times, divergences between some cryptic lineages were significantly greater than those between nominal species (Fig. 3C). Figure 3. View largeDownload slide Patterns of divergence between pairwise lineage-comparisons for (A) $$d_{xy}$$ for silent sites in mitochondrial DNA (mtDNA) and nuclear DNA (nDNA), (B) $$d_{a}$$ at silent sites in mtDNA and nDNA, and (C) branching times in units of substitutions per site per million years and $$d_{a}$$ at mtDNA; branching times were inferred from the tree depicted in Fig. 1. Each pairwise comparison is coded as either being between 1) recognized: two lineages that were already recognized at the species-level, 2) elevated: two lineages that were elevated in the current study, or 3) population: lineages for which there is insufficient evidence to elevate them to species. Arrows identify the three pairwise comparisons for which we have data on reproductive isolation from contact zone studies (Fig. 4), the gray line indicates the transition point at which we first recover evidence for isolation between lineages (i.e., between Carlia rubrigularis N and S). Many of the lineages that we propose to elevate are more diverged than recognised species. Figure 3. View largeDownload slide Patterns of divergence between pairwise lineage-comparisons for (A) $$d_{xy}$$ for silent sites in mitochondrial DNA (mtDNA) and nuclear DNA (nDNA), (B) $$d_{a}$$ at silent sites in mtDNA and nDNA, and (C) branching times in units of substitutions per site per million years and $$d_{a}$$ at mtDNA; branching times were inferred from the tree depicted in Fig. 1. Each pairwise comparison is coded as either being between 1) recognized: two lineages that were already recognized at the species-level, 2) elevated: two lineages that were elevated in the current study, or 3) population: lineages for which there is insufficient evidence to elevate them to species. Arrows identify the three pairwise comparisons for which we have data on reproductive isolation from contact zone studies (Fig. 4), the gray line indicates the transition point at which we first recover evidence for isolation between lineages (i.e., between Carlia rubrigularis N and S). Many of the lineages that we propose to elevate are more diverged than recognised species. Figure 4. View largeDownload slide Evidence for rates of hybridization and introgression at three contact zones between lineages included in this analysis: Carlia rubrigularis N and S, Lampropholis coggeri C and S, and L. coggeri N and C. Contacts are listed in the legend in order of least to most divergent. A) Percent of individuals in the center of contact zones that were identified as hybrid. A hybrid individual was defined as individuals that had $$\geqslant$$10% membership in both parental species as determined by STRUCTURE. B) Distribution of cline widths in the contact zone across an average of 9.5K clines. C) The extent of linkage disequilibrium in each contact zone. Moran”s I measures the autocorrelation in cline widths across the genome, which serves as a proxy for linkage disequilibrium. These genetic estimates of reproductive isolation show evidence for selection against hybrids in C. rubrigularis N and S and L. coggeri C and S. Figure 4. View largeDownload slide Evidence for rates of hybridization and introgression at three contact zones between lineages included in this analysis: Carlia rubrigularis N and S, Lampropholis coggeri C and S, and L. coggeri N and C. Contacts are listed in the legend in order of least to most divergent. A) Percent of individuals in the center of contact zones that were identified as hybrid. A hybrid individual was defined as individuals that had $$\geqslant$$10% membership in both parental species as determined by STRUCTURE. B) Distribution of cline widths in the contact zone across an average of 9.5K clines. C) The extent of linkage disequilibrium in each contact zone. Moran”s I measures the autocorrelation in cline widths across the genome, which serves as a proxy for linkage disequilibrium. These genetic estimates of reproductive isolation show evidence for selection against hybrids in C. rubrigularis N and S and L. coggeri C and S. Our D-statistic tests for historical introgression among non-sister lineages recovered four likely cases of historical introgression between lineage-pairs: C. rubrigularis N and S; C. rubrigularis N and C. rhomboidalis; L. coggeri S and EU; and L. coggeri C and S (Fig. 2, Supplementary Table S6 available on Dryad). No signature of historical introgression was recovered for strongly allopatric populations—e.g., among montane isolates of the ‘L. robertsi’ species group, and C. wundalthini vs. C. rubrigularis. Our previous results from analyses of hybridization and introgression at contact zones show that C. rubrigularis N and S and L. coggeri C and S exhibit 1) a moderate proportion of hybrids in the center of the contact zones, 2) narrow cline widths across the genome, and 3) auto-correlation in cline widths across physical distances, all indicative of extensive disequilibrium in hybrids and substantial RI (Fig. 4). The less divergent lineage-pair, L. coggeri N and C, shows none of the same patterns, with evidence of extensive introgression across the genome and geography. Statistical species delimitation supported all lineages as species. bpp returned $$\geqslant$$95% probability for a speciation event at all nodes in our guide trees, and stacey inferred each species as a unique cluster (Fig. 2). This result was robust to priors and, for the “L. coggeri” group, uncertainty in the topology. Analysis of Morphological Data For all three species groups, PCAs resulted in a PC1 that accounted for most of the variation (67–84%). It was loaded highly and positively by all body measurements and indicated body size (Supplementary Tables S7–S9 available on Dryad). The remaining four PCs in each species group represented variation in body shape (Supplementary Tables S7–S9 available on Dryad). PC2 accounted for 10–16% of the variation in the species groups and was loaded most heavily by relative AG length (positive) in ‘C. rubrigularis’ and both AG length (positive) and relative L2 length (negative) in ‘L. robertsi’ and ‘L. coggeri’ (Supplementary Tables S7–S9 available on Dryad). The three major lineages in the ‘L. coggeri’ group show no differences in body size (‘lineage’ effect on PC1, $$F_{\mathrm{2,58}} = 2.46$$, $$P = 0.09$$), but they do differ in body shape (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{8,10}} = 3.46$$, $$P = 0.04$$) (Fig. 5B, Table 1). This significant variance was driven by CV1 (Roy”s Greatest Root $$F_{\mathrm{4,6}} = 6.05$$, $$P = 0.03$$), which is loaded most heavily by PC2 (0.489; Supplementary Table S10 available on Dryad). Multivariate contrasts revealed that only L. coggeri EU and S differ significantly in shape ($$F_{\mathrm{4,5}} = 5.00$$, $$P = 0.05$$) (Table 1). Neither L. coggeri EU and N/C ($$F_{\mathrm{4,5}} = 3.35$$, $$P = 0.11$$) nor L. coggeri N/C and S ($$F_{\mathrm{4,5}} = 2.56$$, $$P = 0.17$$) differed in shape. Therefore, the only detectible difference was that L. coggeri EU has a relatively longer body and shorter legs than the geographically adjacent, but distantly related, L. coggeri S. Table 1. Testing for morphological differences in body size and shape between the major lineages in each species group $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * Body size was tested using nested ANOVA. Body shape was tested using nested MANOVA (‘lineage’ overall effect and planned contrasts), and using the Wilks’ Lambda as the $$F$$-statistic. Table 1. Testing for morphological differences in body size and shape between the major lineages in each species group $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * Body size was tested using nested ANOVA. Body shape was tested using nested MANOVA (‘lineage’ overall effect and planned contrasts), and using the Wilks’ Lambda as the $$F$$-statistic. Figure 5. View largeDownload slide Scatter plots of PC1 (representing body size) and PC2 (the primary axis of body shape) in each of the species groups: (A) ‘Carlia rubrigularis’ group, (B) ‘Lampropholis coggeri’ group, (C) ‘Lampropholis robertsi’ group. See Supplementary Tables S7–S9 available on Dryad for details on loadings of PC2. Morphological divergence among the lineages in each species group is limited, with the exception of some shape divergence in L. coggeri EU. Figure 5. View largeDownload slide Scatter plots of PC1 (representing body size) and PC2 (the primary axis of body shape) in each of the species groups: (A) ‘Carlia rubrigularis’ group, (B) ‘Lampropholis coggeri’ group, (C) ‘Lampropholis robertsi’ group. See Supplementary Tables S7–S9 available on Dryad for details on loadings of PC2. Morphological divergence among the lineages in each species group is limited, with the exception of some shape divergence in L. coggeri EU. No morphological differences were detected between the N and S lineages of ‘C. rubrigularis’, for either body size (‘lineage’ effect on PC1, $$F_{\mathrm{1,50}} = 1.17$$, $$P = 0.29$$) or body shape (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{4,5}} = 0.45$$, $$P = 0.77$$) (Fig. 5A, Table 1). Similarly, no significant differences were detected between the TU/CU and BK/BFAU lineages of ‘L. robertsi’. Body size was marginally non-significant (‘lineage’ effect on PC1, $$F_{\mathrm{1,28}} = 3.89$$, $$P = 0.06$$), and body shape did not differ (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{3,1}} = 2.68$$, $$P = 0.42$$) (Fig. 5C, Table 1). Discussion Delimitation of Cryptic Species Initial phylogeographic explorations based on mtDNA revealed that each species group contained at least four to five lineages, most of which were deeply divergent. Subsequent sequencing of five to ten nuclear loci confirmed that these phylogeographic lineages were also diverged at the nuclear genome (Dolman and Moritz 2006; Bell et al. 2010). Now, genetic data based on over 2500 exons confirmed that these lineages exhibit genetic divergences of substantial but varying depths. These genetic divergences all fall within the range that comparative data suggest spans the transition from populations to isolated species (Roux et al. 2016). Although some of the lineages are far more divergent than some already recognized species, and although we focused on morphological traits standardly used in lizard taxonomy and eco-evolutionary studies (Ingram 1991; Losos 2011; Hoskin 2014), we found little or no morphological divergence between the major lineages within each of the three species groups. Accordingly, these can mostly be regarded as truly cryptic, rather than pseudo-cryptic, species. Given morphologically cryptic lineages that span a range of divergences, and all of which are delimited using coalescent methods, we are thus faced with the two challenges of species delimitation: how to determine how much genetic divergence is sufficient when divergences are arrayed on a continuum; and how to reconcile when genetic and phenotypic data give conflicting perspectives. As a first step, we can directly assess levels of isolation between these lineages because three of the lineage-pairs in these groups meet in narrow zones of parapatry (Fig. 4). Through these fine-scale contact zone analyses of isolation, we have two major findings. First, we find that, like genetic divergence, RI exists on a continuum, with lineages exhibiting varying degrees of isolation (Fig. 4A, B). For this set of lineages, divergence and isolation appear to scale, although non-linearly. The average cline width at the L. coggeri C and S hybrid zone zone is about 5.5$$\times$$ less than that at the L. coggeri N and C contact. Yet, L. coggeri C and S is only 1.6$$\times$$ more genetically divergent than L. coggeri N and C. Theory does not predict a linear scaling; the cline width of a locus is proportional to the inverse square root of selection on that locus (Barton and Gale 1993). Thus, as selection on a locus increases, cline width can sharpen narrowly and quickly, as seen here. Further evidence of this non-linear accumulation of RI can be seen in patterns of linkage disequilibrium in these contact zones. Data from L. coggeri N and C show no evidence for disequilibrium at introgressing sites, whereas C. rubrigularis N and S and L. coggeri C and S exhibit extensive disequilibrium extending a few kilobases (Fig. 4C). These results confirm theoretical expectations that lineages can quickly transition from acting as populations (i.e., L. coggeri N and C) to acting as genomically isolated species (i.e., L. coggeri C and S and C. rubrigularis N and S) (Turner 1967; Barton 1983). Second, these data show that, despite being nearly identical morphologically, the more genetically divergent lineages have substantial RI. At least for C. rubrigularis N and S, these lineages are not isolated by premating isolation (Dolman 2008), but rather by post-mating selection against hybrids (Phillips et al. 2004). Based on estimates of dispersal length and cline width of the hybrid index, selection against hybrids is strong. Hybrids between C. rubrigularis N and S and L. coggeri C and S are estimated to be 50–70% and 10–65% less fit than their parents, respectively (Phillips et al. 2004; Singhal and Moritz 2012). Estimated selection on hybrids between L. coggeri N and C, on the other hand, is negligible. The selection against hybrids seen in C. rubrigularis N and S and L. coggeri C and S is comparable (if not greater) than that seen between morphologically distinct hybridizing taxon-pairs. (Barton and Gale 1993; Singhal and Moritz 2012). Such strong selection suggests that these lineages will remain evolutionary distinct in the future, despite the high potential for gene flow. As such, we propose to identify C. rubrigularis N and S as separate species, and likewise for L. coggeri S and L. coggeri C. Because we found no evidence for RI between L. coggeri N and C, we retain them as distinct populations within one species (L. coggeri N/C). However, how should we diagnose those lineages for which we cannot indirectly or directly assay RI? For example, L. coggeri EU is geographically isolated from L. coggeri N/C and S, and the lineages in the “L. robertsi” group are isolated on different mountaintops. These lineages do not meet in contact zones, and because of both practical and ethical reasons, cannot be easily kept in the laboratory for experimental trials. Instead, we extrapolate our estimates of RI from the three lineage-pairs that do meet in parapatry (Fig. 4) to the species group as a whole. This extrapolation assumes that the tempo and mode at which RI evolves is similar across this clade. The few comparative data on the rate at which RI evolves suggests that it can vary across broad clades (Rabosky and Matute 2013). However, for a clade like this, which consists of broadly related, morphologically and ecologically-similar lizards found in a similar biogeographic context, we suspect there is likely to be less variation. Indeed, RI and divergence time correlate closely across five sister-species comparisons in Carlia, Lampropholis, and a closely-related genus, Saproscincus (Singhal and Moritz 2013; Singhal and Bi 2017). Further, and importantly, these lineages likely resulted from similar speciation processes—i.e., these deep, cryptic lineages evolved due to very long periods of isolation in environmentally similar refugia (see below). Therefore, we believe we can sensibly extrapolate our results across other cryptic congeneric lineages that diverged under similar processes. To extrapolate, we use the divergence between C. rubrigularis N and S as our cutoff because they are the youngest lineage-pair for which we have solid evidence of RI. Divergence estimates are highly correlated across the three metrics for genetic divergence ($$r \sim 0.7-9$$; Fig. 3). Still, we take a conservative approach and only elevate those lineages that show greater divergence than what is seen for C. rubrigularis N and S across all metrics (with one exception, below). Notably, this cutoff is greater than that seen among several comparisons between nominal taxa (Fig. 3). Divergences between L. robertsi CU and TU, L. robertsi BK and BFAU, and C. rhomboidalis N and S all fall below this cutoff (Supplementary Table S5 available on Dryad), so we recognize these as phylogeographic lineages within species. However, the divergence between L. robertsi BK/BFAU and L. robertsi CU/TU is greater than this cutoff. Accordingly, we propose to diagnose the BK/BFAU and CU/TU allopatric lineages as species; this deep divergence suggests they are likely to exhibit RI should they ever come into contact. We use these groupings for morphological comparisons, as well. The case of L. coggeri EU is more ambiguous, when considering genetic data alone. Lampropholis coggeri EU is most closely related to L. coggeri N/C, and divergence falls just below our proposed cutoff for raw divergence and branching time (Fig. 3A, C) but just above the cutoff for net divergence (Fig. 3B). However, L. coggeri EU sits isolated off the far south end of the Wet Tropics, geographically closest to L. coggeri S. Therefore, it is much more likely to interact with L. coggeri S in future, with which it shows much greater divergence (Fig. 1, Fig. 3; Supplementary Table S5 available on Dryad). Further, alone among the comparisons made here, L. coggeri EU is morphologically (and perhaps ecologically; see below) distinct from L. coggeri S (Table 1). Because the L. coggeri S and EU comparison is more salient than the L. coggeri C and EU comparison, we further propose to elevate L. coggeri EU as a separate species. Taxonomic Outcomes The formal taxonomic revisions of these lineages are presented in Appendix 1 available on Dryad. To summarize, we revise the ‘C. rubrigularis’ group to retain C. wundalthini and C. rhomboidalis, retain C. rubrigularis S as C. rubrigularis, and elevate C. rubrigularis N to C. crypta sp. nov. For the ‘L. coggeri’ group, we retain L. coggeri N/C as L. coggeri, elevate L. coggeri S to L. similis sp. nov., and elevate L. coggeri EU to L. elliotensis sp. nov. For the ‘L. robertsi’ group, we retain L. robertsi CU/TU as L. robertsi and elevate L. robertsi BFAU/BK to L. bellendenkerensis sp. nov. The key components of the new species descriptions are as follows. Carlia crypta sp. nov. Holotype: QM J75457 Diagnosis: Carlia crypta sp. nov. is distinguished from all other Carlia spp., except members of the ‘C. rhomboidalis’ group, in possessing an interparietal fused to the frontoparietals. As with C. rubrigularis, adult males possess a red throat. It is reliably distinguished from this species by four nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in three amino acid differences (Table A3). Zoobank ID: 3B116172-035D-41EF-958E-F0A35BB15F0D Lampropholis similis sp. nov. Holotype: QM J91380 Diagnosis: Lampropholis similis sp. nov. is a small, dark-sided rainforest skink with pentadactyl limbs (overlapping or very narrowly separated when adpressed) and a movable lower eyelid containing a transparent disc. It is reliably distinguished from its sibling species (L. coggeri and L. elliotensis sp. nov.) by 17 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in 15 amino acid differences (Table A1). Zoobank ID: 23E36086-C261-4D2A-958B-547A17239529 Lampropholis elliotensis sp. nov. Holotype: QM J91382 Diagnosis: Lampropholis elliotensis sp. nov. is a small, dark-sided rainforest skink with pentadactyl limbs (usually separated by several scales rows when adpressed) and a movable lower eyelid containing a transparent disc. It is reliably distinguished from its sibling species (L. similis sp. nov. and L. coggeri) by 17 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in 15 amino acid differences among the species (Table A1). Zoobank ID: 6CE54452-528E-4B3E-9655-AC2262A26849 Lampropholis bellendenkerensis sp. nov. Holotype: QM J39855 Diagnosis: A large Lampropholis with dark flanks and prominent spotting on the posterior ventral surfaces, a row of dark edged pale spots on underside of tail. This species is reliably distinguished from its closest congener (L. robertsi) by 13 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in nine amino acid differences between the species (Table A2). Zoobank ID: C721A396-F4A2-4205-A7FE-ED8FC8F 90FAD Speciation Processes and Cryptic Species Extrapolating a divergence cutoff as done here works, in part, because the Wet Tropics lineages likely share similar divergence histories and processes. These taxa diverged through extended periods of isolation in environmentally-similar, climatically-stable rainforest refugia, resulting in genetic divergence but eco-morphological conservatism (Graham et al. 2006; Moritz et al. 2009). The widespread Wet Tropics lineages (C. rubrigularis N and S; L. coggeri N/C and S) occupy similar habitats, and the lineages of ‘L. robertsi’ occupy mountaintop habitats that appear broadly similar. The morphological stasis in these cryptic species therefore likely reflects a lack of divergent selection across similar environments (Moritz et al. 2009; Hoskin et al. 2011). Only those lineages peripheral to the main block of the Wet Tropics rainforests exhibit morphological divergence. For example, L. coggeri EU differs subtly in body shape from other species in the group; it also occupies a subset of the habitats occupied by the other lineages in the ‘L. coggeri’ group. Whereas the other lineages are found across a broad range of rainforest types and habitat edges, L. coggeri EU is found only in mesic pockets of rocky, upland rainforest. Greater morphological divergence is seen in the two lineages in the ‘C. rubrigularis’ group found outside the Wet Tropics region: C. wundalthini and C. rhomboidalis. Despite being of similar age as other lineages (Fig. 2), these species are distinct for male breeding coloration and some aspects of body size and shape. As seen across other Wet Tropics taxa, phenotypic divergence is only recovered in association with environmentally-driven selection – e.g., across ecotones (Schneider and Moritz 1999), in peripheral isolates (Hoskin et al. 2011), or reinforcing selection in hybrid zones (Hoskin et al. 2005). These cryptic species thus underline the importance of “non-ecological” mechanisms in speciation (Schluter 2001), in particular, mutation-order speciation (Mani and Clarke 1990; Nosil and Flaxman 2011). This work also helps further outline the distinction between historical and current patterns of hybridization and introgression (Edwards et al. 2016). The biological species concept requires current introgression to be limited, but both theoretical models and empirical data show that species can diverge and remain distinct in the presence of historical gene flow (Nosil 2008; Pinho and Hey 2010). Similarly, our tests for introgression suggest there has been introgression between at least four lineage-pairs (Fig. 2), all of which are either previously-recognized or now elevated nominal species. Yet, for at least two of these lineage-pairs (C. rubrigularis N and S and L. coggeri C and S), our analysis of current patterns shows that hybridization occurs but is geographically highly-restricted (Fig. 4). This distinction between historical and current patterns illustrates the complicated relationship between gene flow and species borders. The Practicality of Delimiting Cryptic Species Genetic divergences for almost all our lineage comparisons fall in the so-called “gray zone” of speciation, in which lineages transition from behaving as populations to species (Roux et al. 2016). Defined by net silent divergence ($$d_{a})$$ at coding nuclear genes, this gray zone spans divergences from 0.5% to 2%. Given that many cryptic lineages originated during glacial cycles over the last few million years (Hewitt 2000), many of them should fall within this four-fold range of divergence. This underlines the challenge in delimiting cryptic lineages—many of them have a biogeographic and divergence history that places them in an ambiguous zone of divergence, where lineages are as likely to merge or remain distinct upon secondary contact. Given this ambiguity, identifying strong phylogeographic structure within species should be just the first step in diagnosing species boundaries across cryptic boundaries (Fig. 6, Table 2). Additional validation is required, which we loosely group into four categories: 1) statistical species delimitation, 2) post hoc discovery of phenotypic differences that delimit lineages (i.e., integrative taxonomy; Padial et al. 2010), 3) indirect or direct estimates of evolutionary isolation between lineages, or 4) calibration-based approaches. Herein, we outline these approaches briefly and explain their conceptual and practical benefits and limitations. Note that after applying any of these approaches, researchers must still formally revise the taxonomy for these validations to be recognized. However, all too often, these taxonomic revisions are not done (Carstens et al. 2013; Pante et al. 2014). Table 2. A survey of approaches that can be used to validate putative cryptic lineages and the benefits and limitations of each approach Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). We highlight, which approaches were used in this study and identify other examples from the broader literature; this list of examples is meant to be illustrative not exhaustive. Many of the studies that validated putative cryptic species did not also conduct a formal taxonomic revision. Table 2. A survey of approaches that can be used to validate putative cryptic lineages and the benefits and limitations of each approach Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). We highlight, which approaches were used in this study and identify other examples from the broader literature; this list of examples is meant to be illustrative not exhaustive. Many of the studies that validated putative cryptic species did not also conduct a formal taxonomic revision. Figure 6. View largeDownload slide A flowchart outlining a possible research approach to validating cryptic lineages. Statistical species delimitation is often used to both define putative lineages and to validate them. Width and contrast of validation arrows indicate how robust we believe these approaches are. Benefits and limitations of each approach are further described in Table 2. Figure 6. View largeDownload slide A flowchart outlining a possible research approach to validating cryptic lineages. Statistical species delimitation is often used to both define putative lineages and to validate them. Width and contrast of validation arrows indicate how robust we believe these approaches are. Benefits and limitations of each approach are further described in Table 2. First, perhaps the most common approach currently used is statistical species delimitation, which applies coalescent-based methods to determine which lineages are genetically unique (Ence and Carstens 2011; Yang and Rannala 2014). Often, statistical approaches are used as an early step to more quantitatively assess visual clusters, and in other studies, they are used as a final stage of analysis to diagnose species (Fig. 6). Across our three species groups, statistical species delimitation diagnosed all putative lineages as species (Fig. 2). Yet, this statistical diagnosis contrasts to our biological understanding of species boundaries. For example, although L. coggeri N and C are statistically distinct, introgression between them is widespread genomically and geographically (Fig. 4). This disconnect reflects the emerging consensus that, although statistical species delimitation methods can robustly identify populations, these populations are not always equivalent to species (Rosenblum et al. 2012; Dynesius and Jansson 2014; Sukumaran and Knowles 2017). In other words, the genetic distinctiveness of a population does not necessarily confer robust evolutionary distinctiveness as envisaged under the BSC. Thus, we take a deliberately conservative approach and refrain from elevating morphologically cryptic lineages whose sole support is from statistical analyses of genetic data (Oliver et al. 2015). However, statistical approaches to species delimitation are the most efficient and flexible method across taxonomic groups and are likely to remain an attractive option for many study systems, even if these approaches alone provide insufficient evidence to denote robust taxa. Second, while mostly not true in the present study (Fig. 5), researchers often discover post-hoc phenotypic differences after further investigation of putative cryptic lineages, leading to so-called pseudo-cryptic species (Knowlton 1993). For example, cryptic lineages might vary in traits that facilitate RI between lineages (e.g., mating calls, (Barber 1951)) or ecological co-existence (e.g., divergence in life-history; Leys et al. 2017). Often, however, these phenotypic differences might only relate trivially to how distinctive a lineage is—i.e., minor differences in scale counts. In such cases, phenotypic differentiation alone might not necessarily confer evolutionary distinctiveness, and further validation would be required. Third, researchers can assay strength of RI between lineages either through indirect studies of contact zones, observation of sympatry among cryptic lineages, or laboratory- or field-based tests of mate choice and hybrid fitness (Blair 1972; Hoskin et al. 2005; Dolman 2008). As shown in this work (Fig. 4), these approaches offer detailed data on how likely lineages are to remain distinct if and when they overlap with their congenerics. However, these approaches are often quite practically limited—generating these data can require contiguous ranges, high population densities, organisms amenable to experimentation, and substantial investment in both time and money. These practical limitations surfaced in the present work; we were unable to test for RI between allopatric lineages, nor could we bring them into a laboratory setting. Barring this, alternative approaches could be used to identify cases of abrupt genetic boundaries across dense sampling or test for geographically-extensive introgression across lineage boundaries using large numbers of markers (Melville et al. 2017). For those lineages not amenable to indirect or direct testing for RI, we used a fourth general approach: using calibrations to determine how much divergence is sufficient to elevate a species. In the DNA barcoding literature, these calibrations are typically informed by patterns of divergence among nominal taxa, although they are used across broad swaths of the tree of life—i.e., all birds or all butterflies (Hebert et al. 2004b; Janzen et al. 2005). Another option is to use calibrations informed by data on the tempo at which RI evolves in closely-related taxa, as done in this study. While this approach still requires the expensive and time-consuming generation of data on isolation between lineages, the cutoff is more principled than barcode gaps. Unlike barcode gaps, which are applied widely across taxonomic groups, our cutoff is based on the observed isolation between lineages within a given clade, which likely share a common mode of lineage divergence and speciation. Such an approach allows us to tackle cryptic diversity while reflecting the variable nature of the speciation process. Moving forward, we will explore whether this cutoff could contribute to diagnosing species boundaries among phylogeographic lineages of other species of Carlia (Potter et al. 2016). Across all these approaches, a potential flaw is insufficient sampling of species ranges. In particular, by sampling two ends of an array of populations, we can infer distinction between lineages where gene flow is actually continuous throughout (Pante et al. 2015). Because all approaches to cryptic species validation either have conceptual or practical weaknesses (Table 2), the ideal approach will likely vary across taxonomic groups. For example, in the Wet Tropics, phylogeographic studies have recovered deep splits in other taxa outside of lizards, including frogs and mammals (Moritz et al. 2009). For the frog lineages, many of which occur in dense numbers, meet in narrow contact zones, and are amenable to breeding experiments, data on hybridization patterns and mate choice have helped validate putative cryptic lineages and led to the formal revision of lineages (Hoskin et al. 2005; Hoskin 2007). However, for mammals, density is too low to allow indirect or direct tests of isolation, so other approaches will be required. Importantly, the framework we applied here—generating initial descriptions of within-species phylogeographic diversity (Dolman and Moritz 2006; Bell et al. 2010), confirming these patterns with multilocus data sets, and then inferring fine-scale patterns of current RI (Phillips et al. 2004; Singhal and Moritz 2012, 2013; Singhal and Bi 2017)—represents a significant investment in time and money. For many systems, this approach is simply not tenable, and further, often too slow where decisions on species limits have immediate conservation consequences (Hedin 2015). That said, given the dangers of ‘taxonomic overinflation’ (Isaac et al. 2004; Frankham et al. 2012), we advocate that researchers validate putative cryptic lineages by both considering statistical delimitation approaches and other data on the biological reality of lineages, whether that be direct or indirect evidence for isolation. Conclusion Cryptic species challenge traditional notions of species, because the discrepancy between morphological and genetic axes of divergence can make them hard to categorize. Yet, other data suggest that many cryptic species are phenotypically divergent, but on axes of variation that are harder to measure (e.g., mating pheromones in lizards). In cases where we cannot identify phenotypic differences, like the lizards of the Wet Tropics, we can test the validity of these lineages through other means, such as looking at interactions between cryptic lineages in parapatry. Often, these richer, more integrative datasets complement genetic data and show that cryptic lineages are independently evolving units. However, as we also see in these taxa, despite marked genetic differentiation, some cryptic lineages might just be ephemera, destined to be lost to hybridization with congeneric lineages, if they meet in the future. These data remind us that species boundaries are hypotheses (Fujita et al. 2012), our best estimate of the fate of these lineages and a recognition of the ever-evolving nature of species (Darwin 1859; Mallet 1995). Acknowledgments We thank J. Bragg and also staff at UMich Arc TS Flux for advice and logistical support. In addition, we gratefully acknowledge the members of the C. Moritz & S. Williams research groups, both past and present, who have contributed samples, environmental data, and their expertise to the Wet Tropics research program. Funding Funding for this project came from a National Science Foundation Post-doctoral Research Fellowship in Biology (NSF #1519732 to SS), grants from the Australian Research Council (ARC), the Australian Biological Resources Study (ABRS) and the National Science Foundation to CM, and grants from the ARC and ABRS to CH. Data Accessibility Raw sequencing data: PRJNA448788 Pseudo-reference files and variant sets: DataDryad accession 10.5061/dryad.g7v1b. Code used to generate data: https://github.com/singhal/AWT_delimit and https://github.com/singhal/SqCL. Specimens used for morphological analyses & species descriptions: all are deposited at Queensland Museum for research use; accession numbers available in Supplementary Appendix 1 available on Dryad. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.g7v1b. References Adams M. , Raadik T.A. , Burridge C.P. , Georges A. 2014 . Global biodiversity assessment and hyper-cryptic species complexes: more than one species of elephant in the room? Syst. Biol. 63 : 518 – 533 . Google Scholar CrossRef Search ADS PubMed Amato A. , Kooistra W.H. , Ghiron J.H.L. , Mann D.G. , Pröschold T. , Montresor M. 2007 . Reproductive isolation among sympatric cryptic species in marine diatoms . Protist 158 : 193 – 207 . Google Scholar CrossRef Search ADS PubMed Bálint M. , Domisch S. , Engelhardt C. , Haase P. , Lehrian S. , Sauer J. , Theissinger K. , Pauls S. , Nowak C. 2011 . Cryptic biodiversity loss linked to global climate change . Nat. Clim. Change 1 : 313 – 318 . Google Scholar CrossRef Search ADS Barber H.S. 1951 . North American fireflies of the genus Photuris . Smithson. Misc. Collect. 117 : 1 – 58 . Barton N. 1983 . Multilocus clines . Evolution 37 : 454 – 471 . Google Scholar CrossRef Search ADS PubMed Barton N.H. , Gale K.S. 1993 . Genetic analysis of hybrid zones. In: Hybrid zones and the evolutionary process ., (ed. Harrison RG ), pp. 13 – 45 . New York : Oxford University Press . Bell R.C. , Parra J.L. , Tonione M. , Hoskin C.J. , Mackenzie J.B. , Williams S.E. , Moritz C. 2010 . Patterns of persistence and isolation indicate resilience to climate change in montane rainforest lizards . Mol. Ecol. 19 : 2531 – 2544 . Google Scholar PubMed Bickford D. , Lohman D.J. , Sodhi N.S. , Ng P.K. , Meier R. , Winker K. , Ingram K.K. , Das I. 2007 . Cryptic species as a window on diversity and conservation . Trends Ecol. Evol. 22 : 148 – 155 . Google Scholar CrossRef Search ADS PubMed Blair W.F. 1972 . Evolution in the genus Bufo . Austin (TX) : University of Texas Press . Bolger A.M. , Lohse M. , Usadel B. 2014 . Trimmomatic: a flexible trimmer for Illumina sequence data . Bioinformatics 30 : 2114 – 2120 . Google Scholar CrossRef Search ADS PubMed Bolnick D.I. , Near T.J. , Wainwright P.C. 2006 . Body size divergence promotes post-zygotic reproductive isolation in centrarchids . Evol. Ecol. Res. 8 : 903 – 913 . Bragg J.G. , Potter S. , Afonso Silva A.C. , Hoskin C.J. , Bai B.Y.H. , Moritz C. 2018 . Phylogenomics of a rapid radiation: the Australian rainbow skinks . BMC Evol. Biol. 18 : 15 . Google Scholar CrossRef Search ADS PubMed Bragg J.G. , Potter S. , Bi K. , Moritz C. 2016 . Exon capture phylogenomics: efficacy across scales of divergence . Mol. Ecol. Resour. 16 : 1059 – 1068 . Google Scholar CrossRef Search ADS PubMed Brown J.M. , Hedtke S.M. , Lemmon A.R. , Lemmon E.M. 2009 . When trees grow too long: investigating the causes of highly inaccurate Bayesian branch-length estimates . Syst. Biol. 59 : 145 – 161 . Google Scholar CrossRef Search ADS PubMed Carstens B.C. , Pelletier T.A. , Reid N.M. , Satler J.D. 2013 . How to fail at species delimitation . Mol. Ecol. 22 : 4369 – 4383 . Google Scholar CrossRef Search ADS PubMed Carstens B.C. , Satler J.D. 2013 . The carnivorous plant described as Sarracenia alata contains two cryptic species . Biol. J. Linnean Soc. 109 : 737 – 746 . Google Scholar CrossRef Search ADS Coyne J.A. , Orr H.A. 1989 . Patterns of speciation in Drosophila . Evolution 43 : 362 – 381 . Google Scholar CrossRef Search ADS PubMed Darwin C. 1859 . The origin of species by means of natural selection: or, the preservation of favored races in the struggle for life . London, UK : John Murray . de León G.P.-P. , Poulin R. 2016 . Taxonomic distribution of cryptic diversity among metazoans: not so homogeneous after all . Biol. Lett. 12 : 20160371 . Google Scholar CrossRef Search ADS PubMed De Queiroz K. 2007 . Species concepts and species delimitation . Syst. Biol. 56 : 879 – 886 . Google Scholar CrossRef Search ADS PubMed Dolman G. 2008 . Evidence for differential assortative female preference in association with refugial isolation of rainbow skinks in Australia”s tropical rainforests . PLoS One 3 : e3499 . Google Scholar CrossRef Search ADS PubMed Dolman G. , Moritz C. 2006 . A multilocus perspective on refugial isolation and divergence in rainforest skinks (Carlia) . Evolution 60 : 573 – 582 . Google Scholar CrossRef Search ADS PubMed Durand E.Y. , Patterson N. , Reich D. , Slatkin M. 2011 . Testing for ancient admixture between closely related populations . Mol. Biol. Evol. 28 : 2239 – 2252 . Google Scholar CrossRef Search ADS PubMed Dynesius M. , Jansson R. 2014 . Persistence of within species lineages: a neglected control of speciation rates . Evolution 68 : 923 – 934 . Google Scholar CrossRef Search ADS PubMed Eaton D.A. , Ree R.H. 2013 . Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae) . Syst. Biol. 62 : 689 – 706 . Google Scholar CrossRef Search ADS PubMed Edwards S.V. , Potter S. , Schmitt C.J. , Bragg J.G. , Moritz C. 2016 . Reticulation, divergence, and the phylogeography–phylogenetics continuum . Proc. Natl. Acad. Sci. U.S.A 113 : 8025 – 8032 . Google Scholar CrossRef Search ADS PubMed Ence D.D. , Carstens B.C. 2011 . SpedeSTEM: a rapid and accurate method for species delimitation . Mol. Ecol. Resour. 11 : 473 – 480 . Google Scholar CrossRef Search ADS PubMed Fitzpatrick B.M. 2002 . Molecular correlates of reproductive isolation . Evolution 56 : 191 – 198 . Google Scholar CrossRef Search ADS PubMed Fouquet A. , Gilles A. , Vences M. , Marty C. , Blanc M. , Gemmell N.J. 2007 . Underestimation of species richness in Neotropical frogs revealed by mtDNA analyses . PLoS One 2 : e1109 . Google Scholar CrossRef Search ADS PubMed Frankham R. , Ballou J.D. , Dudash M.R. , Eldridge M.D. , Fenster C.B. , Lacy R.C. , Mendelson J.R. , Porton I.J. , Ralls K. , Ryder O.A. 2012 . Implications of different species concepts for conserving biodiversity . Biol. Conserv. 153 : 25 – 31 . Google Scholar CrossRef Search ADS Fregin S. , Haase M. , Olsson U. , Alström P. 2012 . Pitfalls in comparisons of genetic distances: a case study of the avian family Acrocephalidae . Mol. Phylogenet. Evol. 62 : 319 – 328 . Google Scholar CrossRef Search ADS PubMed Fujita M.K. , Leaché A.D. , Burbrink F.T. , McGuire J.A. , Moritz C. 2012 . Coalescent-based species delimitation in an integrative taxonomy . Trends Ecol. Evol. 27 : 480 – 488 . Google Scholar CrossRef Search ADS PubMed Funk D.J. , Nosil P. , Etges W.J. 2006 . Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa . Proc. Natl. Acad. Sci. U.S.A 103 : 3209 – 3213 . Google Scholar CrossRef Search ADS PubMed Gómez A. , Serra M. , Carvalho G.R. , Lunt D.H. 2002 . Speciation in ancient cryptic species complexes: evidence from the molecular phylogeny of Brachionus plicatilis (Rotifera) . Evolution 56 : 1431 – 1444 . Google Scholar CrossRef Search ADS PubMed Gomez A. , Wright P.J. , Lunt D.H. , Cancino J.M. , Carvalho G.R. , Hughes R.N. 2007 . Mating trials validate the use of DNA barcoding to reveal cryptic speciation of a marine bryozoan taxon . Proc. R. Soc. Lond. B Biol. Sci. 274 : 199 – 207 . Google Scholar CrossRef Search ADS Grabherr M.G. , Haas B.J. , Yassour M. , Levin J.Z. , Thompson D.A. , Amit I. , Adiconis X. , Fan L. , Raychowdhury R. , Zeng Q. 2011 . Full-length transcriptome assembly from RNA-Seq data without a reference genome . Nat. Biotechnol. 29 : 644 – 652 . Google Scholar CrossRef Search ADS PubMed Graham C.H. , Moritz C. , Williams S.E. 2006 . Habitat history improves prediction of biodiversity in rainforest fauna . Proc. Natl. Acad. Sci. U.S.A 103 : 632 – 636 . Google Scholar CrossRef Search ADS PubMed Greer A. 1997 . A new species of Lampropholis (Squamata: Scincidae) with a restricted, high altitude distribution in eastern Australia . Aust. Zool. 30 : 360 – 368 . Google Scholar CrossRef Search ADS Hebert P.D. , Penton E.H. , Burns J.M. , Janzen D.H. , Hallwachs W. 2004a . Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly Astraptes fulgerator . Proc. Natl. Acad. Sci. U.S.A 101 : 14812 – 14817 . Google Scholar CrossRef Search ADS Hebert P.D. , Stoeckle M.Y. , Zemlak T.S. , Francis C.M. 2004b . Identification of birds through DNA barcodes . PLoS Biol. 2 : e312 . Google Scholar CrossRef Search ADS Hedgecock D. , Ayala F.J. 1974 . Evolutionary divergence in the genus Taricha (Salamandridae) . Copeia 3 : 738 – 747 . Google Scholar CrossRef Search ADS Hedin M. 2015 . High stakes species delimitation in eyeless cave spiders (Cicurina, Dictynidae, Araneae) from central Texas . Mol. Ecol. 24 : 346 – 361 . Google Scholar CrossRef Search ADS PubMed Hewitt G. 2000 . The genetic legacy of the Quaternary ice ages . Nature 405 : 907 – 913 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. 2007 . Description, biology and conservation of a new species of Australian tree frog (Amphibia: Anura: Hylidae: Litoria) and an assessment of the remaining populations of Litoria genimaculata Horst, 1883: systematic and conservation implications of an unusual speciation event . Biol. J. Linnean Soc. 91 : 549 – 563 . Google Scholar CrossRef Search ADS Hoskin C.J. 2014 . A new skink (Scincidae: Carlia) from the rainforest uplands of Cape Melville, north-east Australia . Zootaxa 3869 : 224 – 236 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. , Higgie M. , McDonald K.R. , Moritz C. 2005 . Reinforcement drives rapid allopatric speciation . Nature 437 : 1353 – 1356 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. , Tonione M. , Higgie M. , MacKenzie J.B. , Williams S.E. , VanDerWal J. , Moritz C. 2011 . Persistence in peripheral refugia promotes phenotypic divergence and speciation in a rainforest frog . Am. Nat. 178 : 561 – 578 . Google Scholar CrossRef Search ADS PubMed Hotaling S. , Foley M.E. , Lawrence N.M. , Bocanegra J. , Blanco M.B. , Rasoloarison R. , Kappeler P.M. , Barrett M.A. , Yoder A.D. , Weisrock D.W. 2016 . Species discovery and validation in a cryptic radiation of endangered primates: coalescent based species delimitation in Madagascar”s mouse lemurs . Mol. Ecol. 25 : 2029 – 2045 . Google Scholar CrossRef Search ADS PubMed Ingram G. 1991 . Five new skinks from Queensland rainforests . Mem. Queensl. Mus. 30 : 443 – 453 . Ingram G. , Covacevich J. 1989 . Revision of the genus Carlia (Reptilia, Scincidae) in Australia with comments on Carlia bicarinata of New Guinea . Mem. Queensl. Mus. 27 : 443 – 490 . Isaac N.J. , Mallet J. , Mace G.M. 2004 . Taxonomic inflation: its influence on macroecology and conservation . Trends Ecol. Evol. 19 : 464 – 469 . Google Scholar CrossRef Search ADS PubMed Janzen D.H. , Hajibabaei M. , Burns J.M. , Hallwachs W. , Remigio E. , Hebert P.D. 2005 . Wedding biodiversity inventory of a large and complex Lepidoptera fauna with DNA barcoding . Philos. Trans. R. Soc. Lond. B Biol. Sci. 360 : 1835 – 1845 . Google Scholar CrossRef Search ADS PubMed Jolliffe I.T. 2002 . Principal component analysis . New York (NY) : Springer-Verlag . Jones G. 2017 . Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent . J. Math. Biol. 74 : 447 – 467 . Google Scholar CrossRef Search ADS PubMed Katoh K. , Misawa K. , Kuma K.I. , Miyata T. 2002 . MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform . Nucleic Acids Res. 30 : 3059 – 3066 . Google Scholar CrossRef Search ADS PubMed Kent W.J. 2002 . BLAT—the BLAST-like alignment tool . Genome Res. 12 : 656 – 664 . Google Scholar CrossRef Search ADS PubMed King R.A. , Tibble A.L. , Symondson W.O. 2008 . Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms . Mol. Ecol. 17 : 4684 – 4698 . Google Scholar CrossRef Search ADS PubMed Knowlton N. 1993 . Sibling species in the sea . Ann. Rev. Ecol. Evol. Syst. 24 : 189 – 216 . Google Scholar CrossRef Search ADS Ladner J.T. , Palumbi S.R. 2012 . Extensive sympatry, cryptic diversity and introgression throughout the geographic distribution of two coral species complexes . Mol. Ecol. 21 : 2224 – 2238 . Google Scholar CrossRef Search ADS PubMed Lanfear R. , Frandsen P.B. , Wright A.M. , Senfeld T. , Calcott B. 2016 . PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses . Mol. Biol. Evol. 34 : 772 – 773 . Leaché A.D. , Fujita M.K. 2010 . Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus) . Proc. R. Soc. Lond. B Biol. Sci. 277 1697 : rspb20100662 . Leys M. , Keller I. , Robinson C.T. , Räsänen K. 2017 . Cryptic lineages of a common alpine mayfly show strong life history divergence . Mol. Ecol. 26 : 1670 – 1686 . Google Scholar CrossRef Search ADS PubMed Li H. , Durbin R. 2009 . Fast and accurate short read alignment with Burrows–Wheeler transform . Bioinformatics 25 : 1754 – 1760 . Google Scholar CrossRef Search ADS PubMed Losos J.B. 2011 . Lizards in an evolutionary tree: ecology and adaptive radiation of anoles . Berkeley (CA) : University of California Press . Mallet J. 1995 . A species definition for the modern synthesis . Trends Ecol. Evol. 10 : 294 – 299 . Google Scholar CrossRef Search ADS PubMed Mani G. , Clarke B. 1990 . Mutational order: a major stochastic process in evolution . Proc. R. Soc. Lond. B Biol. Sci. 240 : 29 – 37 . Google Scholar CrossRef Search ADS PubMed Mayr E. 1942 . Systematics and the origin of species, from the viewpoint of a zoologist . Cambridge (MA) : Harvard University Press . McDaniel S.F. , Shaw A.J. 2003 . Phylogeographic structure and cryptic speciation in the trans-Antarctic moss Pyrrhobryum mnioides . Evolution 57 : 205 – 215 . Google Scholar CrossRef Search ADS PubMed McKenna A. , Hanna M. , Banks E. , Sivachenko A. , Cibulskis K. , Kernytsky A. , Garimella K. , Altshuler D. , Gabriel S. , Daly M. 2010 . The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res. 20 : 1297 – 1303 . Google Scholar CrossRef Search ADS PubMed Melville J. , Haines M.L. , Boysen K. , Hodkinson L. , Kilian A. , Date K.L.S. , Potvin D.A. , Parris K.M. 2017 . Identifying hybridization and admixture using SNPs: application of the DArTseq platform in phylogeographic research on vertebrates . R. Soc. Open Sci. 4 : 161061 . Google Scholar CrossRef Search ADS PubMed Meyer M. , Kircher M. 2010 . Illumina sequencing library preparation for highly multiplexed target capture and sequencing . Cold Spring Harbor Protocols 2010 : pdb. prot5448 . Moritz C. , Cicero C. 2004 . DNA barcoding: promise and pitfalls . PLoS Biol. 2 : e354 . Google Scholar CrossRef Search ADS PubMed Moritz C. , Hoskin C. , MacKenzie J.B. , Phillips B. , Tonione M. , Silva N. , VanDerWal J. , Williams S.E. , Graham C. 2009 . Identification and dynamics of a cryptic suture zone in tropical rainforest . Proc. R. Soc. Lond. B Biol. Sci . rspb. 2008.1622 . Nater A. , Mattle-Greminger M.P. , Nurcahyo A. , Nowak M.G. , de Manuel M. , Desai T. , Groves C. , Pybus M. , Sonay T.B. , Roos C. 2017 . Morphometric, behavioral, and genomic evidence for a new Orangutan species . Curr. Biol. 27 : 3487 – 3498.e3410 . Google Scholar CrossRef Search ADS PubMed Nei M. , Li W.-H. 1979 . Mathematical model for studying genetic variation in terms of restriction endonucleases . Proc. Natl. Acad. Sci. U.S.A 76 : 5269 – 5273 . Google Scholar CrossRef Search ADS PubMed Niemiller M.L. , Near T.J. , Fitzpatrick B.M. 2012 . Delimiting species using multilocus data: diagnosing cryptic diversity in the southern cavefish, Typhlichthys subterraneus (Teleostei: Amblyopsidae) . Evolution 66 : 846 – 866 . Google Scholar CrossRef Search ADS PubMed Nosil P. 2008 . Speciation with gene flow could be common . Mol. Ecol. 17 : 2103 – 2106 . Google Scholar CrossRef Search ADS PubMed Nosil P. , Flaxman S.M. 2011 . Conditions for mutation-order speciation . Proc. R. Soc. Lond. B Biol. Sci. 278 : 399 – 407 . Google Scholar CrossRef Search ADS Ogilvie H.A. , Bouckaert R.R. , Drummond A.J. 2017 . StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates . Mol. Biol. Evol. 34 : 2101 – 2114 . Google Scholar CrossRef Search ADS PubMed Olave M. , Solà E. , Knowles L.L. 2014 . Upstream analyses create problems with DNA-based species delimitation . Syst. Biol. 63 : 263 – 271 . Google Scholar CrossRef Search ADS PubMed Oliver P. , Keogh J.S. , Moritz C. 2015 . New approaches to cataloguing and understanding evolutionary diversity: a perspective from Australian herpetology . Aust. J. Zool. 62 : 417 – 430 . Oliver P.M. , Adams M. , Doughty P. 2010 . Molecular evidence for ten species and Oligo-Miocene vicariance within a nominal Australian gecko species (Crenadactylus ocellatus, Diplodactylidae) . BMC Evol. Biol. 10 : 386 . Google Scholar CrossRef Search ADS PubMed Olsson U. , Alström P. , Ericson P.G. , Sundberg P. 2005 . Non-monophyletic taxa and cryptic species—evidence from a molecular phylogeny of leaf-warblers (Phylloscopus, Aves) . Mol. Phylogenet. Evol. 36 : 261 – 276 . Google Scholar CrossRef Search ADS PubMed Padial J.M. , Miralles A. , De la Riva I. , Vences M. 2010 . The integrative future of taxonomy . Front. Zool. 7 : 16 . Google Scholar CrossRef Search ADS PubMed Pante E. , Puillandre N. , Viricel A. , Arnaud Haond S. , Aurelle D. , Castelin M. , Chenuil A. , Destombe C. , Forcioli D. , Valero M. 2015 . Species are hypotheses: avoid connectivity assessments based on pillars of sand . Mol. Ecol. 24 : 525 – 544 . Google Scholar CrossRef Search ADS PubMed Pante E. , Schoelinck C. , Puillandre N. 2014 . From integrative taxonomy to species description: one step beyond . Syst. Biol. 64 : 152 – 160 . Google Scholar CrossRef Search ADS PubMed Pease J.B. , Hahn M.W. 2015 . Detection and polarization of introgression in a five-taxon phylogeny . Syst. Biol. 64 : 651 – 662 . Google Scholar CrossRef Search ADS PubMed Perkins S.L. 2000 . Species concepts and malaria parasites: detecting a cryptic species of Plasmodium . Proc. R. Soc. Lond. B Biol. Sci. 267 : 2345 – 2350 . Google Scholar CrossRef Search ADS Phillips B.L. , Baird S.J. , Moritz C. 2004 . When vicars meet: a narrow contact zone between morphologically cryptic phylogeographic lineages of the rainforest skink, Carlia rubrigularis . Evolution 58 : 1536 – 1548 . Google Scholar CrossRef Search ADS PubMed Pinho C. , Hey J. 2010 . Divergence with gene flow: models and data . Ann. Rev. Ecol. Evol. Syst. 41 : 215 – 230 . Google Scholar CrossRef Search ADS Potter S. , Bragg J.G. , Peter B.M. , Bi K. , Moritz C. 2016 . Phylogenomics at the tips: inferring lineages and their demographic history in a tropical lizard, Carlia amax . Mol. Ecol. 25 : 1367 – 1380 . Google Scholar CrossRef Search ADS PubMed Rabosky D.L. , Matute D.R. 2013 . Macroevolutionary speciation rates are decoupled from the evolution of intrinsic reproductive isolation in Drosophila and birds . Proc. Natl. Acad. Sci. U.S.A 110 : 15354 – 15359 . Google Scholar CrossRef Search ADS PubMed Rannala B. 2015 . The art and science of species delimitation . Curr. Zool. 61 : 846 – 853 . Google Scholar CrossRef Search ADS Richardson B.J. , Baverstock P.R. , Adams M. 1986 . Allozyme electrophoresis: a handbook for animal systematics and population studies . San Diego (CA) : Academic Press . Rissler L.J. , Apodaca J.J. 2007 . Adding more ecology into species delimitation: ecological niche models and phylogeography help define cryptic species in the black salamander (Aneides flavipunctatus) . Syst. Biol. 56 : 924 – 942 . Google Scholar CrossRef Search ADS PubMed Ronquist F. , Teslenko M. , Van Der Mark P. , Ayres D.L. , Darling A. , Höhna S. , Larget B. , Liu L. , Suchard M.A. , Huelsenbeck J.P. 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Syst. Biol. 61 : 539 – 542 . Google Scholar CrossRef Search ADS PubMed Rosenblum E.B. , Sarver B.A. , Brown J.W. , Des Roches S. , Hardwick K.M. , Hether T.D. , Eastman J.M. , Pennell M.W. , Harmon L.J. 2012 . Goldilocks meets Santa Rosalia: an ephemeral speciation model explains patterns of diversification across time scales . J. Evol. Biol. 39 : 255 – 261 . Google Scholar CrossRef Search ADS Rosindell J. , Cornell S.J. , Hubbell S.P. , Etienne R.S. 2010 . Protracted speciation revitalizes the neutral theory of biodiversity . Ecol. Lett. 13 : 716 – 727 . Google Scholar CrossRef Search ADS PubMed Roux C. , Fraisse C. , Romiguier J. , Anciaux Y. , Galtier N. , Bierne N. 2016 . Shedding light on the grey zone of speciation along a continuum of genomic divergence . PLoS Biol. 14 : e2000234 . Google Scholar CrossRef Search ADS PubMed Sasa M.M. , Chippindale P.T. , Johnson N.A. 1998 . Patterns of postzygotic isolation in frogs . Evolution 52 : 1811 – 1820 . Google Scholar CrossRef Search ADS PubMed Schluter D. 2001 . Ecology and the origin of species . Trends Ecol. Evol. 16 : 372 – 380 . Google Scholar CrossRef Search ADS PubMed Schneider C. , Moritz C. 1999 . Rainforest refugia and Australia”s Wet Tropics . Proc. R. Soc. Lond. B Biol. Sci. 266 : 191 – 196 . Google Scholar CrossRef Search ADS Schonrogge K. , Barr B. , Wardlaw J.C. , Napper E. , Gardner M.G. , Breen J. , Elmes G.W. , Thomas J.A. 2002 . When rare species become endangered: cryptic speciation in myrmecophilous hoverflies . Biol. J. Linnean Soc. 75 : 291 – 300 . Google Scholar CrossRef Search ADS Seehausen O. , Takimoto G. , Roy D. , Jokela J. 2008 . Speciation reversal and biodiversity dynamics with hybridization in changing environments . Mol. Ecol. 17 : 30 – 44 . Google Scholar CrossRef Search ADS PubMed Silva A.C.A. , Santos N. , Ogilvie H.A. , Moritz C. 2017 . Validation and description of two new north-western Australian Rainbow skinks with multispecies coalescent methods and morphology . PeerJ 5 : e3724 . Google Scholar CrossRef Search ADS PubMed Singhal S. 2013 . De novo transcriptomic analyses for non model organisms: an evaluation of methods across a multi species data set . Mol. Ecol. Resour. 13 : 403 – 416 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Bi K. 2017 . History cleans up messes: the impact of time in driving divergence and introgression in a tropical suture zone . Evolution 71 : 1888 – 1899 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Moritz C. 2012 . Strong selection against hybrids maintains a narrow contact zone between morphologically cryptic lineages in a rainforest lizard . Evolution 66 : 1474 – 1489 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Moritz C. 2013 . Reproductive isolation between phylogeographic lineages scales with divergence . Proc. R. Soc. Lond. B Biol. Sci. 280 : 20132246 . Google Scholar CrossRef Search ADS Sites J.W. , Barton N.H. , Reed K.M. 1995 . The genetic structure of a hybrid zone between two chromosome races of the Sceloporus grammicus complex (Sauria, Phrynosomatidae) in Central Mexico . Evolution 49 : 9 – 36 . Google Scholar CrossRef Search ADS PubMed Skinner A. , Hugall A.F. , Hutchinson M.N. 2011 . Lygosomine phylogeny and the origins of Australian scincid lizards . J. Biogeogr. 38 : 1044 – 1058 . Google Scholar CrossRef Search ADS Slater G.S.C. , Birney E. 2005 . Automated generation of heuristics for biological sequence comparison . BMC Bioinformatics 6 : 31 . Google Scholar CrossRef Search ADS PubMed Smith M.A. , Rodriguez J.J. , Whitfield J.B. , Deans A.R. , Janzen D.H. , Hallwachs W. , Hebert P.D. 2008 . Extreme diversity of tropical parasitoid wasps exposed by iterative integration of natural history, DNA barcoding, morphology, and collections . Proc. Natl. Acad. Sci. U.S.A 105 : 12359 – 12364 . Google Scholar CrossRef Search ADS PubMed Smith M.A. , Woodley N.E. , Janzen D.H. , Hallwachs W. , Hebert P.D. 2006 . DNA barcodes reveal cryptic host-specificity within the presumed polyphagous members of a genus of parasitoid flies (Diptera: Tachinidae) . Proc. Natl. Acad. Sci. U.S.A 103 : 3657 – 3662 . Google Scholar CrossRef Search ADS PubMed Stewart K.A. , Lougheed S.C. 2013 . Testing for intraspecific postzygotic isolation between cryptic lineages of Pseudacris crucifer . Ecol. Evol. 3 : 4621 – 4630 . Google Scholar CrossRef Search ADS PubMed Struck T.H. , Feder J.L. , Bendiksby M. , Birkeland S. , Cerca J. , Gusarov V.I. , Kistenich S. , Larsson K.-H. , Liow L.H. , Nowak M.D. 2018 . Finding evolutionary processes hidden in cryptic species . Trends Ecol. Evol. 33 : 153 – 163 . Google Scholar CrossRef Search ADS PubMed Stuart B.L. , Inger R.F. , Voris H.K. 2006 . High level of cryptic species diversity revealed by sympatric lineages of Southeast Asian forest frogs . Biol. Lett. 2 : 470 – 474 . Google Scholar CrossRef Search ADS PubMed Suatoni E. , Vicario S. , Rice S. , Snell T. , Caccone A. 2006 . An analysis of species boundaries and biogeographic patterns in a cryptic species complex: the rotifer—Brachionus plicatilis . Mol. Phylogenet. Evol. 41 : 86 – 98 . Google Scholar CrossRef Search ADS PubMed Sukumaran J. , Knowles L.L. 2017 . Multispecies coalescent delimits structure, not species . Proc. Natl. Acad. Sci. U.S.A 114 : 1607 – 1612 . Google Scholar CrossRef Search ADS PubMed Sunnucks P. , Hales D.F. 1996 . Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae) . Mol. Biol. Evol. 13 : 510 – 524 . Google Scholar CrossRef Search ADS PubMed Tellier F. , Vega J.A. , Broitman B.R. , Vasquez J.A. , Valero M. , Faugeron S. 2011 . The importance of having two species instead of one in kelp management: the Lessonia nigrescens species complex . Cah. Biol. Mar. 52 : 455 – 465 . Turner J.R. 1967 . Why does the genotype not congeal? Evolution 21 : 645 – 656 . Google Scholar CrossRef Search ADS PubMed Vilas R. , Criscione C. , Blouin M. 2005 . A comparison between mitochondrial DNA and the ribosomal internal transcribed regions in prospecting for cryptic species of platyhelminth parasites . Parasitology 131 : 839 – 846 . Google Scholar CrossRef Search ADS PubMed Winger B.M. , Bates J.M. 2015 . The tempo of trait divergence in geographic isolation: avian speciation across the Marañon Valley of Peru . Evolution 69 : 772 – 787 . Google Scholar CrossRef Search ADS PubMed Witt J.D. , Threloff D.L. , Hebert P.D. 2006 . DNA barcoding reveals extraordinary cryptic diversity in an amphipod genus: implications for desert spring conservation . Mol. Ecol. 15 : 3073 – 3082 . Google Scholar CrossRef Search ADS PubMed Yang Z. , Rannala B. 2014 . Unguided species delimitation using DNA sequence data from multiple loci . Mol. Biol. Evol. 31 : 3125 – 3135 . Google Scholar CrossRef Search ADS PubMed Zhang C. , Zhang D.-X. , Zhu T. , Yang Z. 2011 . Evaluation of a Bayesian coalescent method of species delimitation . Syst. Biol. 60 : 747 – 761 . Google Scholar CrossRef Search ADS PubMed Zhang J. , Kobert K. , Flouri T. , Stamatakis A. 2013 . PEAR: a fast and accurate Illumina Paired-End reAd mergeR . Bioinformatics 30 : 614 – 620 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press

# A Framework for Resolving Cryptic Species: A Case Study from the Lizards of the Australian Wet Tropics

, Volume Advance Article – Apr 4, 2018
15 pages

/lp/ou_press/a-framework-for-resolving-cryptic-species-a-case-study-from-the-nQVRD7MnFl
Publisher
Oxford University Press
ISSN
1063-5157
eISSN
1076-836X
D.O.I.
10.1093/sysbio/syy026
Publisher site
See Article on Publisher Site

### Abstract

Abstract As we collect range-wide genetic data for morphologically-defined species, we increasingly unearth evidence for cryptic diversity. Delimiting this cryptic diversity is challenging, both because the divergences span a continuum and because the lack of overt morphological differentiation suggests divergence has proceeded heterogeneously. Herein, we address these challenges as we diagnose and describe species in three co-occurring species groups of Australian lizards. By integrating genomic and morphological data with data on hybridization and introgression from contact zones, we explore several approaches—and their relative benefits and weaknesses—for testing the validity of cryptic lineages. More generally, we advocate that genetic delimitations of cryptic diversity must consider whether these lineages are likely to be durable and persistent through evolutionary time. Exome capture, cryptic species, phylogeography, species delimitation, squamates, taxonomy Cryptic species, or taxa that are morphologically similar but genetically divergent, exemplify the two major challenges of species delimitation. First, species form on a continuum (Darwin 1859; Mayr 1942; Mallet 1995; De Queiroz 2007). As populations differentiate across space and time, they gradually become more divergent. As reflected in the debate over defining operational taxonomic units via DNA barcoding (Moritz and Cicero 2004), deciding how much divergence is sufficient to name lineages as species can be arbitrary. As with morphologically distinct lineages, diagnosing cryptic lineages presents this challenge because they fall through the full range of the divergence continuum (Hedgecock and Ayala 1974; Gómez et al. 2002; McDaniel and Shaw 2003). Second, speciation proceeds heterogeneously across many axes. We typically recover correlations across axes—e.g., rates of trait evolution such as song and mitochondrial divergence (Winger and Bates 2015). However, when axes of differentiation are discordant—for example, when phenotypic disparity is high and genetic divergence is low – the status of lineages becomes ambiguous. By definition, cryptic lineages have diverged heterogeneously (Bickford et al. 2007)—they are genetically-distinct groups that exhibit little or no morphological divergence. The taxonomic process of naming a species is a binary exercise—either a lineage is recognized as a species or not—and accounting for heterogeneity in this binary framework remains a challenge. Biodiversity researchers increasingly face these challenges because we are increasingly discovering new cryptic lineages (Bickford et al. 2007). A confluence in genetic advances and broader geographic sampling has led to rapid increases in the number of cryptic species, such that a quarter of articles published in Zoological Record Plus mention cryptic species (Bickford et al. 2007; see also de León and Poulin 2016). Cryptic species comprise a significant proportion of the diversity in some regions (e.g., tropics; Smith et al. 2008) and taxonomic groups (e.g., reptiles; Oliver et al. 2010), and recently, putative cryptic species have been identified in high-profile threatened species like orangutans (Nater et al. 2017). These findings, and their implications for evolutionary biology and conservation (Frankham et al. 2012), emphasize the need for a more rigorous framework to assess the taxonomic status of cryptic lineages (Adams et al. 2014; Struck et al. in press). In this work, we propose a framework that diagnoses those cryptic lineages that are expected to be sufficiently durable to contribute to build-up of diversity over time and space. Because speciation is a continuum, we expect that many nascent species are lost to hybridization and extinction as part of the protracted speciation process (Rosindell et al. 2010; Dynesius and Jansson 2014). As such, we adopt the biological species concept (BSC, (Mayr 1942)), which defines species as units that exhibit barriers to reproduction and are thus more likely to persist through time. Although some might find this definition overly restrictive, we apply it here in hopes of avoiding “taxonomic over-inflation” (Isaac et al. 2004). However, how do we diagnose populations that are likely to have substantial reproductive isolation (RI)? In allopatry, the degree of morphological difference is expected to correlate with the extent of RI (Mayr 1942; Bolnick et al. 2006; Funk et al. 2006). Thus, when genetic and phenotypic divergence concur, species delimitation is typically uncontroversial. For cryptic species, where we cannot use phenotypic divergence as a proxy for RI, we must instead use multiple lines of evidence to assess likelihood of strong RI. A popular approach to assess cryptic lineages is to apply statistical species delimitation to multilocus genetic data (Fujita et al. 2012; Carstens et al. 2013), which some argue makes species delimitation more objective (Rannala 2015). An illuminating line of research has explored the parameter space under which these coalescent-based methods are expected to return statistically robust results (Zhang et al. 2011; Olave et al. 2014). What remains to be seen is if these statistically robust lineages are also biologically robust (Sukumaran and Knowles 2017). That is, will newly delimited lineages remain distinct through changing geographies and environments, or will they be mere evolutionary ephemera lost to hybridization and/or extinction (Seehausen et al. 2008; Rosenblum et al. 2012; Dynesius and Jansson 2014)? A stronger and more direct approach to delimitation is to test for strongly restricted gene flow in sympatry or parapatry (Richardson et al. 1986; Adams et al. 2014). This can be done either by sampling a modest number of individuals in sympatry or by assessing the extent of genetic introgression through intensive analysis of contact zones between parapatric taxa. However, when candidate taxa are allopatric, assessing the likelihood of strong RI is even more challenging, as long recognized under operational versions of the BSC. To the extent that RI is time-dependent (Coyne and Orr 1989; Sasa et al. 1998; Fitzpatrick 2002; Roux et al. 2016), another approach is to extrapolate from closely-related taxa where the relationship between divergence and extent of RI has already been determined. In contrast to a “bar-coding gap” (Hebert et al. 2004b), this approach uses genome-scale evidence to assess the likelihood of strong RI rather than patterns of genetic divergence across nominal species vs. populations. We explore these multiple approaches to species delimitation through the study of morphologically cryptic, phylogeographic lineages in the lizards of the Australian wet tropics (AWT). The AWT is a narrow region of rainforest in northeast Queensland, Australia (Fig. 1). During repeated glacial cycles in the Quaternary (Graham et al. 2006), the rainforest and, accordingly, the species endemic to these rainforests, were split across two major refugia. Populations of these rainforest species diverged across these refugia; comparative data have recovered deep phylogeographic splits within species across more than twenty taxa (Moritz et al. 2009). Morphological analyses show limited phenotypic divergence among phylogeographic lineages (Schneider and Moritz 1999; Hoskin et al. 2005; Hoskin et al. 2011). Subsequent contact zone studies showed that some of these phylogeographic lineages are reproductively isolated (Phillips et al. 2004; Hoskin et al. 2005; Singhal and Moritz 2013). Figure 1. View largeDownload slide Geographic distribution of lineages in the three groups; boundaries were inferred after sequencing an average of 27 individuals for mtDNA and 12 for nDNA (Supplementary Table S1, Fig. S2 available on Dryad). Crosses represent localities for samples included in this study. A) The five lineages in the ‘Carlia rubrigularis’ group: C. wundalthini, C. rubrigularis N, C. rubrigularis S, C. rhomboidalis N, and C. rhomboidalis S, B) the four lineages in the ‘Lampropholis coggeri’ group: L. coggeri N, L. coggeri C, L. coggeri S, and L.coggeri EU, and C) the four lineages in the ‘L. robertsi’ group: L. robertsi CU, L. robertsi BFAU, L.robertsi TU, L. robertsi BK. In the inset map, the distribution of the rainforest is shown in light green. Like many phylogeographic lineages, these lineages are geographically circumscribed, and their ranges either are geographically proximate or narrowly overlapping. Pictures courtesy of B. Phillips, C. Peng, and S. Zozaya (L-R). Figure 1. View largeDownload slide Geographic distribution of lineages in the three groups; boundaries were inferred after sequencing an average of 27 individuals for mtDNA and 12 for nDNA (Supplementary Table S1, Fig. S2 available on Dryad). Crosses represent localities for samples included in this study. A) The five lineages in the ‘Carlia rubrigularis’ group: C. wundalthini, C. rubrigularis N, C. rubrigularis S, C. rhomboidalis N, and C. rhomboidalis S, B) the four lineages in the ‘Lampropholis coggeri’ group: L. coggeri N, L. coggeri C, L. coggeri S, and L.coggeri EU, and C) the four lineages in the ‘L. robertsi’ group: L. robertsi CU, L. robertsi BFAU, L.robertsi TU, L. robertsi BK. In the inset map, the distribution of the rainforest is shown in light green. Like many phylogeographic lineages, these lineages are geographically circumscribed, and their ranges either are geographically proximate or narrowly overlapping. Pictures courtesy of B. Phillips, C. Peng, and S. Zozaya (L-R). Herein, we focus on three species groups—the ‘Carlia rubrigularis’, ‘Lampropholis coggeri’, and ‘Lampropholis robertsi’ groups—which are part of the broader radiation of Eugongylus lizards (family: Scincidae) (Skinner et al. 2011). These groups are ecologically-similar; all are small (30–55 mm) skinks of the leaf-litter in the rainforests of the AWT. Previous multilocus analyses revealed several phylogeographic lineages within each of these nominal taxa (Dolman and Moritz 2006; Bell et al. 2010). Like many phylogeographic units, these lineages are mostly morphologically cryptic and geographically circumscribed, and their ranges are either geographically proximate or narrowly overlapping. With an eye to integrative taxonomy (Padial et al. 2010), we synthesize data on genetics, morphology, and reproductive isolation assessed in contact zones to resolve the species status of lineages within these three groups and to formally revise their taxonomy. More generally, we use these lizards as a data-rich case study to explore the challenges of delimiting species among cryptic lineages that are parapatric or allopatric. Methods Sampling, Data Collection, and Data Processing In this study, we analyze genetic data for individuals across three species groups across five nominal species and 13 putative lineages. The ‘Carlia rubrigularis’ group consists of five lineages: C. rubrigularis, northern Wet Tropics (N); C. rubrigularis, southern Wet Tropics (S); C. rhomboidalis, northern mid-east Queensland (N); C. rhomboidalis, southern mid-east Queensland (S); and C. wundalthini at Cape Melville (Fig. 1). The ‘Lampropholis coggeri’ group consists of four lineages: L. coggeri, northern Wet Tropics (N); L. coggeri, central Wet Tropics (C); L. coggeri, southern Wet Tropics (S); and L. coggeri in the Mt Elliot uplands (EU). The montane ‘Lampropholis robertsi’ group consists of four allopatric lineages: L. robertsi, Carbine Tableland uplands (CU); L. robertsi, Thornton Peak uplands (TU); L. robertsi, Mt Bellenden Ker (BK); and L. robertsi, Mt Bartle Frere and southern Atherton Tablelands (BFAU). These lineages had been previously described through the sequencing of an average of 27 geographically dispersed individuals per lineage for mtDNA and 12 for multilocus nDNA (Dolman and Moritz 2006; Bell et al. 2010). Because these lineages were sampled extensively in previous genetic analyses, and because they are circumscribed geographically (Fig. 1), we limited our sampling to a few individuals per lineage (Fig. 1; Supplementary Tables S1 and S2 available on Dryad at http://dx.doi.org/10.5061/dryad.g7v1b). Further, across these lineages, three lineage-pairs meet in hybrid zones, two of which are very narrow ($$<$$1 km) and the other of which is wider ($$<$$10 km) (Singhal and Moritz 2013; Singhal and Bi 2017). Given that introgression is geographically limited relative to lineage ranges, we selected individuals well removed from contemporary contact zones to estimate species trees, test for lineage-wide introgression, and to apply statistical delimitation methods. We further included two outgroups: C. storri for the ‘C. rubrigularis’ group and L. amicula for the ‘L. coggeri’ and ‘L. robertsi’ groups. For L. coggeri N and C, we supplemented our sampling by including previously-published transcriptome data (Singhal 2013), which we analyzed using the same approach outlined below. In total, we sampled 25 individuals (Supplementary Tables S1 and S2 available on Dryad). We sequenced a homologous set of 3320 loci across all individuals using an exome capture approach previously described (Bragg et al. 2016). Briefly, we identified homologous exons across the Eugongylus skink clade from transcriptome data of three species. After filtering the exons for GC-content and length, we included probes specific to each of the three species in the capture array. For each sample, we extracted DNA (Sunnucks and Hales 1996) and generated uniquely-barcoded libraries (Meyer and Kircher 2010). Individuals were then pooled in equimolar amounts along with other Eugongylus group taxa (Bragg et al. 2018), and we captured our target exons using a SeqCap EZ Developer Probes following manufacturer”s instructions. The captured libraries were subsequently sequenced along with samples for other projects on the Illumina HiSeq2000 and 2500 for 100 bp paired-end reads. Our bioinformatics pipeline is designed for both population genetic and phylogenetic analyses and thus generates both variant data for each lineage and locus alignments across all haplotypes in a given group. The basic approach follows, with further details available in Singhal et al. (2017). Following de-multiplexing, we trimmed low quality sequence and adaptor sequence from reads using Trimmomatic v0.36 and merged overlapping reads using pear v0.9.10 (Zhang et al. 2013; Bolger et al. 2014). We generated an assembly for each individual with Trinity v2.3.2 and identified assembled contigs homologous to our original targets with blat v36x1 (Kent 2002; Grabherr et al. 2011). For each exon, we picked the best matching contig across all individuals in the lineage to generate a pseudo-reference genome. We then annotated the coding sequence for these exons using exonerate v2.2.0 (Slater and Birney 2005). To generate variant data, we mapped reads from each individual to the pseudo-reference genome using bwa v0.7.12 (Li and Durbin 2009), called variant and invariant sites using gatk v3.6 Unifiedgenotyper, filtered out sites with quality $$<$$20, quality depth (QD) $$<$$5, and coverage $${<}20{\times}$$, and then determined haplotypes using gatkReadBackedPhasing (McKenna et al. 2010). Eight percent of sites were unphased; for these sites, we randomly phased them with respect to other phased blocks. We then generated multi-lineage alignments across all haplotypes with mafft v7.294 (Katoh et al. 2002). Like many phylogeographic studies, we first identified these lineages by sequencing mitochondrial loci. For these taxa, mtDNA and nDNA are highly congruent except within narrow contact zones, and mtDNA provides our most complete understanding of geographic limits (Dolman and Moritz 2006; Bell et al. 2010). Accordingly, we downloaded all geo-referenced NADH dehydrogenase subunit 4 (ND4) data for these groups from GenBank (Supplementary Table S3 available on Dryad). We aligned the data using mafft and identified the coding sequence using exonerate. Phylogenetic Analyses We reconstructed the evolutionary history of 13 lineages and two outgroups by using starbeast2 v0.13.5, a coalescent multilocus method (Ogilvie et al. 2017). Individuals were assigned to lineages following their ‘putative lineage’ designations (Supplementary Table S2 available on Dryad). We filtered loci to only include those that were $$\geqslant$$75% complete across samples and to remove loci $$\geqslant$$1500 bp because longer loci are more likely to capture recombination events. Because of the computational demands of running starbeast2, we generated three random subsamples of 200 loci each from the remaining loci. We ran starbeast2 on these random subsets for 500e6 generations sampling every 1e5 generations. Each locus was assigned to its own strict clock and a GTR model of molecular evolution. Because we lack robust age constraints for nodes in this species tree, we instead inferred branch lengths in units of substitutions per site. For the mitochondrial phylogenetic analysis, we determined the best-fitting partitioning strategy using PartitionFinder2 (Lanfear et al. 2016). We then inferred the mtDNA gene trees using MrBayes v3.2.6, running two runs of four chains each for 50e6 steps (Ronquist et al. 2012). We set the branch length prior to exponential (100); the default prior overestimates branch lengths when the majority of bifurcations occur within-species (Brown et al. 2009). Population Genetic Analyses Our population genetic analyses were aimed at describing basic patterns of diversity, divergence, and current and historical introgression among the lineages in these groups. For each lineage, we calculated within-lineage genetic diversity ($$\pi$$; Nei and Li 1979) for both the exome and mtDNA data, across silent sites only. For each pairwise-comparison between lineages in a species group, we calculated raw and net divergence ($$d_{xy}$$ and $$d_{a})$$ for the exome and mtDNA data (Nei and Li 1979), again across silent sites only. The patterns of divergence and diversity among individuals within a lineage and across lineages confirm our lineage assignments based on mtDNA and previous multilocus nuclear data (Dolman and Moritz 2006; Bell et al. 2010). To test for historical introgression, we used the D-statistic (Durand et al. 2011). For the topology [(P1, P2), P3), outgroup], the D-statistic distinguishes if P1 and P3 exhibit incomplete lineage sorting or introgression by comparing site patterns across these four tips. We conducted this test across all possible comparisons within species groups, except for sister taxa, which cannot be accommodated in the D-statistic framework. Based on the overall species tree topology (Fig. 2), we determined the appropriate species to use as P2 (Supplementary Table S5 available on Dryad). For all Lampropholis comparisons, we used L. amicula as the outgroup, and for Carlia, we used C. storri. Because of our limited within-lineage sampling, we implemented the version of the test designed for fixed variants. Thus, we removed any site that was missing or polymorphic for any lineage in a species group. For the remaining sites, we calculated the D-statistic across all possible species comparisons. To assess significance, we calculated the standard deviation across 200 bootstraps and used a one-tailed $$z$$-test (Eaton and Ree 2013). For L. robertsi, we further tested for introgression using the D$$_{\mathrm{FOIL}}$$ approach designed for five-taxon symmetric topologies (Pease and Hahn 2015); we could not apply this method to other groups because their topologies are asymmetric. This method confirmed our D-statistic results, so we do not discuss them further. Figure 2. View largeDownload slide Species tree for the included taxa as inferred using STARBEAST2 with 200 randomly selected exons. This topology and branching times are robust across multiple random samples (Supplementary Fig. S2 available on Dryad). Nodes with $$\geqslant$$95% local posterior probability are indicated with filled circles. Nodes for which BPP inferred $$\geqslant$$95% probability of a speciation event are shown by open circles. STACEY results match BPP results and are not shown. Arrows indicate pairwise relationships for which there is significant evidence for historical introgression (Supplementary Table S6 available on Dryad). The phylogeographic lineages included in this study all diverged over a relatively narrow span of time. Figure 2. View largeDownload slide Species tree for the included taxa as inferred using STARBEAST2 with 200 randomly selected exons. This topology and branching times are robust across multiple random samples (Supplementary Fig. S2 available on Dryad). Nodes with $$\geqslant$$95% local posterior probability are indicated with filled circles. Nodes for which BPP inferred $$\geqslant$$95% probability of a speciation event are shown by open circles. STACEY results match BPP results and are not shown. Arrows indicate pairwise relationships for which there is significant evidence for historical introgression (Supplementary Table S6 available on Dryad). The phylogeographic lineages included in this study all diverged over a relatively narrow span of time. Finally, we collated previously-published data on reproductive isolation at three contact zones: L. coggeri N and C, L. coggeri C and S, and C. rubrigularis N and S (Phillips et al. 2004; Dolman 2008; Singhal and Moritz 2012, 2013; Singhal and Bi 2017). These studies sampled densely through each contact zone to infer current rates of hybridization and to determine patterns of introgression across the genome and geography. Statistical Species Delimitation One of the most common ways to validate cryptic lineages is through multilocus coalescent-based (MSC) approaches (Fujita et al. 2012). Accordingly, and as recommended by Rannala (2015), we used two MSC approaches (bpp v3.3a and stacey v1.2.4) to test species boundaries across these groups (Yang and Rannala 2014; Jones 2017). Using the same filtering as in our starbeast2 analyses, we selected three random samples of 100 loci per species group and generated input files for bpp and stacey. We used the species tree inferred from the starbeast2 analyses (see Phylogenetic Analyses) as the guide tree for BPP. For the ‘L. coggeri’ group, we used two topologies of the species tree that reflected uncertainty in the placement of L. coggeri EU (Supplementary Fig. S1 available on Dryad). We then ran bpp for 500,000 generations across three sets of priors to ensure our results were robust to prior specification. These priors were: 1) $$\rm{\theta} \sim$$ (2, 2000), $$\tau \sim$$ (2, 2000), 2) $$\theta \sim$$ (1, 10), $$\tau \sim$$ (1, 10), and 3) $$\theta \sim$$ (1, 10), $$\tau \sim$$(2, 2000). We ensured that we had 20–80% acceptance rate; having too high or too low of acceptance rates can affect results. We ran stacey for 1e7 generations. Each locus was set to have its own clock rate and own substitution model under a HKY model. Priors were: collapse height $$=$$ 0.001, growth rate $$\sim$$ lognormal(mean $$=$$ 5, sd $$=$$ 2); collapse weight $$\sim$$ uniform (0,1); population prior scale $$\sim$$ lognormal (mean $$= -7$$, sd $$=$$ 2), and relative death rate $$\sim \beta (\alpha = 1, \beta = 8$$). Species delimitations were determined using speciesDA with a burn-in of 10% and a collapse height of 0.0001 (Jones 2017). Analyses showed that results were robust across collapse heights from 0.0001 to 0.001. Morphological Data and Analysis We assessed morphological differences in scale and body measurement traits between the major genetic lineages in each of the three species groups in the Wet Tropics region. For the ‘C. rubrigularis’ species group, we excluded C. wundalthini and C. rhomboidalis because they differ in breeding colors and to some degree body size and shape (Dolman 2008; Hoskin 2014). Further, we combined individuals from L. coggeri N and C, from L. robertsi CU and TU, and L. robertsi BK and BFAU because our analyses of the genomic evidence suggested each of these lineage pairs are not distinct and should be collapsed (see the Discussion). A single investigator took scale counts and body measurements on an average of 26 specimens per lineage (Supplementary Appendix 1, Table S4 available on Dryad). The traits measured summarize morphological characteristics that affect how lizards function in their environment (e.g., relative limb length affects locomotion; Losos 2011) and are standardly used to delimit skink species (Ingram and Covacevich 1989; Ingram 1991; Greer 1997). The scale traits counted were the number of supraciliaries, infralabials, supralabials, midbody scale rows, paravertebrals, and lamellae under the fourth toe (Hoskin 2014). We found no to little variation for these scale traits within species groups and thus performed no further analyses (Supplementary Table A4 available on Dryad). We measured five aspects of body size and shape: total head and body length (snout to vent length, SVL), distance between the front and hindlimbs (axilla to groin length, AG), length of the hindlimb (L2), head length (HL), and head width (HW) (Hoskin 2014). Only adults were measured, defined as individuals with SVL greater than 48 mm, 44 mm, and 32 mm in C. rubrigularis, L. robertsi, and L. coggeri, respectively. For each species group, we used multivariate analyses to determine if lineages differed significantly in size or shape. We first used a principal component analysis (PCA) to summarize across all five body measurements (Jolliffe 2002). We then tested if body size (PC1) and body shape (PCs 2–5) varied significantly across lineages. All analyses were nested by region within lineage to account for among-region variation within lineages. “Region” represents different mountain ranges, tablelands, and lowland areas, and roughly matched the subregions defined for the Wet Tropics (Bell et al. 2010). For body size, we used a nested univariate analysis of variance (ANOVAs) of PC1 between lineages within each species group. For body shape, we used a nested multivariate analysis of variance (MANOVA) of PCs 2–5. Significance was assessed using Roy”s Greatest Root. For those species groups where we included more than two lineages, we performed nested univariate (size) or multivariate (shape) planned contrasts to test which lineages differed significantly. All analyses were conducted in sas v9.2. Results Analysis of Genetic Data Our exome capture approach worked well across all lineages. On average, we recovered an average of 2.29 Mb per individual across 2668 loci at an average coverage of 112$$\times$$ (Supplementary Table S2 available on Dryad). The inferred topology is well-supported and is consistent with previous phylogenetic hypotheses based on mtDNA data (Fig. 2, Supplementary Fig. S2 available on Dryad; Bell et al. 2010) and other phylogenomic analyses (Bragg et al. 2018). Branching times and topologies were quantitatively and qualitatively similar across replicate analyses of starbeast2 (Supplementary Fig. S1 available on Dryad). The branching times between these lineages all occur within a narrow range and highlight that two lineages of C. rubrigularis as currently recognized are polyphyletic. Further, the splitting patterns generally agree with the biogeographic relationships between species. For example, L. robertsi BFAU is sister to L. robertsi BK, and the two lineages occur as adjacent montane isolates (Figs. 1 and 2). The exception to this congruence is a leapfrog distribution of L. coggeri EU, which is sister to L. coggeri N and C rather than the geographically adjacent L. coggeri S (Figs. 1 and 2). We inferred a 4$$\times$$ to 9$$\times$$ range of genetic divergences between lineages within nominal species, with nuclear $$d_{xy}$$ at silent sites ranging from 0.5% to 1.95% and nuclear $$d_{a}$$ ranging from 0.18% to 1.68% (Fig. 3, Supplementary Table S5 available on Dryad). Lineages associated with small ranges such as C. wundalthini and the montane lineages of L. robertsi showed lower levels of within-population diversity. Measures of mtDNA and nuclear $$d_{xy}$$ and $$d_{a}$$ were correlated ($$d_{xy}$$; $$r = 0.61$$, $$P$$-value $$=$$ 0.016; $$d_{a}$$; $$r = 0.57$$, $$P$$-value $$=$$ 0.025; Fig. 3A, B). As for inferred branching times, divergences between some cryptic lineages were significantly greater than those between nominal species (Fig. 3C). Figure 3. View largeDownload slide Patterns of divergence between pairwise lineage-comparisons for (A) $$d_{xy}$$ for silent sites in mitochondrial DNA (mtDNA) and nuclear DNA (nDNA), (B) $$d_{a}$$ at silent sites in mtDNA and nDNA, and (C) branching times in units of substitutions per site per million years and $$d_{a}$$ at mtDNA; branching times were inferred from the tree depicted in Fig. 1. Each pairwise comparison is coded as either being between 1) recognized: two lineages that were already recognized at the species-level, 2) elevated: two lineages that were elevated in the current study, or 3) population: lineages for which there is insufficient evidence to elevate them to species. Arrows identify the three pairwise comparisons for which we have data on reproductive isolation from contact zone studies (Fig. 4), the gray line indicates the transition point at which we first recover evidence for isolation between lineages (i.e., between Carlia rubrigularis N and S). Many of the lineages that we propose to elevate are more diverged than recognised species. Figure 3. View largeDownload slide Patterns of divergence between pairwise lineage-comparisons for (A) $$d_{xy}$$ for silent sites in mitochondrial DNA (mtDNA) and nuclear DNA (nDNA), (B) $$d_{a}$$ at silent sites in mtDNA and nDNA, and (C) branching times in units of substitutions per site per million years and $$d_{a}$$ at mtDNA; branching times were inferred from the tree depicted in Fig. 1. Each pairwise comparison is coded as either being between 1) recognized: two lineages that were already recognized at the species-level, 2) elevated: two lineages that were elevated in the current study, or 3) population: lineages for which there is insufficient evidence to elevate them to species. Arrows identify the three pairwise comparisons for which we have data on reproductive isolation from contact zone studies (Fig. 4), the gray line indicates the transition point at which we first recover evidence for isolation between lineages (i.e., between Carlia rubrigularis N and S). Many of the lineages that we propose to elevate are more diverged than recognised species. Figure 4. View largeDownload slide Evidence for rates of hybridization and introgression at three contact zones between lineages included in this analysis: Carlia rubrigularis N and S, Lampropholis coggeri C and S, and L. coggeri N and C. Contacts are listed in the legend in order of least to most divergent. A) Percent of individuals in the center of contact zones that were identified as hybrid. A hybrid individual was defined as individuals that had $$\geqslant$$10% membership in both parental species as determined by STRUCTURE. B) Distribution of cline widths in the contact zone across an average of 9.5K clines. C) The extent of linkage disequilibrium in each contact zone. Moran”s I measures the autocorrelation in cline widths across the genome, which serves as a proxy for linkage disequilibrium. These genetic estimates of reproductive isolation show evidence for selection against hybrids in C. rubrigularis N and S and L. coggeri C and S. Figure 4. View largeDownload slide Evidence for rates of hybridization and introgression at three contact zones between lineages included in this analysis: Carlia rubrigularis N and S, Lampropholis coggeri C and S, and L. coggeri N and C. Contacts are listed in the legend in order of least to most divergent. A) Percent of individuals in the center of contact zones that were identified as hybrid. A hybrid individual was defined as individuals that had $$\geqslant$$10% membership in both parental species as determined by STRUCTURE. B) Distribution of cline widths in the contact zone across an average of 9.5K clines. C) The extent of linkage disequilibrium in each contact zone. Moran”s I measures the autocorrelation in cline widths across the genome, which serves as a proxy for linkage disequilibrium. These genetic estimates of reproductive isolation show evidence for selection against hybrids in C. rubrigularis N and S and L. coggeri C and S. Our D-statistic tests for historical introgression among non-sister lineages recovered four likely cases of historical introgression between lineage-pairs: C. rubrigularis N and S; C. rubrigularis N and C. rhomboidalis; L. coggeri S and EU; and L. coggeri C and S (Fig. 2, Supplementary Table S6 available on Dryad). No signature of historical introgression was recovered for strongly allopatric populations—e.g., among montane isolates of the ‘L. robertsi’ species group, and C. wundalthini vs. C. rubrigularis. Our previous results from analyses of hybridization and introgression at contact zones show that C. rubrigularis N and S and L. coggeri C and S exhibit 1) a moderate proportion of hybrids in the center of the contact zones, 2) narrow cline widths across the genome, and 3) auto-correlation in cline widths across physical distances, all indicative of extensive disequilibrium in hybrids and substantial RI (Fig. 4). The less divergent lineage-pair, L. coggeri N and C, shows none of the same patterns, with evidence of extensive introgression across the genome and geography. Statistical species delimitation supported all lineages as species. bpp returned $$\geqslant$$95% probability for a speciation event at all nodes in our guide trees, and stacey inferred each species as a unique cluster (Fig. 2). This result was robust to priors and, for the “L. coggeri” group, uncertainty in the topology. Analysis of Morphological Data For all three species groups, PCAs resulted in a PC1 that accounted for most of the variation (67–84%). It was loaded highly and positively by all body measurements and indicated body size (Supplementary Tables S7–S9 available on Dryad). The remaining four PCs in each species group represented variation in body shape (Supplementary Tables S7–S9 available on Dryad). PC2 accounted for 10–16% of the variation in the species groups and was loaded most heavily by relative AG length (positive) in ‘C. rubrigularis’ and both AG length (positive) and relative L2 length (negative) in ‘L. robertsi’ and ‘L. coggeri’ (Supplementary Tables S7–S9 available on Dryad). The three major lineages in the ‘L. coggeri’ group show no differences in body size (‘lineage’ effect on PC1, $$F_{\mathrm{2,58}} = 2.46$$, $$P = 0.09$$), but they do differ in body shape (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{8,10}} = 3.46$$, $$P = 0.04$$) (Fig. 5B, Table 1). This significant variance was driven by CV1 (Roy”s Greatest Root $$F_{\mathrm{4,6}} = 6.05$$, $$P = 0.03$$), which is loaded most heavily by PC2 (0.489; Supplementary Table S10 available on Dryad). Multivariate contrasts revealed that only L. coggeri EU and S differ significantly in shape ($$F_{\mathrm{4,5}} = 5.00$$, $$P = 0.05$$) (Table 1). Neither L. coggeri EU and N/C ($$F_{\mathrm{4,5}} = 3.35$$, $$P = 0.11$$) nor L. coggeri N/C and S ($$F_{\mathrm{4,5}} = 2.56$$, $$P = 0.17$$) differed in shape. Therefore, the only detectible difference was that L. coggeri EU has a relatively longer body and shorter legs than the geographically adjacent, but distantly related, L. coggeri S. Table 1. Testing for morphological differences in body size and shape between the major lineages in each species group $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * Body size was tested using nested ANOVA. Body shape was tested using nested MANOVA (‘lineage’ overall effect and planned contrasts), and using the Wilks’ Lambda as the $$F$$-statistic. Table 1. Testing for morphological differences in body size and shape between the major lineages in each species group $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * $$F$$ d.f. numerator d.f. denominator $$P$$-value Sig. ‘C. rubrigularis’ N vs. S Body size (PC1) 1.17 1 50 0.29 Body shape (PC2–PC5) 0.45 4 5 0.77 ‘L. robertsi’ CU/TU vs. BK/BFAU Body size (PC1) 3.89 1 28 0.06 Body shape (PC2–PC5) 2.68 3 1 0.42 ‘L. coggeri’ N/C vs. S vs EU Body size (PC1) 2.46 2 58 0.09 Body shape (PC2–PC5) 3.46 8 10 0.04 * Multivariate contrasts N/C vs. S 2.56 4 5 0.17 N/C vs. EU 3.35 4 5 0.11 S vs. EU 5.00 4 5 0.05 * Body size was tested using nested ANOVA. Body shape was tested using nested MANOVA (‘lineage’ overall effect and planned contrasts), and using the Wilks’ Lambda as the $$F$$-statistic. Figure 5. View largeDownload slide Scatter plots of PC1 (representing body size) and PC2 (the primary axis of body shape) in each of the species groups: (A) ‘Carlia rubrigularis’ group, (B) ‘Lampropholis coggeri’ group, (C) ‘Lampropholis robertsi’ group. See Supplementary Tables S7–S9 available on Dryad for details on loadings of PC2. Morphological divergence among the lineages in each species group is limited, with the exception of some shape divergence in L. coggeri EU. Figure 5. View largeDownload slide Scatter plots of PC1 (representing body size) and PC2 (the primary axis of body shape) in each of the species groups: (A) ‘Carlia rubrigularis’ group, (B) ‘Lampropholis coggeri’ group, (C) ‘Lampropholis robertsi’ group. See Supplementary Tables S7–S9 available on Dryad for details on loadings of PC2. Morphological divergence among the lineages in each species group is limited, with the exception of some shape divergence in L. coggeri EU. No morphological differences were detected between the N and S lineages of ‘C. rubrigularis’, for either body size (‘lineage’ effect on PC1, $$F_{\mathrm{1,50}} = 1.17$$, $$P = 0.29$$) or body shape (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{4,5}} = 0.45$$, $$P = 0.77$$) (Fig. 5A, Table 1). Similarly, no significant differences were detected between the TU/CU and BK/BFAU lineages of ‘L. robertsi’. Body size was marginally non-significant (‘lineage’ effect on PC1, $$F_{\mathrm{1,28}} = 3.89$$, $$P = 0.06$$), and body shape did not differ (‘lineage’ effect on PC2–5, Wilks” Lambda $$F_{\mathrm{3,1}} = 2.68$$, $$P = 0.42$$) (Fig. 5C, Table 1). Discussion Delimitation of Cryptic Species Initial phylogeographic explorations based on mtDNA revealed that each species group contained at least four to five lineages, most of which were deeply divergent. Subsequent sequencing of five to ten nuclear loci confirmed that these phylogeographic lineages were also diverged at the nuclear genome (Dolman and Moritz 2006; Bell et al. 2010). Now, genetic data based on over 2500 exons confirmed that these lineages exhibit genetic divergences of substantial but varying depths. These genetic divergences all fall within the range that comparative data suggest spans the transition from populations to isolated species (Roux et al. 2016). Although some of the lineages are far more divergent than some already recognized species, and although we focused on morphological traits standardly used in lizard taxonomy and eco-evolutionary studies (Ingram 1991; Losos 2011; Hoskin 2014), we found little or no morphological divergence between the major lineages within each of the three species groups. Accordingly, these can mostly be regarded as truly cryptic, rather than pseudo-cryptic, species. Given morphologically cryptic lineages that span a range of divergences, and all of which are delimited using coalescent methods, we are thus faced with the two challenges of species delimitation: how to determine how much genetic divergence is sufficient when divergences are arrayed on a continuum; and how to reconcile when genetic and phenotypic data give conflicting perspectives. As a first step, we can directly assess levels of isolation between these lineages because three of the lineage-pairs in these groups meet in narrow zones of parapatry (Fig. 4). Through these fine-scale contact zone analyses of isolation, we have two major findings. First, we find that, like genetic divergence, RI exists on a continuum, with lineages exhibiting varying degrees of isolation (Fig. 4A, B). For this set of lineages, divergence and isolation appear to scale, although non-linearly. The average cline width at the L. coggeri C and S hybrid zone zone is about 5.5$$\times$$ less than that at the L. coggeri N and C contact. Yet, L. coggeri C and S is only 1.6$$\times$$ more genetically divergent than L. coggeri N and C. Theory does not predict a linear scaling; the cline width of a locus is proportional to the inverse square root of selection on that locus (Barton and Gale 1993). Thus, as selection on a locus increases, cline width can sharpen narrowly and quickly, as seen here. Further evidence of this non-linear accumulation of RI can be seen in patterns of linkage disequilibrium in these contact zones. Data from L. coggeri N and C show no evidence for disequilibrium at introgressing sites, whereas C. rubrigularis N and S and L. coggeri C and S exhibit extensive disequilibrium extending a few kilobases (Fig. 4C). These results confirm theoretical expectations that lineages can quickly transition from acting as populations (i.e., L. coggeri N and C) to acting as genomically isolated species (i.e., L. coggeri C and S and C. rubrigularis N and S) (Turner 1967; Barton 1983). Second, these data show that, despite being nearly identical morphologically, the more genetically divergent lineages have substantial RI. At least for C. rubrigularis N and S, these lineages are not isolated by premating isolation (Dolman 2008), but rather by post-mating selection against hybrids (Phillips et al. 2004). Based on estimates of dispersal length and cline width of the hybrid index, selection against hybrids is strong. Hybrids between C. rubrigularis N and S and L. coggeri C and S are estimated to be 50–70% and 10–65% less fit than their parents, respectively (Phillips et al. 2004; Singhal and Moritz 2012). Estimated selection on hybrids between L. coggeri N and C, on the other hand, is negligible. The selection against hybrids seen in C. rubrigularis N and S and L. coggeri C and S is comparable (if not greater) than that seen between morphologically distinct hybridizing taxon-pairs. (Barton and Gale 1993; Singhal and Moritz 2012). Such strong selection suggests that these lineages will remain evolutionary distinct in the future, despite the high potential for gene flow. As such, we propose to identify C. rubrigularis N and S as separate species, and likewise for L. coggeri S and L. coggeri C. Because we found no evidence for RI between L. coggeri N and C, we retain them as distinct populations within one species (L. coggeri N/C). However, how should we diagnose those lineages for which we cannot indirectly or directly assay RI? For example, L. coggeri EU is geographically isolated from L. coggeri N/C and S, and the lineages in the “L. robertsi” group are isolated on different mountaintops. These lineages do not meet in contact zones, and because of both practical and ethical reasons, cannot be easily kept in the laboratory for experimental trials. Instead, we extrapolate our estimates of RI from the three lineage-pairs that do meet in parapatry (Fig. 4) to the species group as a whole. This extrapolation assumes that the tempo and mode at which RI evolves is similar across this clade. The few comparative data on the rate at which RI evolves suggests that it can vary across broad clades (Rabosky and Matute 2013). However, for a clade like this, which consists of broadly related, morphologically and ecologically-similar lizards found in a similar biogeographic context, we suspect there is likely to be less variation. Indeed, RI and divergence time correlate closely across five sister-species comparisons in Carlia, Lampropholis, and a closely-related genus, Saproscincus (Singhal and Moritz 2013; Singhal and Bi 2017). Further, and importantly, these lineages likely resulted from similar speciation processes—i.e., these deep, cryptic lineages evolved due to very long periods of isolation in environmentally similar refugia (see below). Therefore, we believe we can sensibly extrapolate our results across other cryptic congeneric lineages that diverged under similar processes. To extrapolate, we use the divergence between C. rubrigularis N and S as our cutoff because they are the youngest lineage-pair for which we have solid evidence of RI. Divergence estimates are highly correlated across the three metrics for genetic divergence ($$r \sim 0.7-9$$; Fig. 3). Still, we take a conservative approach and only elevate those lineages that show greater divergence than what is seen for C. rubrigularis N and S across all metrics (with one exception, below). Notably, this cutoff is greater than that seen among several comparisons between nominal taxa (Fig. 3). Divergences between L. robertsi CU and TU, L. robertsi BK and BFAU, and C. rhomboidalis N and S all fall below this cutoff (Supplementary Table S5 available on Dryad), so we recognize these as phylogeographic lineages within species. However, the divergence between L. robertsi BK/BFAU and L. robertsi CU/TU is greater than this cutoff. Accordingly, we propose to diagnose the BK/BFAU and CU/TU allopatric lineages as species; this deep divergence suggests they are likely to exhibit RI should they ever come into contact. We use these groupings for morphological comparisons, as well. The case of L. coggeri EU is more ambiguous, when considering genetic data alone. Lampropholis coggeri EU is most closely related to L. coggeri N/C, and divergence falls just below our proposed cutoff for raw divergence and branching time (Fig. 3A, C) but just above the cutoff for net divergence (Fig. 3B). However, L. coggeri EU sits isolated off the far south end of the Wet Tropics, geographically closest to L. coggeri S. Therefore, it is much more likely to interact with L. coggeri S in future, with which it shows much greater divergence (Fig. 1, Fig. 3; Supplementary Table S5 available on Dryad). Further, alone among the comparisons made here, L. coggeri EU is morphologically (and perhaps ecologically; see below) distinct from L. coggeri S (Table 1). Because the L. coggeri S and EU comparison is more salient than the L. coggeri C and EU comparison, we further propose to elevate L. coggeri EU as a separate species. Taxonomic Outcomes The formal taxonomic revisions of these lineages are presented in Appendix 1 available on Dryad. To summarize, we revise the ‘C. rubrigularis’ group to retain C. wundalthini and C. rhomboidalis, retain C. rubrigularis S as C. rubrigularis, and elevate C. rubrigularis N to C. crypta sp. nov. For the ‘L. coggeri’ group, we retain L. coggeri N/C as L. coggeri, elevate L. coggeri S to L. similis sp. nov., and elevate L. coggeri EU to L. elliotensis sp. nov. For the ‘L. robertsi’ group, we retain L. robertsi CU/TU as L. robertsi and elevate L. robertsi BFAU/BK to L. bellendenkerensis sp. nov. The key components of the new species descriptions are as follows. Carlia crypta sp. nov. Holotype: QM J75457 Diagnosis: Carlia crypta sp. nov. is distinguished from all other Carlia spp., except members of the ‘C. rhomboidalis’ group, in possessing an interparietal fused to the frontoparietals. As with C. rubrigularis, adult males possess a red throat. It is reliably distinguished from this species by four nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in three amino acid differences (Table A3). Zoobank ID: 3B116172-035D-41EF-958E-F0A35BB15F0D Lampropholis similis sp. nov. Holotype: QM J91380 Diagnosis: Lampropholis similis sp. nov. is a small, dark-sided rainforest skink with pentadactyl limbs (overlapping or very narrowly separated when adpressed) and a movable lower eyelid containing a transparent disc. It is reliably distinguished from its sibling species (L. coggeri and L. elliotensis sp. nov.) by 17 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in 15 amino acid differences (Table A1). Zoobank ID: 23E36086-C261-4D2A-958B-547A17239529 Lampropholis elliotensis sp. nov. Holotype: QM J91382 Diagnosis: Lampropholis elliotensis sp. nov. is a small, dark-sided rainforest skink with pentadactyl limbs (usually separated by several scales rows when adpressed) and a movable lower eyelid containing a transparent disc. It is reliably distinguished from its sibling species (L. similis sp. nov. and L. coggeri) by 17 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in 15 amino acid differences among the species (Table A1). Zoobank ID: 6CE54452-528E-4B3E-9655-AC2262A26849 Lampropholis bellendenkerensis sp. nov. Holotype: QM J39855 Diagnosis: A large Lampropholis with dark flanks and prominent spotting on the posterior ventral surfaces, a row of dark edged pale spots on underside of tail. This species is reliably distinguished from its closest congener (L. robertsi) by 13 nucleotide differences in the mitochondrial gene NADH dehydrogenase 4 that result in nine amino acid differences between the species (Table A2). Zoobank ID: C721A396-F4A2-4205-A7FE-ED8FC8F 90FAD Speciation Processes and Cryptic Species Extrapolating a divergence cutoff as done here works, in part, because the Wet Tropics lineages likely share similar divergence histories and processes. These taxa diverged through extended periods of isolation in environmentally-similar, climatically-stable rainforest refugia, resulting in genetic divergence but eco-morphological conservatism (Graham et al. 2006; Moritz et al. 2009). The widespread Wet Tropics lineages (C. rubrigularis N and S; L. coggeri N/C and S) occupy similar habitats, and the lineages of ‘L. robertsi’ occupy mountaintop habitats that appear broadly similar. The morphological stasis in these cryptic species therefore likely reflects a lack of divergent selection across similar environments (Moritz et al. 2009; Hoskin et al. 2011). Only those lineages peripheral to the main block of the Wet Tropics rainforests exhibit morphological divergence. For example, L. coggeri EU differs subtly in body shape from other species in the group; it also occupies a subset of the habitats occupied by the other lineages in the ‘L. coggeri’ group. Whereas the other lineages are found across a broad range of rainforest types and habitat edges, L. coggeri EU is found only in mesic pockets of rocky, upland rainforest. Greater morphological divergence is seen in the two lineages in the ‘C. rubrigularis’ group found outside the Wet Tropics region: C. wundalthini and C. rhomboidalis. Despite being of similar age as other lineages (Fig. 2), these species are distinct for male breeding coloration and some aspects of body size and shape. As seen across other Wet Tropics taxa, phenotypic divergence is only recovered in association with environmentally-driven selection – e.g., across ecotones (Schneider and Moritz 1999), in peripheral isolates (Hoskin et al. 2011), or reinforcing selection in hybrid zones (Hoskin et al. 2005). These cryptic species thus underline the importance of “non-ecological” mechanisms in speciation (Schluter 2001), in particular, mutation-order speciation (Mani and Clarke 1990; Nosil and Flaxman 2011). This work also helps further outline the distinction between historical and current patterns of hybridization and introgression (Edwards et al. 2016). The biological species concept requires current introgression to be limited, but both theoretical models and empirical data show that species can diverge and remain distinct in the presence of historical gene flow (Nosil 2008; Pinho and Hey 2010). Similarly, our tests for introgression suggest there has been introgression between at least four lineage-pairs (Fig. 2), all of which are either previously-recognized or now elevated nominal species. Yet, for at least two of these lineage-pairs (C. rubrigularis N and S and L. coggeri C and S), our analysis of current patterns shows that hybridization occurs but is geographically highly-restricted (Fig. 4). This distinction between historical and current patterns illustrates the complicated relationship between gene flow and species borders. The Practicality of Delimiting Cryptic Species Genetic divergences for almost all our lineage comparisons fall in the so-called “gray zone” of speciation, in which lineages transition from behaving as populations to species (Roux et al. 2016). Defined by net silent divergence ($$d_{a})$$ at coding nuclear genes, this gray zone spans divergences from 0.5% to 2%. Given that many cryptic lineages originated during glacial cycles over the last few million years (Hewitt 2000), many of them should fall within this four-fold range of divergence. This underlines the challenge in delimiting cryptic lineages—many of them have a biogeographic and divergence history that places them in an ambiguous zone of divergence, where lineages are as likely to merge or remain distinct upon secondary contact. Given this ambiguity, identifying strong phylogeographic structure within species should be just the first step in diagnosing species boundaries across cryptic boundaries (Fig. 6, Table 2). Additional validation is required, which we loosely group into four categories: 1) statistical species delimitation, 2) post hoc discovery of phenotypic differences that delimit lineages (i.e., integrative taxonomy; Padial et al. 2010), 3) indirect or direct estimates of evolutionary isolation between lineages, or 4) calibration-based approaches. Herein, we outline these approaches briefly and explain their conceptual and practical benefits and limitations. Note that after applying any of these approaches, researchers must still formally revise the taxonomy for these validations to be recognized. However, all too often, these taxonomic revisions are not done (Carstens et al. 2013; Pante et al. 2014). Table 2. A survey of approaches that can be used to validate putative cryptic lineages and the benefits and limitations of each approach Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). We highlight, which approaches were used in this study and identify other examples from the broader literature; this list of examples is meant to be illustrative not exhaustive. Many of the studies that validated putative cryptic species did not also conduct a formal taxonomic revision. Table 2. A survey of approaches that can be used to validate putative cryptic lineages and the benefits and limitations of each approach Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). Approaches to validate putative cryptic species Benefits Limitations Examples from this study Other examples Statistical species delimitation Efficient and affordable; can be applied to asexual and sexual organisms; many methods can handle ancestral polymorphism. Populations cannot be easily distinguished from true species (Sukumaran and Knowles 2017); different approaches often lead to differing results (Carstens et al. 2013); results can be deceiving in the presence of gene-flow Geckos (Leaché and Fujita 2010), carnivorous plants (Carstens and Satler 2013), cavefish (Niemiller et al. 2012), mouse lemurs (Hotaling et al. 2016). Identification of morphological, behavioral, physiological differences among ‘pseudo’ cryptic lineages Can lead to the identification of divergence in traits that are likely to keep lineages distinct (e.g., ecological differences, phenological shifts, mating calls). Phenotypic differences do not guarantee that lineages will remain distinct if they interact, (e.g., Hoskin et al. 2005). L. coggeri EU lizards (Silva et al. 2017), hoverflies (Schonrogge et al. 2002), mayflies (Leys et al. 2017), wasps (Smith et al. 2008), butterflies (Hebert et al. 2004a), salamanders (Rissler and Apodaca 2007), frogs (Hoskin et al. 2011). Identification of range overlap among cryptic species Offers robust evidence that lineages are not interbreeding. Many cryptic species are parapatric or allopatric, so not applicable to many taxa. Rotifers (Gómez et al. 2002), frogs (Stuart et al. 2006), Plasmodium (Perkins 2000), earthworms (King et al. 2008), fish (Adams et al. 2014). Direct estimates of reproductive isolation through mate choice studies and crossing experiments Offers robust evidence that lineages are not interbreeding. Can only be used for sexually reproducing species that are amenable to lab husbandry; expensive and time-consuming; lab-based estimates of mate choice and hybrid fitness can differ from field-based estimates. Rotifers (Suatoni et al. 2006), bryozoans (Gomez et al. 2007), diatoms (Amato et al. 2007), frogs (Hoskin et al. 2005). Indirect estimates of reproductive isolation from regions of parapatry or sympatry Allows indirect measure of factors structuring extent of gene flow between lineages such as extent of assortative mating, genetic incompatibilities, etc. Can only be used for sexually reproducing species that co-occur and exist at sufficient density for sampling; indirect estimates can be influenced by often uncharacterized demographic factors C. rubrigularis N, C. rubrigularis S, L. coggeri C, L. coggeri S Frogs (Stewart and Lougheed 2013), lizards (Sites et al. 1995), kelp (Tellier et al. 2011), corals (Ladner and Palumbi 2012). Using calibrations of RI vs. genomic divergence based on data from closely-related species Provides a well-informed guideline for species that are likely to evolve reproductive isolation at a similar tempo. Still requires the extensive and expensive collection of data on reproductive isolation in closely-related species. L. robertsi TU/CU, L. robertsi BK/ BFAU Using calibrations based on sequence divergence Efficient and affordable; can be applied to asexual and sexual organisms; can be more clade-specific if informed by patterns of divergence between nominal species in the clade. Not conceptually well-grounded (but see (Roux et al. 2016)); what metric of divergence to use is unclear (Fregin et al. 2012), clades might vary in the rate at which they evolve reproductive barriers (Rabosky and Matute 2013). Caddisflies (Bálint et al. 2011), platyhelminth (Vilas et al. 2005), frogs (Fouquet et al. 2007), amphipods (Witt et al. 2006), birds (Olsson et al. 2005), parasitoid flies (Smith et al. 2006). We highlight, which approaches were used in this study and identify other examples from the broader literature; this list of examples is meant to be illustrative not exhaustive. Many of the studies that validated putative cryptic species did not also conduct a formal taxonomic revision. Figure 6. View largeDownload slide A flowchart outlining a possible research approach to validating cryptic lineages. Statistical species delimitation is often used to both define putative lineages and to validate them. Width and contrast of validation arrows indicate how robust we believe these approaches are. Benefits and limitations of each approach are further described in Table 2. Figure 6. View largeDownload slide A flowchart outlining a possible research approach to validating cryptic lineages. Statistical species delimitation is often used to both define putative lineages and to validate them. Width and contrast of validation arrows indicate how robust we believe these approaches are. Benefits and limitations of each approach are further described in Table 2. First, perhaps the most common approach currently used is statistical species delimitation, which applies coalescent-based methods to determine which lineages are genetically unique (Ence and Carstens 2011; Yang and Rannala 2014). Often, statistical approaches are used as an early step to more quantitatively assess visual clusters, and in other studies, they are used as a final stage of analysis to diagnose species (Fig. 6). Across our three species groups, statistical species delimitation diagnosed all putative lineages as species (Fig. 2). Yet, this statistical diagnosis contrasts to our biological understanding of species boundaries. For example, although L. coggeri N and C are statistically distinct, introgression between them is widespread genomically and geographically (Fig. 4). This disconnect reflects the emerging consensus that, although statistical species delimitation methods can robustly identify populations, these populations are not always equivalent to species (Rosenblum et al. 2012; Dynesius and Jansson 2014; Sukumaran and Knowles 2017). In other words, the genetic distinctiveness of a population does not necessarily confer robust evolutionary distinctiveness as envisaged under the BSC. Thus, we take a deliberately conservative approach and refrain from elevating morphologically cryptic lineages whose sole support is from statistical analyses of genetic data (Oliver et al. 2015). However, statistical approaches to species delimitation are the most efficient and flexible method across taxonomic groups and are likely to remain an attractive option for many study systems, even if these approaches alone provide insufficient evidence to denote robust taxa. Second, while mostly not true in the present study (Fig. 5), researchers often discover post-hoc phenotypic differences after further investigation of putative cryptic lineages, leading to so-called pseudo-cryptic species (Knowlton 1993). For example, cryptic lineages might vary in traits that facilitate RI between lineages (e.g., mating calls, (Barber 1951)) or ecological co-existence (e.g., divergence in life-history; Leys et al. 2017). Often, however, these phenotypic differences might only relate trivially to how distinctive a lineage is—i.e., minor differences in scale counts. In such cases, phenotypic differentiation alone might not necessarily confer evolutionary distinctiveness, and further validation would be required. Third, researchers can assay strength of RI between lineages either through indirect studies of contact zones, observation of sympatry among cryptic lineages, or laboratory- or field-based tests of mate choice and hybrid fitness (Blair 1972; Hoskin et al. 2005; Dolman 2008). As shown in this work (Fig. 4), these approaches offer detailed data on how likely lineages are to remain distinct if and when they overlap with their congenerics. However, these approaches are often quite practically limited—generating these data can require contiguous ranges, high population densities, organisms amenable to experimentation, and substantial investment in both time and money. These practical limitations surfaced in the present work; we were unable to test for RI between allopatric lineages, nor could we bring them into a laboratory setting. Barring this, alternative approaches could be used to identify cases of abrupt genetic boundaries across dense sampling or test for geographically-extensive introgression across lineage boundaries using large numbers of markers (Melville et al. 2017). For those lineages not amenable to indirect or direct testing for RI, we used a fourth general approach: using calibrations to determine how much divergence is sufficient to elevate a species. In the DNA barcoding literature, these calibrations are typically informed by patterns of divergence among nominal taxa, although they are used across broad swaths of the tree of life—i.e., all birds or all butterflies (Hebert et al. 2004b; Janzen et al. 2005). Another option is to use calibrations informed by data on the tempo at which RI evolves in closely-related taxa, as done in this study. While this approach still requires the expensive and time-consuming generation of data on isolation between lineages, the cutoff is more principled than barcode gaps. Unlike barcode gaps, which are applied widely across taxonomic groups, our cutoff is based on the observed isolation between lineages within a given clade, which likely share a common mode of lineage divergence and speciation. Such an approach allows us to tackle cryptic diversity while reflecting the variable nature of the speciation process. Moving forward, we will explore whether this cutoff could contribute to diagnosing species boundaries among phylogeographic lineages of other species of Carlia (Potter et al. 2016). Across all these approaches, a potential flaw is insufficient sampling of species ranges. In particular, by sampling two ends of an array of populations, we can infer distinction between lineages where gene flow is actually continuous throughout (Pante et al. 2015). Because all approaches to cryptic species validation either have conceptual or practical weaknesses (Table 2), the ideal approach will likely vary across taxonomic groups. For example, in the Wet Tropics, phylogeographic studies have recovered deep splits in other taxa outside of lizards, including frogs and mammals (Moritz et al. 2009). For the frog lineages, many of which occur in dense numbers, meet in narrow contact zones, and are amenable to breeding experiments, data on hybridization patterns and mate choice have helped validate putative cryptic lineages and led to the formal revision of lineages (Hoskin et al. 2005; Hoskin 2007). However, for mammals, density is too low to allow indirect or direct tests of isolation, so other approaches will be required. Importantly, the framework we applied here—generating initial descriptions of within-species phylogeographic diversity (Dolman and Moritz 2006; Bell et al. 2010), confirming these patterns with multilocus data sets, and then inferring fine-scale patterns of current RI (Phillips et al. 2004; Singhal and Moritz 2012, 2013; Singhal and Bi 2017)—represents a significant investment in time and money. For many systems, this approach is simply not tenable, and further, often too slow where decisions on species limits have immediate conservation consequences (Hedin 2015). That said, given the dangers of ‘taxonomic overinflation’ (Isaac et al. 2004; Frankham et al. 2012), we advocate that researchers validate putative cryptic lineages by both considering statistical delimitation approaches and other data on the biological reality of lineages, whether that be direct or indirect evidence for isolation. Conclusion Cryptic species challenge traditional notions of species, because the discrepancy between morphological and genetic axes of divergence can make them hard to categorize. Yet, other data suggest that many cryptic species are phenotypically divergent, but on axes of variation that are harder to measure (e.g., mating pheromones in lizards). In cases where we cannot identify phenotypic differences, like the lizards of the Wet Tropics, we can test the validity of these lineages through other means, such as looking at interactions between cryptic lineages in parapatry. Often, these richer, more integrative datasets complement genetic data and show that cryptic lineages are independently evolving units. However, as we also see in these taxa, despite marked genetic differentiation, some cryptic lineages might just be ephemera, destined to be lost to hybridization with congeneric lineages, if they meet in the future. These data remind us that species boundaries are hypotheses (Fujita et al. 2012), our best estimate of the fate of these lineages and a recognition of the ever-evolving nature of species (Darwin 1859; Mallet 1995). Acknowledgments We thank J. Bragg and also staff at UMich Arc TS Flux for advice and logistical support. In addition, we gratefully acknowledge the members of the C. Moritz & S. Williams research groups, both past and present, who have contributed samples, environmental data, and their expertise to the Wet Tropics research program. Funding Funding for this project came from a National Science Foundation Post-doctoral Research Fellowship in Biology (NSF #1519732 to SS), grants from the Australian Research Council (ARC), the Australian Biological Resources Study (ABRS) and the National Science Foundation to CM, and grants from the ARC and ABRS to CH. Data Accessibility Raw sequencing data: PRJNA448788 Pseudo-reference files and variant sets: DataDryad accession 10.5061/dryad.g7v1b. Code used to generate data: https://github.com/singhal/AWT_delimit and https://github.com/singhal/SqCL. Specimens used for morphological analyses & species descriptions: all are deposited at Queensland Museum for research use; accession numbers available in Supplementary Appendix 1 available on Dryad. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.g7v1b. References Adams M. , Raadik T.A. , Burridge C.P. , Georges A. 2014 . Global biodiversity assessment and hyper-cryptic species complexes: more than one species of elephant in the room? Syst. Biol. 63 : 518 – 533 . Google Scholar CrossRef Search ADS PubMed Amato A. , Kooistra W.H. , Ghiron J.H.L. , Mann D.G. , Pröschold T. , Montresor M. 2007 . Reproductive isolation among sympatric cryptic species in marine diatoms . Protist 158 : 193 – 207 . Google Scholar CrossRef Search ADS PubMed Bálint M. , Domisch S. , Engelhardt C. , Haase P. , Lehrian S. , Sauer J. , Theissinger K. , Pauls S. , Nowak C. 2011 . Cryptic biodiversity loss linked to global climate change . Nat. Clim. Change 1 : 313 – 318 . Google Scholar CrossRef Search ADS Barber H.S. 1951 . North American fireflies of the genus Photuris . Smithson. Misc. Collect. 117 : 1 – 58 . Barton N. 1983 . Multilocus clines . Evolution 37 : 454 – 471 . Google Scholar CrossRef Search ADS PubMed Barton N.H. , Gale K.S. 1993 . Genetic analysis of hybrid zones. In: Hybrid zones and the evolutionary process ., (ed. Harrison RG ), pp. 13 – 45 . New York : Oxford University Press . Bell R.C. , Parra J.L. , Tonione M. , Hoskin C.J. , Mackenzie J.B. , Williams S.E. , Moritz C. 2010 . Patterns of persistence and isolation indicate resilience to climate change in montane rainforest lizards . Mol. Ecol. 19 : 2531 – 2544 . Google Scholar PubMed Bickford D. , Lohman D.J. , Sodhi N.S. , Ng P.K. , Meier R. , Winker K. , Ingram K.K. , Das I. 2007 . Cryptic species as a window on diversity and conservation . Trends Ecol. Evol. 22 : 148 – 155 . Google Scholar CrossRef Search ADS PubMed Blair W.F. 1972 . Evolution in the genus Bufo . Austin (TX) : University of Texas Press . Bolger A.M. , Lohse M. , Usadel B. 2014 . Trimmomatic: a flexible trimmer for Illumina sequence data . Bioinformatics 30 : 2114 – 2120 . Google Scholar CrossRef Search ADS PubMed Bolnick D.I. , Near T.J. , Wainwright P.C. 2006 . Body size divergence promotes post-zygotic reproductive isolation in centrarchids . Evol. Ecol. Res. 8 : 903 – 913 . Bragg J.G. , Potter S. , Afonso Silva A.C. , Hoskin C.J. , Bai B.Y.H. , Moritz C. 2018 . Phylogenomics of a rapid radiation: the Australian rainbow skinks . BMC Evol. Biol. 18 : 15 . Google Scholar CrossRef Search ADS PubMed Bragg J.G. , Potter S. , Bi K. , Moritz C. 2016 . Exon capture phylogenomics: efficacy across scales of divergence . Mol. Ecol. Resour. 16 : 1059 – 1068 . Google Scholar CrossRef Search ADS PubMed Brown J.M. , Hedtke S.M. , Lemmon A.R. , Lemmon E.M. 2009 . When trees grow too long: investigating the causes of highly inaccurate Bayesian branch-length estimates . Syst. Biol. 59 : 145 – 161 . Google Scholar CrossRef Search ADS PubMed Carstens B.C. , Pelletier T.A. , Reid N.M. , Satler J.D. 2013 . How to fail at species delimitation . Mol. Ecol. 22 : 4369 – 4383 . Google Scholar CrossRef Search ADS PubMed Carstens B.C. , Satler J.D. 2013 . The carnivorous plant described as Sarracenia alata contains two cryptic species . Biol. J. Linnean Soc. 109 : 737 – 746 . Google Scholar CrossRef Search ADS Coyne J.A. , Orr H.A. 1989 . Patterns of speciation in Drosophila . Evolution 43 : 362 – 381 . Google Scholar CrossRef Search ADS PubMed Darwin C. 1859 . The origin of species by means of natural selection: or, the preservation of favored races in the struggle for life . London, UK : John Murray . de León G.P.-P. , Poulin R. 2016 . Taxonomic distribution of cryptic diversity among metazoans: not so homogeneous after all . Biol. Lett. 12 : 20160371 . Google Scholar CrossRef Search ADS PubMed De Queiroz K. 2007 . Species concepts and species delimitation . Syst. Biol. 56 : 879 – 886 . Google Scholar CrossRef Search ADS PubMed Dolman G. 2008 . Evidence for differential assortative female preference in association with refugial isolation of rainbow skinks in Australia”s tropical rainforests . PLoS One 3 : e3499 . Google Scholar CrossRef Search ADS PubMed Dolman G. , Moritz C. 2006 . A multilocus perspective on refugial isolation and divergence in rainforest skinks (Carlia) . Evolution 60 : 573 – 582 . Google Scholar CrossRef Search ADS PubMed Durand E.Y. , Patterson N. , Reich D. , Slatkin M. 2011 . Testing for ancient admixture between closely related populations . Mol. Biol. Evol. 28 : 2239 – 2252 . Google Scholar CrossRef Search ADS PubMed Dynesius M. , Jansson R. 2014 . Persistence of within species lineages: a neglected control of speciation rates . Evolution 68 : 923 – 934 . Google Scholar CrossRef Search ADS PubMed Eaton D.A. , Ree R.H. 2013 . Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae) . Syst. Biol. 62 : 689 – 706 . Google Scholar CrossRef Search ADS PubMed Edwards S.V. , Potter S. , Schmitt C.J. , Bragg J.G. , Moritz C. 2016 . Reticulation, divergence, and the phylogeography–phylogenetics continuum . Proc. Natl. Acad. Sci. U.S.A 113 : 8025 – 8032 . Google Scholar CrossRef Search ADS PubMed Ence D.D. , Carstens B.C. 2011 . SpedeSTEM: a rapid and accurate method for species delimitation . Mol. Ecol. Resour. 11 : 473 – 480 . Google Scholar CrossRef Search ADS PubMed Fitzpatrick B.M. 2002 . Molecular correlates of reproductive isolation . Evolution 56 : 191 – 198 . Google Scholar CrossRef Search ADS PubMed Fouquet A. , Gilles A. , Vences M. , Marty C. , Blanc M. , Gemmell N.J. 2007 . Underestimation of species richness in Neotropical frogs revealed by mtDNA analyses . PLoS One 2 : e1109 . Google Scholar CrossRef Search ADS PubMed Frankham R. , Ballou J.D. , Dudash M.R. , Eldridge M.D. , Fenster C.B. , Lacy R.C. , Mendelson J.R. , Porton I.J. , Ralls K. , Ryder O.A. 2012 . Implications of different species concepts for conserving biodiversity . Biol. Conserv. 153 : 25 – 31 . Google Scholar CrossRef Search ADS Fregin S. , Haase M. , Olsson U. , Alström P. 2012 . Pitfalls in comparisons of genetic distances: a case study of the avian family Acrocephalidae . Mol. Phylogenet. Evol. 62 : 319 – 328 . Google Scholar CrossRef Search ADS PubMed Fujita M.K. , Leaché A.D. , Burbrink F.T. , McGuire J.A. , Moritz C. 2012 . Coalescent-based species delimitation in an integrative taxonomy . Trends Ecol. Evol. 27 : 480 – 488 . Google Scholar CrossRef Search ADS PubMed Funk D.J. , Nosil P. , Etges W.J. 2006 . Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa . Proc. Natl. Acad. Sci. U.S.A 103 : 3209 – 3213 . Google Scholar CrossRef Search ADS PubMed Gómez A. , Serra M. , Carvalho G.R. , Lunt D.H. 2002 . Speciation in ancient cryptic species complexes: evidence from the molecular phylogeny of Brachionus plicatilis (Rotifera) . Evolution 56 : 1431 – 1444 . Google Scholar CrossRef Search ADS PubMed Gomez A. , Wright P.J. , Lunt D.H. , Cancino J.M. , Carvalho G.R. , Hughes R.N. 2007 . Mating trials validate the use of DNA barcoding to reveal cryptic speciation of a marine bryozoan taxon . Proc. R. Soc. Lond. B Biol. Sci. 274 : 199 – 207 . Google Scholar CrossRef Search ADS Grabherr M.G. , Haas B.J. , Yassour M. , Levin J.Z. , Thompson D.A. , Amit I. , Adiconis X. , Fan L. , Raychowdhury R. , Zeng Q. 2011 . Full-length transcriptome assembly from RNA-Seq data without a reference genome . Nat. Biotechnol. 29 : 644 – 652 . Google Scholar CrossRef Search ADS PubMed Graham C.H. , Moritz C. , Williams S.E. 2006 . Habitat history improves prediction of biodiversity in rainforest fauna . Proc. Natl. Acad. Sci. U.S.A 103 : 632 – 636 . Google Scholar CrossRef Search ADS PubMed Greer A. 1997 . A new species of Lampropholis (Squamata: Scincidae) with a restricted, high altitude distribution in eastern Australia . Aust. Zool. 30 : 360 – 368 . Google Scholar CrossRef Search ADS Hebert P.D. , Penton E.H. , Burns J.M. , Janzen D.H. , Hallwachs W. 2004a . Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly Astraptes fulgerator . Proc. Natl. Acad. Sci. U.S.A 101 : 14812 – 14817 . Google Scholar CrossRef Search ADS Hebert P.D. , Stoeckle M.Y. , Zemlak T.S. , Francis C.M. 2004b . Identification of birds through DNA barcodes . PLoS Biol. 2 : e312 . Google Scholar CrossRef Search ADS Hedgecock D. , Ayala F.J. 1974 . Evolutionary divergence in the genus Taricha (Salamandridae) . Copeia 3 : 738 – 747 . Google Scholar CrossRef Search ADS Hedin M. 2015 . High stakes species delimitation in eyeless cave spiders (Cicurina, Dictynidae, Araneae) from central Texas . Mol. Ecol. 24 : 346 – 361 . Google Scholar CrossRef Search ADS PubMed Hewitt G. 2000 . The genetic legacy of the Quaternary ice ages . Nature 405 : 907 – 913 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. 2007 . Description, biology and conservation of a new species of Australian tree frog (Amphibia: Anura: Hylidae: Litoria) and an assessment of the remaining populations of Litoria genimaculata Horst, 1883: systematic and conservation implications of an unusual speciation event . Biol. J. Linnean Soc. 91 : 549 – 563 . Google Scholar CrossRef Search ADS Hoskin C.J. 2014 . A new skink (Scincidae: Carlia) from the rainforest uplands of Cape Melville, north-east Australia . Zootaxa 3869 : 224 – 236 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. , Higgie M. , McDonald K.R. , Moritz C. 2005 . Reinforcement drives rapid allopatric speciation . Nature 437 : 1353 – 1356 . Google Scholar CrossRef Search ADS PubMed Hoskin C.J. , Tonione M. , Higgie M. , MacKenzie J.B. , Williams S.E. , VanDerWal J. , Moritz C. 2011 . Persistence in peripheral refugia promotes phenotypic divergence and speciation in a rainforest frog . Am. Nat. 178 : 561 – 578 . Google Scholar CrossRef Search ADS PubMed Hotaling S. , Foley M.E. , Lawrence N.M. , Bocanegra J. , Blanco M.B. , Rasoloarison R. , Kappeler P.M. , Barrett M.A. , Yoder A.D. , Weisrock D.W. 2016 . Species discovery and validation in a cryptic radiation of endangered primates: coalescent based species delimitation in Madagascar”s mouse lemurs . Mol. Ecol. 25 : 2029 – 2045 . Google Scholar CrossRef Search ADS PubMed Ingram G. 1991 . Five new skinks from Queensland rainforests . Mem. Queensl. Mus. 30 : 443 – 453 . Ingram G. , Covacevich J. 1989 . Revision of the genus Carlia (Reptilia, Scincidae) in Australia with comments on Carlia bicarinata of New Guinea . Mem. Queensl. Mus. 27 : 443 – 490 . Isaac N.J. , Mallet J. , Mace G.M. 2004 . Taxonomic inflation: its influence on macroecology and conservation . Trends Ecol. Evol. 19 : 464 – 469 . Google Scholar CrossRef Search ADS PubMed Janzen D.H. , Hajibabaei M. , Burns J.M. , Hallwachs W. , Remigio E. , Hebert P.D. 2005 . Wedding biodiversity inventory of a large and complex Lepidoptera fauna with DNA barcoding . Philos. Trans. R. Soc. Lond. B Biol. Sci. 360 : 1835 – 1845 . Google Scholar CrossRef Search ADS PubMed Jolliffe I.T. 2002 . Principal component analysis . New York (NY) : Springer-Verlag . Jones G. 2017 . Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent . J. Math. Biol. 74 : 447 – 467 . Google Scholar CrossRef Search ADS PubMed Katoh K. , Misawa K. , Kuma K.I. , Miyata T. 2002 . MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform . Nucleic Acids Res. 30 : 3059 – 3066 . Google Scholar CrossRef Search ADS PubMed Kent W.J. 2002 . BLAT—the BLAST-like alignment tool . Genome Res. 12 : 656 – 664 . Google Scholar CrossRef Search ADS PubMed King R.A. , Tibble A.L. , Symondson W.O. 2008 . Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms . Mol. Ecol. 17 : 4684 – 4698 . Google Scholar CrossRef Search ADS PubMed Knowlton N. 1993 . Sibling species in the sea . Ann. Rev. Ecol. Evol. Syst. 24 : 189 – 216 . Google Scholar CrossRef Search ADS Ladner J.T. , Palumbi S.R. 2012 . Extensive sympatry, cryptic diversity and introgression throughout the geographic distribution of two coral species complexes . Mol. Ecol. 21 : 2224 – 2238 . Google Scholar CrossRef Search ADS PubMed Lanfear R. , Frandsen P.B. , Wright A.M. , Senfeld T. , Calcott B. 2016 . PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses . Mol. Biol. Evol. 34 : 772 – 773 . Leaché A.D. , Fujita M.K. 2010 . Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus) . Proc. R. Soc. Lond. B Biol. Sci. 277 1697 : rspb20100662 . Leys M. , Keller I. , Robinson C.T. , Räsänen K. 2017 . Cryptic lineages of a common alpine mayfly show strong life history divergence . Mol. Ecol. 26 : 1670 – 1686 . Google Scholar CrossRef Search ADS PubMed Li H. , Durbin R. 2009 . Fast and accurate short read alignment with Burrows–Wheeler transform . Bioinformatics 25 : 1754 – 1760 . Google Scholar CrossRef Search ADS PubMed Losos J.B. 2011 . Lizards in an evolutionary tree: ecology and adaptive radiation of anoles . Berkeley (CA) : University of California Press . Mallet J. 1995 . A species definition for the modern synthesis . Trends Ecol. Evol. 10 : 294 – 299 . Google Scholar CrossRef Search ADS PubMed Mani G. , Clarke B. 1990 . Mutational order: a major stochastic process in evolution . Proc. R. Soc. Lond. B Biol. Sci. 240 : 29 – 37 . Google Scholar CrossRef Search ADS PubMed Mayr E. 1942 . Systematics and the origin of species, from the viewpoint of a zoologist . Cambridge (MA) : Harvard University Press . McDaniel S.F. , Shaw A.J. 2003 . Phylogeographic structure and cryptic speciation in the trans-Antarctic moss Pyrrhobryum mnioides . Evolution 57 : 205 – 215 . Google Scholar CrossRef Search ADS PubMed McKenna A. , Hanna M. , Banks E. , Sivachenko A. , Cibulskis K. , Kernytsky A. , Garimella K. , Altshuler D. , Gabriel S. , Daly M. 2010 . The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res. 20 : 1297 – 1303 . Google Scholar CrossRef Search ADS PubMed Melville J. , Haines M.L. , Boysen K. , Hodkinson L. , Kilian A. , Date K.L.S. , Potvin D.A. , Parris K.M. 2017 . Identifying hybridization and admixture using SNPs: application of the DArTseq platform in phylogeographic research on vertebrates . R. Soc. Open Sci. 4 : 161061 . Google Scholar CrossRef Search ADS PubMed Meyer M. , Kircher M. 2010 . Illumina sequencing library preparation for highly multiplexed target capture and sequencing . Cold Spring Harbor Protocols 2010 : pdb. prot5448 . Moritz C. , Cicero C. 2004 . DNA barcoding: promise and pitfalls . PLoS Biol. 2 : e354 . Google Scholar CrossRef Search ADS PubMed Moritz C. , Hoskin C. , MacKenzie J.B. , Phillips B. , Tonione M. , Silva N. , VanDerWal J. , Williams S.E. , Graham C. 2009 . Identification and dynamics of a cryptic suture zone in tropical rainforest . Proc. R. Soc. Lond. B Biol. Sci . rspb. 2008.1622 . Nater A. , Mattle-Greminger M.P. , Nurcahyo A. , Nowak M.G. , de Manuel M. , Desai T. , Groves C. , Pybus M. , Sonay T.B. , Roos C. 2017 . Morphometric, behavioral, and genomic evidence for a new Orangutan species . Curr. Biol. 27 : 3487 – 3498.e3410 . Google Scholar CrossRef Search ADS PubMed Nei M. , Li W.-H. 1979 . Mathematical model for studying genetic variation in terms of restriction endonucleases . Proc. Natl. Acad. Sci. U.S.A 76 : 5269 – 5273 . Google Scholar CrossRef Search ADS PubMed Niemiller M.L. , Near T.J. , Fitzpatrick B.M. 2012 . Delimiting species using multilocus data: diagnosing cryptic diversity in the southern cavefish, Typhlichthys subterraneus (Teleostei: Amblyopsidae) . Evolution 66 : 846 – 866 . Google Scholar CrossRef Search ADS PubMed Nosil P. 2008 . Speciation with gene flow could be common . Mol. Ecol. 17 : 2103 – 2106 . Google Scholar CrossRef Search ADS PubMed Nosil P. , Flaxman S.M. 2011 . Conditions for mutation-order speciation . Proc. R. Soc. Lond. B Biol. Sci. 278 : 399 – 407 . Google Scholar CrossRef Search ADS Ogilvie H.A. , Bouckaert R.R. , Drummond A.J. 2017 . StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates . Mol. Biol. Evol. 34 : 2101 – 2114 . Google Scholar CrossRef Search ADS PubMed Olave M. , Solà E. , Knowles L.L. 2014 . Upstream analyses create problems with DNA-based species delimitation . Syst. Biol. 63 : 263 – 271 . Google Scholar CrossRef Search ADS PubMed Oliver P. , Keogh J.S. , Moritz C. 2015 . New approaches to cataloguing and understanding evolutionary diversity: a perspective from Australian herpetology . Aust. J. Zool. 62 : 417 – 430 . Oliver P.M. , Adams M. , Doughty P. 2010 . Molecular evidence for ten species and Oligo-Miocene vicariance within a nominal Australian gecko species (Crenadactylus ocellatus, Diplodactylidae) . BMC Evol. Biol. 10 : 386 . Google Scholar CrossRef Search ADS PubMed Olsson U. , Alström P. , Ericson P.G. , Sundberg P. 2005 . Non-monophyletic taxa and cryptic species—evidence from a molecular phylogeny of leaf-warblers (Phylloscopus, Aves) . Mol. Phylogenet. Evol. 36 : 261 – 276 . Google Scholar CrossRef Search ADS PubMed Padial J.M. , Miralles A. , De la Riva I. , Vences M. 2010 . The integrative future of taxonomy . Front. Zool. 7 : 16 . Google Scholar CrossRef Search ADS PubMed Pante E. , Puillandre N. , Viricel A. , Arnaud Haond S. , Aurelle D. , Castelin M. , Chenuil A. , Destombe C. , Forcioli D. , Valero M. 2015 . Species are hypotheses: avoid connectivity assessments based on pillars of sand . Mol. Ecol. 24 : 525 – 544 . Google Scholar CrossRef Search ADS PubMed Pante E. , Schoelinck C. , Puillandre N. 2014 . From integrative taxonomy to species description: one step beyond . Syst. Biol. 64 : 152 – 160 . Google Scholar CrossRef Search ADS PubMed Pease J.B. , Hahn M.W. 2015 . Detection and polarization of introgression in a five-taxon phylogeny . Syst. Biol. 64 : 651 – 662 . Google Scholar CrossRef Search ADS PubMed Perkins S.L. 2000 . Species concepts and malaria parasites: detecting a cryptic species of Plasmodium . Proc. R. Soc. Lond. B Biol. Sci. 267 : 2345 – 2350 . Google Scholar CrossRef Search ADS Phillips B.L. , Baird S.J. , Moritz C. 2004 . When vicars meet: a narrow contact zone between morphologically cryptic phylogeographic lineages of the rainforest skink, Carlia rubrigularis . Evolution 58 : 1536 – 1548 . Google Scholar CrossRef Search ADS PubMed Pinho C. , Hey J. 2010 . Divergence with gene flow: models and data . Ann. Rev. Ecol. Evol. Syst. 41 : 215 – 230 . Google Scholar CrossRef Search ADS Potter S. , Bragg J.G. , Peter B.M. , Bi K. , Moritz C. 2016 . Phylogenomics at the tips: inferring lineages and their demographic history in a tropical lizard, Carlia amax . Mol. Ecol. 25 : 1367 – 1380 . Google Scholar CrossRef Search ADS PubMed Rabosky D.L. , Matute D.R. 2013 . Macroevolutionary speciation rates are decoupled from the evolution of intrinsic reproductive isolation in Drosophila and birds . Proc. Natl. Acad. Sci. U.S.A 110 : 15354 – 15359 . Google Scholar CrossRef Search ADS PubMed Rannala B. 2015 . The art and science of species delimitation . Curr. Zool. 61 : 846 – 853 . Google Scholar CrossRef Search ADS Richardson B.J. , Baverstock P.R. , Adams M. 1986 . Allozyme electrophoresis: a handbook for animal systematics and population studies . San Diego (CA) : Academic Press . Rissler L.J. , Apodaca J.J. 2007 . Adding more ecology into species delimitation: ecological niche models and phylogeography help define cryptic species in the black salamander (Aneides flavipunctatus) . Syst. Biol. 56 : 924 – 942 . Google Scholar CrossRef Search ADS PubMed Ronquist F. , Teslenko M. , Van Der Mark P. , Ayres D.L. , Darling A. , Höhna S. , Larget B. , Liu L. , Suchard M.A. , Huelsenbeck J.P. 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Syst. Biol. 61 : 539 – 542 . Google Scholar CrossRef Search ADS PubMed Rosenblum E.B. , Sarver B.A. , Brown J.W. , Des Roches S. , Hardwick K.M. , Hether T.D. , Eastman J.M. , Pennell M.W. , Harmon L.J. 2012 . Goldilocks meets Santa Rosalia: an ephemeral speciation model explains patterns of diversification across time scales . J. Evol. Biol. 39 : 255 – 261 . Google Scholar CrossRef Search ADS Rosindell J. , Cornell S.J. , Hubbell S.P. , Etienne R.S. 2010 . Protracted speciation revitalizes the neutral theory of biodiversity . Ecol. Lett. 13 : 716 – 727 . Google Scholar CrossRef Search ADS PubMed Roux C. , Fraisse C. , Romiguier J. , Anciaux Y. , Galtier N. , Bierne N. 2016 . Shedding light on the grey zone of speciation along a continuum of genomic divergence . PLoS Biol. 14 : e2000234 . Google Scholar CrossRef Search ADS PubMed Sasa M.M. , Chippindale P.T. , Johnson N.A. 1998 . Patterns of postzygotic isolation in frogs . Evolution 52 : 1811 – 1820 . Google Scholar CrossRef Search ADS PubMed Schluter D. 2001 . Ecology and the origin of species . Trends Ecol. Evol. 16 : 372 – 380 . Google Scholar CrossRef Search ADS PubMed Schneider C. , Moritz C. 1999 . Rainforest refugia and Australia”s Wet Tropics . Proc. R. Soc. Lond. B Biol. Sci. 266 : 191 – 196 . Google Scholar CrossRef Search ADS Schonrogge K. , Barr B. , Wardlaw J.C. , Napper E. , Gardner M.G. , Breen J. , Elmes G.W. , Thomas J.A. 2002 . When rare species become endangered: cryptic speciation in myrmecophilous hoverflies . Biol. J. Linnean Soc. 75 : 291 – 300 . Google Scholar CrossRef Search ADS Seehausen O. , Takimoto G. , Roy D. , Jokela J. 2008 . Speciation reversal and biodiversity dynamics with hybridization in changing environments . Mol. Ecol. 17 : 30 – 44 . Google Scholar CrossRef Search ADS PubMed Silva A.C.A. , Santos N. , Ogilvie H.A. , Moritz C. 2017 . Validation and description of two new north-western Australian Rainbow skinks with multispecies coalescent methods and morphology . PeerJ 5 : e3724 . Google Scholar CrossRef Search ADS PubMed Singhal S. 2013 . De novo transcriptomic analyses for non model organisms: an evaluation of methods across a multi species data set . Mol. Ecol. Resour. 13 : 403 – 416 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Bi K. 2017 . History cleans up messes: the impact of time in driving divergence and introgression in a tropical suture zone . Evolution 71 : 1888 – 1899 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Moritz C. 2012 . Strong selection against hybrids maintains a narrow contact zone between morphologically cryptic lineages in a rainforest lizard . Evolution 66 : 1474 – 1489 . Google Scholar CrossRef Search ADS PubMed Singhal S. , Moritz C. 2013 . Reproductive isolation between phylogeographic lineages scales with divergence . Proc. R. Soc. Lond. B Biol. Sci. 280 : 20132246 . Google Scholar CrossRef Search ADS Sites J.W. , Barton N.H. , Reed K.M. 1995 . The genetic structure of a hybrid zone between two chromosome races of the Sceloporus grammicus complex (Sauria, Phrynosomatidae) in Central Mexico . Evolution 49 : 9 – 36 . Google Scholar CrossRef Search ADS PubMed Skinner A. , Hugall A.F. , Hutchinson M.N. 2011 . Lygosomine phylogeny and the origins of Australian scincid lizards . J. Biogeogr. 38 : 1044 – 1058 . Google Scholar CrossRef Search ADS Slater G.S.C. , Birney E. 2005 . Automated generation of heuristics for biological sequence comparison . BMC Bioinformatics 6 : 31 . Google Scholar CrossRef Search ADS PubMed Smith M.A. , Rodriguez J.J. , Whitfield J.B. , Deans A.R. , Janzen D.H. , Hallwachs W. , Hebert P.D. 2008 . Extreme diversity of tropical parasitoid wasps exposed by iterative integration of natural history, DNA barcoding, morphology, and collections . Proc. Natl. Acad. Sci. U.S.A 105 : 12359 – 12364 . Google Scholar CrossRef Search ADS PubMed Smith M.A. , Woodley N.E. , Janzen D.H. , Hallwachs W. , Hebert P.D. 2006 . DNA barcodes reveal cryptic host-specificity within the presumed polyphagous members of a genus of parasitoid flies (Diptera: Tachinidae) . Proc. Natl. Acad. Sci. U.S.A 103 : 3657 – 3662 . Google Scholar CrossRef Search ADS PubMed Stewart K.A. , Lougheed S.C. 2013 . Testing for intraspecific postzygotic isolation between cryptic lineages of Pseudacris crucifer . Ecol. Evol. 3 : 4621 – 4630 . Google Scholar CrossRef Search ADS PubMed Struck T.H. , Feder J.L. , Bendiksby M. , Birkeland S. , Cerca J. , Gusarov V.I. , Kistenich S. , Larsson K.-H. , Liow L.H. , Nowak M.D. 2018 . Finding evolutionary processes hidden in cryptic species . Trends Ecol. Evol. 33 : 153 – 163 . Google Scholar CrossRef Search ADS PubMed Stuart B.L. , Inger R.F. , Voris H.K. 2006 . High level of cryptic species diversity revealed by sympatric lineages of Southeast Asian forest frogs . Biol. Lett. 2 : 470 – 474 . Google Scholar CrossRef Search ADS PubMed Suatoni E. , Vicario S. , Rice S. , Snell T. , Caccone A. 2006 . An analysis of species boundaries and biogeographic patterns in a cryptic species complex: the rotifer—Brachionus plicatilis . Mol. Phylogenet. Evol. 41 : 86 – 98 . Google Scholar CrossRef Search ADS PubMed Sukumaran J. , Knowles L.L. 2017 . Multispecies coalescent delimits structure, not species . Proc. Natl. Acad. Sci. U.S.A 114 : 1607 – 1612 . Google Scholar CrossRef Search ADS PubMed Sunnucks P. , Hales D.F. 1996 . Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae) . Mol. Biol. Evol. 13 : 510 – 524 . Google Scholar CrossRef Search ADS PubMed Tellier F. , Vega J.A. , Broitman B.R. , Vasquez J.A. , Valero M. , Faugeron S. 2011 . The importance of having two species instead of one in kelp management: the Lessonia nigrescens species complex . Cah. Biol. Mar. 52 : 455 – 465 . Turner J.R. 1967 . Why does the genotype not congeal? Evolution 21 : 645 – 656 . Google Scholar CrossRef Search ADS PubMed Vilas R. , Criscione C. , Blouin M. 2005 . A comparison between mitochondrial DNA and the ribosomal internal transcribed regions in prospecting for cryptic species of platyhelminth parasites . Parasitology 131 : 839 – 846 . Google Scholar CrossRef Search ADS PubMed Winger B.M. , Bates J.M. 2015 . The tempo of trait divergence in geographic isolation: avian speciation across the Marañon Valley of Peru . Evolution 69 : 772 – 787 . Google Scholar CrossRef Search ADS PubMed Witt J.D. , Threloff D.L. , Hebert P.D. 2006 . DNA barcoding reveals extraordinary cryptic diversity in an amphipod genus: implications for desert spring conservation . Mol. Ecol. 15 : 3073 – 3082 . Google Scholar CrossRef Search ADS PubMed Yang Z. , Rannala B. 2014 . Unguided species delimitation using DNA sequence data from multiple loci . Mol. Biol. Evol. 31 : 3125 – 3135 . Google Scholar CrossRef Search ADS PubMed Zhang C. , Zhang D.-X. , Zhu T. , Yang Z. 2011 . Evaluation of a Bayesian coalescent method of species delimitation . Syst. Biol. 60 : 747 – 761 . Google Scholar CrossRef Search ADS PubMed Zhang J. , Kobert K. , Flouri T. , Stamatakis A. 2013 . PEAR: a fast and accurate Illumina Paired-End reAd mergeR . Bioinformatics 30 : 614 – 620 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

### Journal

Systematic BiologyOxford University Press

Published: Apr 4, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations