Materials and methods

Sampling sites

Puget Sound (Fig. 1) is part of the Salish Sea, an estuarine system in Washington State, USA, linked to Pacific Ocean via three connections to the Strait Juan de Fuca. Four marinas were studied, three (Edmonds, Shilshole and Des Moines marinas) located in the central basin, and Shelton Yacht Club, located in the southern basin. With 36,000 recreational boats, approximately 24,000 actively cruising, the Salish Sea is a famous destination for recreational boating (City of Des Moines Marina, http://www.desmoinesmarina.com/uploads/7/2/2/4/72248139/city_of
has ~1500 mooring slips, and is the biggest marina considered in this study. The City of Des Moines Marina has anywhere between 400 - 800 boats moored, with a constant flow of 10-30 daily visitors (Tara Reilly, Office Specialist of City of Des Moines Marina, pers. comm.). Edmonds Marina has 668 mooring slips and receives around 3100 visiting boats per year. The remote Shelton Yacht Club offers just 109 mooring slips, and is surrounded by extensive shellfish aquaculture.

Sampling procedures, DNA extraction and microsatellite typing

Botryllus schlosseri samples were collected every two years during the period 1999 to 2007, and in 2013 and 2018 (7 collection sessions; Suppl. Table A.1). Colonies, or colony fragments separated by at least one meter, were removed using single edge-razor blades and put into separate 1.5ml vials. The DNA extraction followed the protocol in Paz, Douek, Mo, Goren and Rinkevich (2003).
Eight B. schlosseri microsatellites, BS-2, BS-8, BS-9 (Abdoullaye et al., 2010), BS-811 (Pancer, Gershon, & Rinkevich, 1994), PB-29, PB-41, PB-49 and PBC-1 (Stoner, Quattro, & Weissman, 1997) were amplified by polymerase chain reaction (PCR) using specific primers under the following conditions: 95°C, 4min; 35x [94°C, 60s; 48-56°C, 60s; 72°C, 60s], 72°C, 45min; 10°C. The PCR-mix (10µl) containing 5µl 2x Taq PCR MasterMix (KT201, Tiangen Biotech Beijing.), 3.9µl DDW, 0.1µl primer mix and 1µl DNA. The PCR products were run on 1.5% agar in gel electrophoresis to determine amplification success. When unsuccessful, PCR was repeated several times using deviant protocols. 1µl of PCR product mix of four primers were added to 0.4µl LIZ 500 size standard and 8.6µl formamide per well and were analyzed using the ABI-PRISM 310 sequencer. The length of the microsatellite alleles was evaluated using GeneMapper 5.0 (Applied Biosystems) software package. Where in doubt, peak sizes were shared and discussed with other researchers and amplifications were repeated when necessary.

Data analysis

The analysis considered four geographic scales. On the local scale, each site (Des Moines, Edmonds, Shilshole and Shelton), was analyzed separately. On the regional scale, all four sites were simultaneously compared. On the US west coast scale, the results on the Puget Sound area were compared with two southern sites in central California: Santa Cruz (Reem et al., 2013a) and Moss Landing (Karahan et al., 2016). Finally, allelic patterns from worldwide sites were compared on an international scale.
The alleles were binned using AutoBin v 0.9 (Guichoux et al., 2011) and checked for scoring errors using Micro-Checker v. 2.2.3 (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004). The frequency of null alleles was calculated in FreeNA using the EM algorithm (Chapuis & Estoup, 2007). Null alleles might bias the outcome of the analysis, therefore pairwise population differentiation (Fst) was calculated with and without the excluding null alleles (ENA) correction in FreeNA using 50,000 bootstrap iterations. A two-tailed t-test was used to compare these two Fst estimates.
Allele frequencies, private alleles (for each location, for each sampling year and for the entire list of samples), observed and expected heterozygosity, population differentiation (Dest), population structure (Fst) and the inbreeding coefficient (Fis) (both using the G-statistic approach), Nei genetic distances and Mantel tests were calculated in GenAlEx v.6.51b2 (Peakall & Smouse, 2012). The allelic richness was calculated in FSTAT (Goudet, 2003) using Hurlbert’s (1971) rarefaction approach as modified by Petit, Mousadik and Ponst (1998), while HP-Rare (Kalinowski, 2005) was used to calculate private allelic richness. The significance of the deviance in heterozygosity was tested using the GENEPOP online version (http://www.genepop.curtin.edu.au). Bayesian Analysis of Population Structure BAPS (Corander, Waldmann, Marttinen, & Sillanpää, 2004) and STRUCTURE v.2.1 (Pritchard, Stephens, & Donnelly, 2000) were used for the Bayesian structure analysis, and BAPS to determine the gene flow between different populations. In BAPS, we used the clustering groups of individual option for the analysis and in STRUCTURE, the length of the Burn-in period was set to 100,000 and the number of Markov-chain-Monte-Carlo repetitions to 106. The number of clusters, K , was allowed to range between one and ten and five iterations for each K were carried out. The optimal K was chosen using the Structure harvester (Earl & vonHoldt, 2012) and the different runs aligned using CLUMPP (Jakobsson & Rosenberg, 2007). This step was repeated three times to evaluate if the number of clusters is consistent. Population clustering should be used with caution due to the different assumptions of the programs (Sinai et al., 2019), therefore NetStruct was used as a third method for analyzing the population structure using the network theory (Greenbaum, Templeton, & Bar-David, 2015). After an initial run, the threshold was set between 0.2 and 0.4 and 999 permutations were carried out. The genetic drift (Fs’) was calculated in TempoFS (Jorde & Ryman, 2007). Pairwise Fst, a measure of population structure that is commonly used for showing population differentiation, was calculated in Arlequin (Excoffier, Laval, & Schneider, 2005).