Hard edges, soft edges, and species range evolution: A genomic analysis of the Cumberland Plateau salamander

Description

Notes

Methods

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.