Hard edges, soft edges, and species range evolution: A genomic analysis of the Cumberland Plateau salamander
Data files
May 20, 2024 version files 193.30 MB
Abstract
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.
README: Hard edges, soft edges, and species range evolution: A genomic analysis of the Cumberland Plateau salamander
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 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):
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.
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.