Single nucleotide polymorphism (SNPs) data for Scurria scurra, Scurria variabilis, Scurria ceciliana and Scurria araucana

The distribution of genetic diversity is often heterogeneous in space, and it usually correlates with environmental transitions or historical processes that affect demography. The coast of Chile encompasses two biogeographic provinces and spans a broad environmental gradient, together with oceanographic processes linked to coastal topography that can affect species' genetic diversity. Here, we evaluated the genetic connectivity and historical demography of four Scurria limpets, S. scurra, S. variabilis, S. ceciliana, and S. araucana, between ca. 19° S and 53° S in the Chilean coast using genome-wide SNP markers. Genetic structure varied among species, which was evidenced by species-specific breaks together with two shared breaks. One of the shared breaks was located at 22–25° S and was observed in S. araucana and S. variabilis, while the second break around 31–34° S was shared by three Scurria species. Interestingly, the identified genetic breaks are also shared with other low-dispersing invertebrates. Demographic histories show bottlenecks in S. scurra and S. araucana populations and recent population expansion in all species. The shared genetic breaks can be linked to oceanographic features acting as soft barriers to dispersal and also to historical climate, evidencing the utility of comparing multiple and sympatric species to understand the influence of a particular seascape on genetic diversity.

Methods
These data were retrieved by restriction-site associated DNA sequencing of muscle (foot) tissue using the enzyme pstI. 

Raw reads were demultiplexed using Stacks. Scurria scurra, Scurria cecilian, and Scurria araucana RAD loci were mapped to the reference genome of Scurria scurria (Saenz-Agudelo P., unpublished) using the ref_map.pl pipeline of Stacks, while for Scurria variabilis loci were called de novo using the denovo_map. Pipeline of Stacks.

VCF files were generated in populations from Stacks, keeping only 80% as the minimum percentage of all individuals to process a locus. Further filtering steps were made in VCFtools, namely mean minimum read depth per locus (15), maximum mean depth per site (62-73, depending on the species, es), and a final genotype call rate of 90%. No individuals with more than 20% missing data were kept. The final filter consisted of removing loci with evidence of linkage disequilibrium, which was estimated using the function snpgdsLDpruning from the R package SNPRelate.

For detecting putatively outlier loci, BayeScan and pcadapt were used, and the common loci detected by both analyses were removed from the "NoOuts" datasets, and only these outlier loci were kept in the "Outs" datasets.