Species are the fundamental unit of biodiversity studies. However, many species complexes are difficult to delimit, especially those characterized by complicated patterns of population structure. Salamanders in the family Plethodontidae often form species by slowly fragmenting across a landscape over space and time. They thus provide many examples of species complexes in which gradual Darwinian evolution has resulted in multiple units of varying degrees of differentiation, including incompletely separated lineages. Here we report on a molecular systematic investigation of woodland salamanders in the Plethodon wehrlei group, which has recently been split from two species into five. To quantify patterns of genetic variation, we collected genetic samples from 24 individuals from 20 populations, including all species and representing a carefully selected subset of previous work. From these samples, we obtained genomic data using anchored hybrid enrichment, which resulted in 319 loci averaging 1300 base pairs in length. Biallelic single nucleotide polymorphisms (SNPs) were randomly selected from 316 of these loci for some analyses. We examined patterns of genetic structure using PCA, DAPC, FEEMS, and STRUCTURE, and found that all of the recognized species formed genetic clusters; however, P. wehrlei and P. punctatus were relatively weakly differentiated, and STRUCTURE identified three separate clusters within P. jacksoni. Species trees inferred using wASTRAL, BPP, and TreeMix all recovered the same topology, with P. dixi sister to the other taxa, which included a northern clade (P. wehrlei, P. punctatus, P. pauleyi) and a southern clade (P. jacksoni, with three separate groups). TreeMix only inferred one gene flow event. We evaluated the candidate species using BPP and the genealogical divergence index (gdi). While BPP delimited all candidate species with strong support (all posterior probabilities = 1.0), the gdi only strongly supported P. dixi and P. pauleyi, both of which have only been recently described. We discuss the difficult problem of species delimitation in groups that form species via range fragmentation. We also provide a vision for future research with the aim of better testing and diagnosing the species diversity within the P. wehrlei group.

Blood samples and tail tips were collected from 24 individuals from 20 populations of the P. wehrlei complex. Genomic DNA was extracted using Qiagen DNeasy Blood and Tissue Kits (Qiagen Corp., Valencia, CA). All samples were processed at the Florida State University Center for Anchored Phylogenomics using a hybrid enrichment protocol (Lemmon et al. 2012). Briefly, each sample was sonicated to an average fragment size of 150–300 bp using a Covaris E220 focused ultrasonicator. Library preparation and indexing were done on a Beckman-Coulter Biomex FXp liquid-handling robot according to Meyer & Kircher (2010). Pooled libraries were enriched with an Agilent Custom SureSelect kit that contained amphibian-specific probes (Heinicke et al. 2018). After enrichment, samples were pooled for sequencing on one PE150 lane of an Illumina HiSeq 2000. Raw sequencing reads were filtered for quality, de-multiplexed, and assembled following Lemmon et al. (2012). Sequences were aligned using MAFFT (Katoh and Standley 2013), and the loci were phased using read overlap information in a Bayesian approach (Pyron et al. 2016; Booker et al. 2022).