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).