Aim: Gene flow from central to edge populations is thought to limit population growth at range edges by constraining local adaptation. In this study, we explore the thesis that range edges can differ in their dynamics and be either “hard” (e.g. a river) or “soft” (e.g. ecological gradients). We hypothesize that soft edge populations will have smaller effective population sizes than central populations and that gene flow will be greater from the center to the edge than vice versa. Conversely, we hypothesize that hard edge populations should have similar effective population sizes to central populations and that gene flow will be equal between the two.

Location: Kentucky, West Virginia, and Virginia, USA. Taxon: Plethodon kentucki (Caudata: Plethodontidae).

Methods: We evaluated landscape suitability using an ecological niche model, then we compared gene flow and effective population sizes between edge and central populations and quantified gene flow between populations. Finally, we characterized landscape genetic variation, testing for isolation by distance and isolation by environment. Results: We found continuously decreasing habitat quality along soft edges, with hard edges more variable. Additionally, we found that soft edges had lower effective population sizes than central populations and that gene flow was greater from the center of the range to the soft edges than the reverse. In hard edges, by contrast, we found effective population sizes in edge populations were similar to central populations, with relatively equal gene flow in both directions.

Main conclusions: Understanding why species have range limits is central to investigations of the structure of biodiversity, yet the evolutionary dynamics of range edges remain poorly understood. We show that within a single species with a small range, the evolutionary dynamics operating at range boundaries may depend on the nature of the boundary.

We gathered blood samples and tail tips from 39 individuals in 31 localities of Plethodon kentucki. Plethodon glutinosus, a morphologically cryptic relative, occurs sympatrically, so a known P. glutinosus sample was included to identify species. Genomic DNA was extracted using Qiagen DNeasy Blood and Tissue Kits (Qiagen Corp., Valencia CA) following the manufacturer’s protocol.

We quantified genetic variation using single nucleotide polymorphisms (SNPs). To obtain SNPs, we used double-digest restriction-site associated DNA sequencing (ddRAD) (Peterson et al., 2012), including a protocol that has been optimized for use in salamanders (Jones and Weisrock, 2018). Briefly, we double-digested extracted DNA using equal amounts of the restriction enzymes EcoRI and SphI (New England BioLabs). DNA fragments were indexed for each individual and pooled for size selection of 420 bp +/- 10% on a Pippin Prep (Sage Science, Beverly MA). The resulting libraries were amplified by PCR for 12 cycles with Phusion high-fidelity DNA polymerase (New England Biolabs, Ipswich MA) and cleaned with Dynabeads (ThermoFisher, Waltham MA) and AMPure XP beads (Beckman Coulter, Inc., Brea CA). They were sequenced with 150 bp paired-end reads using a 1% PhiX DNA spike-in on an Illumina (Illumina, San Diego CA) HiSeq 2500 at Novogene.

To filter the raw sequencing data, we first checked the quality of reads with fastQC v. 0.11.9 (Andrews, 2010). Then, we demultiplexed raw sequences in Stacks v. 2.61 (Catchen et al., 2011; Catchen et al. 2013) following Rochette and Catchen (2017). We built a custom pipeline based on Hime et al. (2019) (Appendix A1-10). We demultiplexed the raw, stitched reads by individual with the process_radtags function, allowing for one mismatch between observed and expected barcodes. We retained reads with both restriction enzyme cut sites and had a mean Phred quality score greater than 20 over 45 bp sliding-window intervals (Hime et al., 2019).

We excluded two individuals with fewer than 900,000 reads (DH_64627 and SRK_2997). We optimized parameters by testing M 1-12 using the R80 method following Rochette and Catchen (2017) and found M12 to retain the greatest number of reads. We then assembled the sequences de novo. We used ustacks to build loci within individuals, cstacks to assemble a catalog of loci across individuals, sstacks to match samples to catalog, tsv2bam to convert files, gstacks to genotype individuals, and populations to compute summary statistics and export files. We then further excluded two individuals with > 95% missing data (EFW_0006 and SRK_3200), resulting in a final sample of 35 individuals from 31 localities. When allowing 10% missing data per SNP, we obtained 30,155 SNPs, and when we did not allow any missing data, we retained 6,803 SNPs.

For other methods, see the Materials and Methods section of Watts et al. 2024.