Hybrid zone movement
For each transect-year we fit hybrid zone equilibrium cline models using population ancestry proportions to estimate the shape and center of the hybrid zone. First, for each transect, the distance (in km) of each trap from the respective westernmost trap for that transect (trap T00 [42.6714° N, 73.0145° W] for the Massachusetts transect, and trap CT01 [41.0798° N, 73.7054° W] for the Connecticut transect), was calculated using the Latitude/Longitude Distance Calculator available at (https://www.nhc.noaa.gov/gccalc.shtml). Then, for each transect-year, the population coefficient of assignment (Q ) to Bruce spanworm at each trap was estimated using the software program Structure v.2.3.2 (Falush, Stephens, & Pritchard 2003; Pritchard, Stephens, & Donnelly, 2000) based on two population (K= 2) analyses using the admixture model, correlated allele frequencies, and default settings, with random starting values, runtimes of 200,000 generations, and burn-in periods of 20,000 generations, for each year-transect combination. Finally, using the R package ‘hzar’ v 0.2-5 (Derryberry, Derryberry, Maley & Brumfield, 2014), a null model, plus three different hybrid-zone models: 1) minimum and maximum frequencies fixed to 0 and 1, and with no exponential decay tails, 2) minimum and maximum frequencies as free parameters, with no exponential decay tails, and 3) minimum and maximum frequencies as free parameters, with both tails as independent parameters, were used to estimate the center of the hybrid zone based on three independent analyses with chain lengths of 100,000 generations, burn-in periods of 10,000 generations, using random starting variables for each analysis. The maximum likelihood values for each run were then compared to determine which model provided the best fit to the observed dataset, from which the center of the hybrid zone and the hybrid cline was estimated for each transect-year combination.