Population dynamics within species at the edge of their distributional range, including the formation of genetic structure during range expansion, are difficult to study when they have had limited time to evolve. Western Fence Lizards (Sceloporus occidentalis) have a patchy distribution at the northern edge of their range around the Puget Sound, Washington, where they almost exclusively occur on imperiled coastal habitats. The entire region was covered by Pleistocene glaciation as recently as 16,000 years ago, suggesting that populations must have colonized these habitats relatively recently. We tested for population differentiation across this landscape using genome-wide SNPs and morphological data. A time-calibrated species tree supports the hypothesis of a post-glacial establishment and subsequent population expansion into the region. Despite a strong signal for fine-scale population genetic structure across the Puget Sound with as many as 8–10 distinct subpopulations supported by the SNP data, there is minimal evidence for morphological differentiation at this same spatiotemporal scale. Historical demographic analyses suggest that populations expanded and diverged across the region as the Cordilleran Ice Sheet receded. Population isolation, lack of dispersal corridors, and strict habitat requirements are the key drivers of population divergence in this system. These same factors may prove detrimental to the future persistence of populations as they cope with increasing shoreline development associated with urbanization.

We extracted genomic DNA from liver biopsies using salt-extraction and then conducted double digest restriction-site associated DNA sequencing (ddRADseq). We double-digested each sample using the digestion enzymes SbfI and MspI in CutSmart Buffer (New England Biolabs) for 7 hours at 37 C. For fragment purification, we used Sera-Mag SpeedBeads. We then prepared a master mix for eight distinct barcodes to be ligated to the cut sites of the fragmented DNA. The libraries were size-selected (between 415 and 515 bp after accounting for adapter length) on a Blue Pippin Prep size fractionator (Sage Science). For the final library amplification, we used Phusion Hi-Fidelity DNA Polymerase and Illumina’s indexed primers. We determined the concentration and size distribution of each indexed pool using an Agilent 2200 TapeStation. Lastly, we sent the quantified pools to QB3-Berkeley Genomics, UC Berkeley for qPCR to determine sequenceable library concentrations before multiplexing equimolar amounts of each pool for sequencing on one Illumina HiSeq 4000 lane (51-bp, single-end reads; 11 pools containing up to 8 samples each). The demultiplexed sequences are deposited at the Sequence Read Archive (SRA; BioProject ID: PRJNA757434; Table S1).

We demultiplexed each sample from their respective pool using their unique barcode sequence using iPyRAD v.0.9.50. We conducted a reference-based assembly of the RAD loci using a draft of the S. occidentalis genome from Yosemite National Park, California. A sequence similarity threshold of 90% was used to cluster reads within samples and loci between samples. We removed consensus sequences with low coverage (< 6 reads), excessive undetermined or heterozygous sites (> 5%), too many alleles for a sample (> 2 for diploids), and an excess of shared heterozygosity among samples (paralog filter = 0.5). For the final alignments, we generated output files containing 0% missing data (1,037 loci) and 50% missing data (3,491 loci; Table S2).