North‐facing slopes and elevation shape asymmetric genetic structure in the range‐restricted salamander Plethodon shenandoah

Abstract Species with narrow environmental tolerances are often distributed within fragmented patches of suitable habitat, and dispersal among these subpopulations can be difficult to directly observe. Genetic data can help quantify gene flow between localities, which is especially important for vulnerable species with a disjunct range. The Shenandoah salamander (Plethodon shenandoah) is a federally endangered species known only from three mountaintops in Virginia, USA. To reconstruct the evolutionary history and population connectivity of this species, we generated both mitochondrial and nuclear data using sequence capture from individuals collected across all three mountaintops. Applying population and landscape genetic methods, we found strong population structure that was independent of geographic distance. Both the nuclear markers and mitochondrial genomes indicated a deep split between the most southern population and the genetically similar central and northern populations. Although there was some mitochondrial haplotype‐splitting between the central and northern populations, there was admixture in nuclear markers. This is indicative of either a recent split or current male‐biased dispersal among mountain isolates. Models of landscape resistance found that dispersal across north‐facing slopes at mid‐elevation levels best explain the observed genetic structure among populations. These unexpected results highlight the importance of incorporating landscape features in understanding and predicting the movement and fragmentation of this range‐restricted salamander species across space.

2 Figure S3: Resistance surfaces for the same landscape variables generated for Figure S2, but restricted for elevation by only allowing cells above 500m. This was done to create more biologically relevant dispersal paths for P. shenandoah. Locality centers indicated by white dots and named in the first panel.

Landscape resistance scores:
All resistance surfaces were scaled from 0-100 and in order to choose the direction of the scale we used a combination of biological knowledge of the species and locality information. We calculated the mean values for all GPS-locations and the 8 surrounding cells and compared that to the mean of all values on the elevation-restricted raster. Unless noted below we chose the direction that would cause the GPS-localities to have a lower resistance value, assuming that individuals would prefer to disperse over similar habitat as they currently inhabit. For some, surfaces values were close and we used biological knowledge of the species to make a final call (see Table S1).  We chose five samples each of Plethodon shenandoah and Plethodon cinereus for RADtag sequencing to develop probes that could both differentiate between species as well as between populations. This included two individuals from Mountain N, and four each from Mountain C and S. Six samples were digested with restriction enzyme SbfI and following the ligation of the first adapter they were sheared to around 300 bp for sequencing. Four samples were digested separately with BspqI and sheared to the same size. Libraries were quantified using qPCR and pooled at equimolar concentrations for sequencing.
Libraries were sequenced on two MiSeq runs at the Dana Farber Cancer institute using paired end 150 bp reads for all 10 samples. We mapped the reads onto two bacterial genomes (Helicobacter pylori and Delftia acidovorans) to remove potential bacterial contaminant DNA. Resulting sequences were run through stacks 1.05 using default settings (Catchen et al., 2013). As the majority of sequences were not stacking properly we ran assemblies per individual using Geneious 7.0.1 (Kearse et al., 2012).
The contigs showed that sequences did align, but not all starting at the same bp as would be expected following a successful restriction enzyme digest. For this reason, we used several different approaches to develop 100 bp capture bait sequences to increase our chances of obtaining valid SNPs.
We developed 12,614 baits using the discoSNP software that finds SNPs during the de Bruijn stage of assembly (Uricaru et al., 2015). A further 1575 baits were developed from the Geneious assemblies by manually finding single SNPs in high quality contigs. All assembled contigs were also blasted to GenBank using Blast+ to identify sequences that had known amphibian annotations and all NONmitochondrial hits (660 total) were also used. These baits were filtered for overlapping sequences and bad hybridization factors to end up with a total of 10,000 baits that were ordered using the MyBaits set from MycroArray (Ann Harbor, Michigan, USA) These probes were tested on 63 samples on one HiSeq run to check for enrichment. Following preliminary analyses, we restricted our targets to the 875 loci that showed medium enrichment in order to target only sequences that enriched well but were less likely to be repetitive elements. These 875 loci were supplemented with 334 UCEs, 43 conserved genes and 15kb of the mitochondrial genome (see main text) and were enriched with a new MyBaits probe-set for all samples.