DOI: 10.5061/dryad.xksn02vq1
title: Hard edges, soft edges, and species range evolution: A genomic analysis of the Cumberland Plateau salamander

description: 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.
description: We gathered blood samples and tail tips from 39 individuals in 31
 localities of <em>Plethodon kentucki</em>. <em>Plethodon
 glutinosus</em>, a morphologically cryptic relative, occurs
 sympatrically, so a known <em>P. glutinosus</em> 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
 <em>EcoRI </em>and <em>SphI </em>(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 <em>process_radtags</em> 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 <em>ustacks</em> to build loci
 within individuals, <em>cstacks</em> to assemble a catalog of
 loci across individuals, <em>sstacks</em> to match samples to
 catalog, <em>tsv2bam</em> to convert files,
 <em>gstacks</em> to genotype individuals, and
 <em>populations</em> to compute summary statistics and export
 files. We then further excluded two individuals with &gt; 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.
description: # Hard edges, soft edges, and species range evolution: A genomic analysis
 of the Cumberland Plateau salamander
 [https://doi.org/10.5061/dryad.xksn02vq1](https://doi.org/10.5061/dryad.xksn02vq1) This repository contains the demultiplexed sequence data, SNPs, and haplotypes from 39 individuals in 31 localities of *Plethodon kentucki* and *Plethodon glutinosus*. The sequencing method we used here was double-digest restriction-site associated (ddRAD) sequencing. The samples were collected as part of our study on range limit evolution in *P. kentucki*, and as the two species are almost morphologically indistinguishable, we collected specimens from both species. It also contains the Supplementary Information associated with this publication and the Overview, Data, Model, Assessment and Prediction (ODMAP) Protocol used to create the published Environmental Niche Model (ENM). Scripts used to process these data can be found on the [companion GitHub](https://github.com/efwatts/Plethodonkentucki_RangeLimits) repository for this manuscript.  ## Data Availability ### Sequence data "P_kentucki_ddRAD_0missingdata.fa" - FASTA file N = 39, 0% missing data per SNP, 6,803 SNPs ### SNP data "P_kentucki_ddRAD_0missingdata.structure" - SNPs N = 39, 0% missing data per SNP, 6,803 SNPs "P_kentucki_ddRAD_10missingdata.structure" - SNPs N = 39, 10% missing data per SNP, 30,155 SNPs ### Haplotype data "P_kentucki_ddRAD_0missingdata_haplotypes.tsv" - haplotypes N = 39, 0% missing data per SNP, 6,803 SNPs, (adenine (A), guanine (G), cytosine (C), and thymine (T)) "P_kentucki_ddRAD_10missingdata_haplotypes.tsv" - haplotypes N = 39, 10% missing data per SNP, 30,155 SNPs, (adenine (A), guanine (G), cytosine (C), and thymine (T)) ### Usage notes FASTA (.fa) files can be viewed in a standard text editor  .structure files can be viewed in a standard text editor .tsv files can be viewed in Excel,imported to Google Sheets, or in a standard text editor ### Sharing/Access information The sequence data deposited here are also available in fastq format on the Sequence Read Archive (SRA): * [https://www.ncbi.nlm.nih.gov/sra/PRJNA1109840](https://www.ncbi.nlm.nih.gov/sra/PRJNA1109840) ## Supplementary Information (uploaded to Zenodo) ### Supplementary Tables "Table.S1.xlsx" - WorldClim variables used in the construction of the environmental niche model (ENM) for *Plethodon kentucki*. "Table.S2.xlsx" - Sample numbers and locality information for all *P. kentucki* localities. "Table.S3.xlsx" - Pairwise estimates of FST below the diagonal and pairwise geographic Euclidean distances above the diagonal. Population numbers correspond with Supplemental Table 2. "Table.S4.xlsx" - Genetic connectivity information between localities of *P. kentucki*. "Table.S5.xlsx" - Mutation-scaled proportion of migrants per generation (*M*) determined by MIGRATE-N between pairs of populations of *P. kentucki* populations in the eastern group. The left population is the source population, and the right population is the recipient population. Arrows between populations them indicate direction of migration. "Table.S6.xlsx" - Mutation-scaled proportion of migrants per generation (*M*) determined by MIGRATE-N between pairs of populations of *P. kentucki* populations in the western group. The left population is the source population, and the right population is the recipient population. Arrows between populations them indicate direction of migration. "Table.S7.xlsx" - Effective population size (θNe) estimates of *P. kentucki* populations in the eastern group. "Table.S8.xlsx" - Effective population size (θNe) estimates of *P. kentucki* populations in the western group. ### Supplementary Figures "Fig.S1.png" - Omission rate and predicted area as a function of the cumulative threshold. "Fig.S2.png" - 'tess3r' and ‘sNMF’ (LEA)  cross-entropy scores for 1-20 ancestral populations used to distinguish *P. kentucki* and *P. glutinosus* samples. "Fig.S3.png" - 'tess3r' and 'sNMF' (LEA) bar plots used to distinguish between *P. kentucki* and *P. glutinosus* samples using K = 2. Population numbers correspond to population numbers used in Figure 1. The individual labeled “G” is the known *P. glutinosus* individual included to help with identification. "Fig.S4.png" - ‘tess3r’ and ‘sNMF’ (LEA) cross-entropy scores for 1-20 ancestral populations "Fig.S5.png" - Discriminant Analysis of Principal Components (DAPC) of the *P. kentucki* individuals used in this study. This DAPC confirmed group membership to the sNMF groups. "Fig.S6.png" - Map of *P. kentucki* localities with FST values between pairs of populations. Lines between localities have FST values labeled. ### Overview, Data, Model, Assessment and Prediction (ODMAP) Protocol "ODMAP_Watts_EtAl_2024-05-06.csv" - WorldClim variables used along with settings used for parameterization and creation of ENM.