Biological invasions have become more frequent, driven by increasing globalization, anthropogenic activities, and climate change (Hulme, 2009;Van Kleunen et al., 2015). They cause severe negative impacts on local biota while also imposing economic consequences (Gallardo et al., 2016;Walsh et al., 2016). The footprints of biological invasions can be traced to changes in species diversity, dramatic alterations in communities and habitats, top-down and bottom-up control modifications, shifts in food chains, changes in nutrient cycling, and the attenuation of ecosystem services (Anton et al., 2019;David et al., 2017;Simberloff and Rejmánek, 2011). Similar to terrestrial ecosystems, biological invasions in the marine or oceanic realm are major drivers of ecological and evolutionary shifts, altering community structures and restructuring ecosystem functions with direct and indirect impacts on ecosystem services (Carlton and Geller, 1993;Darling et al., 2017;Katsanevakis et al., 2014). With the increasing impact of biological invasions in marine environments, more interest is being shown in the inclusion of genetic, phylogenetic, and evolutionary aspects of research, with parameters that may improve the resolution and cost effectiveness of monitoring biological invasions (Darling et al., 2017;Rius et al., 2015) and describe changes in population genetics and adaptation properties of important invasive species in greater detail (Barrett, 2015;Tepolt, 2015).
Botryllus schlosseri (Pallas, 1766) is a common colonial ascidian species in the Mediterranean Sea and European Atlantic (Reem et al., 2022) that has been introduced to temperate zones worldwide (Freeman et al., 2016;Lin and Zhan, 2016;Lord, 2017;Reem and Rinkevich, 2014). It is currently found on all continents except Antarctica, including the coasts of Japan, New Zealand, India, South Africa, Chile, (Reem et al., 2013a) was relatively isolated, under high genetic drift, and further characterized by high mutation rates. The Moss Landing population, studied between 1996 and 2008 (Karahan et al., 2016), was affected by episodic freshwater floods and subsequent recovery.
These long-term studies on Botryllus schlosseri populations (Karahan et al., 2016;Reem et al., 2013a) detailed the population genetic parameters of two populations that had already been established for several decades. It is therefore of great interest to compare their parameters with those of recently established B. schlosseri populations, such as those on the US North Pacific coast (Puget Sound, Washington), where in the late 1980s, this species was missing or only just introduced. We studied the genetic parameters of four B. schlosseri populations residing in marinas in the Seattle, WA area, separated by up to 124 km coastline, during a period of up to 19 years (exact patterns of introduction are unknown). Thus, a major aim was to gain a perspective on possible fluctuations in genetic structures on a regional scale in recently introduced populations of this globally invasive species.
The Puget Sound (Fig. 1) is part of the Salish Sea, an estuarine system in Washington State, USA, which is linked to the Pacific Ocean via three connections to the Strait Juan de Fuca. Four marinas were studied: three are located in the central basin (Edmonds (47 • 48 ′ 26 ′′ N, 122 • 23 ′ 34 ′′ W), Shilshole (47 • 40 ′ 52 ′′ N, 122 • 24 ′ 18 ′′ W), and Des Moines (47 • 24 ′ 2 ′′ N, 122 • 19 ′ 47 ′′ W) marinas), as well as in Shelton Yacht Club (47 • 12 ′ 52 ′′ N, 123 • 5 ′ 5 ′′ W) 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_ des_moines_marina_assessment_draft_3-15-19.pdf).
The Shilshole Marina has approximately 1500 mooring slips, and it is the largest marina in this study. The City of Des Moines Marina usually has anywhere between 400 and 800 boats moored, with a constant flow of 10-30 daily visitors (Tara Reilly, Office Specialist of the City of Des Moines Marina, pers. comm.). Edmonds Marina has 668 mooring slips and receives approximately 3100 visiting boats per year. The remote Shelton Yacht Club offers only 109 mooring slips and is surrounded by extensive shellfish aquaculture.
Botryllus schlosseri samples were collected every two years from 1999 to 2007 and in 2013 and 2018 (seven collection sessions; Suppl. Table A1). Colonies or colony fragments separated by at least 1 m were removed using single-edge razor blades and placed into separate 1.5 ml vials with 100 μL lysis buffer. The DNA extraction was carried out according to the protocol in (Paz et al., 2003). In brief, 20 μl NaClO 4 and 120 μl phenol/chloroform/isoamyl alcohol (25:24:1 v:v:v) were added to the homogenized samples. After mixing and centrifugation (10,000×g for 5 min), the aqueous solution was transferred into another tube, and 120 μl of chloroform-isoamyl alcohol (24:1 v/v) was added. Mixing and centrifugation were repeated, and DNA was precipitated by adding 100% ethanol and centrifuging (10,000×g for 15 min at 4 • C). The DNA was washed with 70% ethanol, allowed to dry, and then dissolved and diluted (1:20) with double distilled water.
Eight B. schlosseri microsatellites (Suppl. Table A2), BS-2, BS-8, BS-9 (Abdoullaye et al., 2010), BS-811 (Pancer et al., 1994), PB-29, PB-41, PB-49, and PBC-1 (Stoner et al., 1997), were separately amplified using polymerase chain reaction (PCR) with specific primers under the following conditions: 95 μl double-distilled water 0.1 μl primer mix, and 1 μl DNA. The PCR products were run on 1.5% agar for gel electrophoresis in order to determine amplification success. When unsuccessful, PCR was repeated several times using deviant protocols with different annealing temperatures (within the range of 48-56 • C) or the protocol described in (Bock et al., 2011). Next, 1 μl of PCR product mix of four primers was added to 0.4 μl LIZ 500 size standard and 8.6 μl formamide per well and analyzed using the ABI-PRISM 310 sequencer at the Technion's Biomedical Core Facility, Haifa, Israel. The length of the microsatellite alleles was evaluated using GeneMapper 5.0 (Applied Biosystems) software package. When in doubt (<5% of the cases), peak sizes were shared and discussed with other researchers, and amplifications were repeated when necessary. The majority of DNA extractions and genotyping were carried out in spring 2019 to ensure consistent amplification and genotyping success.
The analysis considered four geographical scales. At the local scale, each site (Des Moines, Edmonds, Shilshole, and Shelton) was analyzed separately. At the regional scale, all four sites were compared simultaneously. On the US west coast scale, the results for the Puget Sound area were compared with those of two southern sites in central California, Santa Cruz (Reem et al., 2013a) and Moss Landing (Karahan et al., 2016), both of which were subjected to the same methodology as that of the current study. Finally, allelic patterns from sites worldwide were compared at the international scale. Here, a "population" was considered to be the collected individuals from any single sampling location in a given year, unless otherwise stated.
The raw allele sizes were binned using AutoBin v 0.9 (Guichoux et al., 2011), a program that sorts alleles by size to detect relevant size changes for the proposed alleles and checks for scoring errors using Micro-Checker v. 2.2.3 (Van Oosterhout et al., 2004). The frequency of null alleles was calculated using FreeNA with the expectation maximization (EM) algorithm (Chapuis and Estoup, 2007). Null alleles might bias the outcome of the analysis; therefore, pairwise population differentiation (Fst) was calculated with and without the exclusion of null allele (ENA) correction in FreeNA using 50,000 bootstrap iterations (Aglieri et al., 2014). A two-tailed t-test was used to compare the 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), and Nei's genetic distances were calculated and Mantel tests were performed using GenAlEx v.6.51b2 (Peakall and Smouse, 2012). For the Mantel tests, the temporal distance (years) was used for within-site analyses, and the geographical distance (km) was used for between-marina analyses. All of these were computed with both Fst and Nei's genetic distance indices, as both have been used previously (for example, Pineda et al., 2016;Reem et al., 2013a). Allelic richness was calculated in FSTAT (Goudet, 2003) using the rarefaction approach (Hurlbert, 1971) as modified by Petit et al. (1998), while HP-Rare (Kalinowski, 2005) was used to calculate private allelic richness. The significance of deviation in heterozygosity was tested using the online version of GENEPOP (http://www.genepop.curtin.edu.au). Pairwise Fst, a measure of population structure commonly used to show population differentiation, was calculated in Arlequin (Excoffier et al., 2005). We further grouped samples from the earlier phase (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and then the later (2013-2018) phase of collection to reveal differences between early and more recently established populations. Bayesian analysis of population structure (BAPS) (Corander et al., 2004) andSTRUCTURE v.2.1 (Pritchard et al., 2000) were used for the Bayesian structure analysis, and BAPS was used to determine the gene flow between different populations, both on a temporal and spatial scale. In BAPS, we used the clustering groups of individual options for the analysis, and in STRUCTURE, the length of the burn-in period was set to 100, 000, while the number of Markov chain Monte Carlo repetitions was set to 10 6 . The number of clusters (K) was allowed to range between one and ten, and five iterations were carried out for each K. The optimal K was chosen using STRUCTURE HARVESTER (Earl and vonHoldt, 2012), and the different runs were aligned using CLUMPP (Jakobsson and Rosenberg, 2007). This step was repeated thrice to evaluate whether the number of clusters was consistent. Population clustering should be used with caution owing 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 network theory (Greenbaum et al., 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') between different years and locations was calculated using TempoFS (Jorde and Ryman, 2007).
Altogether, 533 samples were collected (521 analyzed) over a period of 19 years, between 1999 and 2018, in four locations in the Puget Sound area, WA, USA. The detailed allele frequencies of all eight microsatellite loci are summarized in Suppl. Table A3. A two-tailed t-test between the Fst values (with and without ENA correction) did not show significant results (p = 0.463), indicating that null alleles did not influence the results. The inbreeding coefficient (Fis) was highly significant in all populations.
In total, 190 tissue samples were collected from seven sampling dates (years 1999-2018; Suppl. Table A4). All eight microsatellite loci were polymorphic, with 71 alleles in total. About one third of the total alleles (23/71) appeared in only a single sampling period: eight private alleles were found in 2018, three in each of the sampling years 1999, 2003, and 2013 and two in 2001, 2005, and 2007. Null alleles were detected at all loci, with high occurrence rates (>0.1) in PB-41 (0.284), BS-811 (0.255), PB-49 (0.229), BS-8 (0.18), PB-29 (0.175), and BS-9 (0.126). The genetic diversity indices are summarized in Table 1 andreveal fluctuating allelic richness and significant (p < 0.05) heterozygote deficiency, indicating that the populations were not in Hardy-Weinberg equilibrium (HWE). The highest peak in allelic richness occurred in 2018, followed by another peak in 2003, both of which coincided with the peaks in private allelic richness.
The Mantel test showed no significant correlation (p = 0.438) between either Fst and Nei's genetic distance or the time between the different sampling periods (Suppl. Table A5). Overall, the annual Des Moines population differed significantly (p = 0.001) (Table 2). Interestingly, the Dest of the populations from 1999 to 2007 exceeded that of the 2013-2018 populations by 0.064 (p = 0.001) and 0.024 (p = 0.062). The estimated genetic drift (Fs') was even stronger between 1999 and 2007 (0.076) than between 2013 and 2018 (0.048). Pairwise population differentiation (Fst) values showed that samples from all but three pairs (the years 2005/2013, 2007/2013, and 2005/2018) were significantly different (p < 0.05) (Suppl. Table A6A). Despite the lack of similarities between the years, the BAPS analysis revealed only a single cluster. However, two clusters were identified using STRUCTURE (Fig. 2A). In 1999 and 2001, cluster 1 (red in Fig. 2A) was dominant, whereas the later years were more mixed, with a slight dominance of cluster (yellow in Fig. 2A) in 2007 and 2013.
In total, 152 tissue samples were analyzed from the Shilshole marina. Sample collections were attempted seven times from 1999 to 2018, but 1999 was not included in the analysis because only a single colony was observed (not collected). It should also be noted that the 2005 collection yielded only 13 samples, despite increased collection efforts. All microsatellites were polymorphic, representing only a single BS-8 allele in five of the sampling periods (2001, 2003, 2005, 2007, and 2018) and only one allele of BS-9 in 2005. A total of 63 microsatellite alleles were found, with 15 private alleles occurring during a single sampling period (Suppl. Table A4). Four private alleles occurred in 2018, three in and 2007, two in 2005 and 2003, and one in 2001. The null alleles were frequent (>0.1) at BS-811 (0.302), PB-41 (0.270), PB-49 (0.141), PB-29 (0.127), PBC-1 (0.126), and BS-9 (0.101). The summarized genetic indices (Table 1) revealed a gradual decline in gene diversity (He) and that the populations were not in HWE. The highest value of allelic richness was reached in 2005, which was slightly higher than that in the two preceding sampling years, whereas the peak in private allelic richness was reached in 2007. Both Fst and Nei's genetic distance indices and the time between the sampling periods were significantly correlated (p = 0.02 and p = 0.041, respectively) (Suppl. Table A5). Each run in STRUCTURE resulted in four clusters (Fig. 2C). Cluster 1 (blue in
Four sampling periods yielded a total of 101 tissue samples. In an additional sampling session (1999), no colonies were found on the marina's hard-bottom shallow substrates. All microsatellite loci were polymorphic; however, only a single allele was found in 2003 and 2018 at locus BS-9. Among the 55 alleles (Suppl. Table A4), 17 were private alleles: two private alleles in 2003 and five in each of the other sampling years. The null alleles were frequent at PB-41 (0.351), BS-811 (0.279), PB-29 (0.172), BS-8 (0.166), and PB-49 (0.108). Genetic indices (Table 1) revealed an increase in allelic richness and gene diversity from 2003 to 2005, with a peak in allelic richness and private allelic richness in 2005. There was no correlation between Nei's genetic distance or Fst and the time between the sampling periods (Suppl. Table A5). Pairwise Fst values suggested a significant difference between the samples from the years 2003/2005, 2003/2018, and 2005/2018 (p < 0.05) (Suppl. Table A6C). The weak Dest was significant (p = 0.006), and the Fs was moderate (Table 2). The BAPS analysis suggested two clusters: one for 2003 and the second for 2005-2018. Nevertheless, the inconsistent K number suggested by STRUCTURE showed a very weak population structure. Two clusters were finally assigned owing to the constant (but not always tall) Delta K peaks at K = 2 (Fig. 2D). The BAPS results were confirmed using STRUCTURE, with 2003 being assigned to a single cluster, whereas the other three sampling dates failed to group into a specific cluster.
Shelton marina was visited during only three sessions, yielding 78 tissue samples. All microsatellite loci were polymorphic, although just a single allele was represented each in 1999 at loci BS-8, BS-9, and PB-41, in 2003 at locus BS-9, and in 2018 at locus BS-8. Forty-six alleles were recorded, of which 19 were private (e.g., found on a single sampling date) (Suppl. Table A4). Most private alleles occurred in 2003 (n = 9), followed by 2018 (n = 7) and 1999 (n = 3). Null alleles were frequent at four loci (PB-41, PB-49, BS-811, and PB-29), ranging from 0.220 to 0.268. The genetic indices (Table 1) revealed peaks in allelic richness and private allelic richness in 2003. Fst and Nei's genetic distance indices were not correlated with the time elapsed between the sampling events (Suppl. Table A5). STRUCTURE suggested three clusters (Fig. 2B); however, no meaningful trend was detected. BAPS, on the other hand, grouped the years 1999 and 2003 into one cluster and 2018 as the second. The different populations were weakly (p = 0.014) separated according to Dest (Table 2), and Fs was moderate. The pairwise Fst showed significant differentiation between the samples from 1999/2018 and 2003/2018 (Suppl. Table A6D), which is in line with the clusters suggested by BAPS.
At a regional scale, 104 alleles were recorded during the 19-year sampling period. Nineteen private alleles were found in Des Moines (26.8% of local alleles), eight (17.4%) in Shelton, nine (14.5%) in Shilshole, and eight (14.5%) in Edmonds. Analyses of the three most frequent alleles per locus (Fig. 3) revealed that only a single allele at loci BS-8 and BS-9 (181 bp and 194 bp, respectively) was dominant every year in each location, whereas loci BS-811, PB-29, and PBC-1 showed variable allele frequencies. In Edmonds, the third most frequent allele at locus BS-811 was frequent in 2003 and 2018 but absent in the years between. Only eight alleles were present in every sampled population, alleles 186 bp and 189 bp at BS-2, alleles 152 bp and 156 bp at PB-29, and alleles 199 bp, 202 bp, 205 bp, and 210 bp at PBC-1. Some alleles reflected site-specific distributions. For example, an allele of 178 bp at microsatellite BS-2 was present in all Shilshole and Shelton populations but was absent in Des Moines and Edmonds populations. In contrast, allele 192 bp at microsatellite BS-8 was present during every sampling period in Des Moines and Edmonds populations but completely absent in Shilshole and Shelton populations.
The Des Moines populations were almost twice as differentiated (Dest) as the populations of the other three marinas (Table 2), a result further supported by the population structure Fst, revealing the highest values for Des Moines. For both the indices, the differences between the four locations (Dest and Fst) were higher during the earlier period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) than between 2013 and 2018. For the inbreeding (Fis) index, only Shelton showed remarkably low values. The overall genetic drift values were similar in Des Moines, Edmonds, and Shelton and less than half in Shilshole. For allelic richness, the highest numbers were assigned to Des Moines (Table 1), followed by Shilshole, Edmonds, and Shelton. Des Moines also showed the highest He, which is also a measure of the evenness of allelic frequencies, similar to the number of effective alleles (Brown and Weir, 1983). These high levels were followed by those of Edmonds, Shilshole, and Shelton. The physical distance between marinas and Nei's genetic distance or Fst between the sites was not correlated (Suppl. Table A5).
Analysis of all 20 populations simultaneously resulted in four BAPS individual clusters, where each location was assigned to a separate cluster (Suppl. Fig. A1). STRUCTURE, however, suggested two clusters (Fig. 4A). The Des Moines samples were mainly assigned to cluster 1 (yellow in Fig. 4A), whereas the Shilshole samples belonged primarily to cluster 2 (red in Fig. 4A), with individuals from Shelton and Edmonds mixed together. Netstruct suggested three highly significant clusters (Fig. 4B), with Des Moines belonging mainly to cluster one (blue in Fig. 4B), Shilshole to cluster two (red in Fig. 4B), and Shelton to cluster three (green in Fig. 4B), whereas Edmonds appeared to be well mixed. Gene flow between sites during the study period was somewhat restricted (Suppl. Fig. A2), as Shelton received gene flow just from Shilshole and thus emerged as the most isolated site, with only 0.02 of gene flow from other locations, compared to 0.06 at Edmonds and Shilshole and 0.07 in Des Moines. The strongest gene flow was from Shilshole to Edmonds (0.031), followed by that from Des Moines to Shilshole (0.023). Edmonds and Shelton were not connected at all, while there was a unidirectional connection from Shelton to Des Moines (0.022) and bidirectional exchange with Shilshole. Comparing only the sites without yearly divisions showed a highly significant pairwise Fst for all pairs (Suppl. Table A6E).
Des Moines and Shilshole exhibited a similar number of private alleles on most sampling dates and generally the highest numbers (Table 3). However, in 2005, Edmonds had more private alleles than Des Moines and Shilshole. Interestingly, Shelton had fewer private alleles than Des Moines in 1999 and all other sites in 2018. On a yearly basis, pairwise Fst between years for the whole Puget Sound area showed that samples from all combinations except 2001/2003, 2003/2005, and 2007/2013 were significantly differentiated (Suppl. Table A6F). In both 2003 and 2018, the geographical distance was not correlated with either Nei's genetic distance or Fst between the marinas (Suppl. Table A5).
The analyses included the Puget Sound sites, and two Californian sites, Santa Cruz (Reem et al., 2013a) and Moss Landing (Karahan et al., 2016), involving the five shared microsatellite loci (BS-811, PB-29, PB-41, PB-49, and PBC-1) used for analyses in these locations.
Allelic richness and expected heterozygosity in the Puget Sound populations were lower than those in the two Californian populations (Table 4). Despite similar allelic richness and expected heterozygosity, the Moss Landing and Santa Cruz populations were clustered into two different groups (orange and red in Fig. 5A) in STRUCTURE, and all four Puget Sound populations formed a third cluster (yellow in Fig. 5A). BAPS created three distinct clusters when the four Puget Sound sites were considered as a single location and no gene flow between any of the sites was discernible (Suppl. Fig. A4). Three highly significant clusters were selected using Netstruct for the west coast analysis (Fig. 5B). Strikingly, the Santa Cruz populations clustered together as a single entity, whereas the Moss Landing populations were mostly associated with the remote Seattle populations, showing minimal similarities to the closer Santa Cruz cluster.
During four sampling years (1999, 2001, 2005, and 2007), DNA samples were taken in at least one Californian and two Puget Sound sites. Analyses for each year separately revealed only a single genetic connection between California and Puget Sound ( 2007) and a substantial gene flow (0.039) from Des Moines to Moss Landing, whereas Santa Cruz, Shilshole, and Edmonds remained isolated (Fig. 6).
For global scale analyses, we used ten sites (Table 4) and five shared microsatellite loci (BS-811, PB-29, PB-41, PB-49, and PBC-1). Locus BS-811 had the highest number of alleles at all the locations (20-65; Table 4). When considering the total number of individuals analyzed, Puget Sound emerged as the region with the lowest number of alleles per individual at either locus and at the combined loci. This places the whole Mediterranean area ( (Reem et al., 2017), Moss Landing (Karahan et al., 2016), and South America (Ben-Shlomo et al., 2010) as containing the highest numbers of allele per 100 colonies (N = 57-66). The number of alleles in the Puget Sound was notoriously low (N = 16), primarily at locus PBC-1, presenting only 61.5% of alleles as compared to those of the second lowest site, Moss Landing. The European and Mediterranean locations (without Scandinavia) showed the highest number of alleles at PB-29, PB-41, and PBC-1 and were also leading in the number of alleles among the other two loci. The average expected heterozygosity among the Puget Sound populations was by far the lowest across the globe, suggesting a limited evenness of allele frequencies. Furthermore, allelic richness was among the lowest in Puget Sound, whereas the highest values were observed in Israel and Santa Cruz.
Analysis of frequent alleles (>0.1 in at least a single population; Suppl. Table A7) revealed some disparities among sites worldwide. For example, all the frequent alleles on locus PB-29 on the west coast were <160 bp, while in the other regions, at least one allele was >160 bp. At the most polymorphic locus, BS-811, most of the European, South American, and New Zealand frequent alleles were <260 bp in contrast to the North American alleles (<300 + bp). Locus PB-41 did not display different patterns among the locations, except for Israel (Paz et al., 2003). At most sites, the majority of frequent alleles were <180 bp; however, most alleles in the Israeli populations were >180 bp. For locus PB-49, allele sizes >240 bp were just found in Santa Cruz (Reem et al., 2013a), the US east coast, (Stoner et al., 2002), and the Mediterranean Sea (Reem et al., 2017), whereas in Israel (Paz et al., 2003), Scandinavia, (Reem et al., 2013b)
Botryllus schlosseri is a marine invertebrate (Berrill, 1950;López-Legentil et al., 2006;Paz et al., 2003;Reem et al., 2017) found in the Mediterranean sea and European Atlantic shallow waters (down to 200 m depth) (Berrill, 1950;López-Legentil et al., 2006;Paz et al., 2003;Reem et al., 2017) that has become a cosmopolitan alien species on man-made submerged hard-bottom substrates, primarily inhabiting marinas and harbors all over the temperate zones of the northern and southern hemispheres (Berrill, 1950;Ben-Shlomo et al., 2001, 2010;Lambert, 2001;Reem et al., 2013aReem et al., ,b, 2017)). This species is also considered to be a fouling pest that is associated with economic costs in aquaculture (Arens et al., 2011;Carver et al., 2006). Traits such as fast adaptation to man-made environmental conditions (Lambert and
A summary for five microsatellite loci in 10 assigned populations worldwide. Area: Location of the study; N: Number of studied samples; BS-811, PB-29, PB-41, PB-49, PBC-1: number of alleles found at each locus (number alleles/N*100); Total: number of alleles over all loci; NA/100 specimens: total number of alleles per 100 specimens; AR: allelic richness (where provided); He: expected heterozygosity; Sources: 1: current study; 2: Karahan et al. (2016);3: Reem et al. (2013a);4: Stoner et al. (2002); 5: Reem et al. (2017);6: Paz et al. (2003);7: Ben-Shlomo et al. (2006);8: Reem et al. (2013b);9: Ben-Shlomo et al. (2010);10: Ben-Shlomo et al. (2001). Lambert, 2003;Lambert, 2001) and enhanced mutation rates that increase genetic variability in newly established populations (Reem et al., 2013a) further contribute to B. schlosseri's invasiveness in a way that, upon the establishment of pioneering colonies, this species quickly spreads to become one of the most common species in the hard-bottom invertebrate consortia (Lambert andLambert, 1998, 2003;Locke et al., 2009;Martin et al., 2011). Many B. schlosseri populations beyond the Mediterranean and European Atlantic coasts have been established for more than five decades (e.g., in Africa (Millar, 1955) and South America (Orensanz et al., 2002)) and up to a century and more, as in the US east coast (Gould, 1841), California, (Van Name, 1945), and New Zealand (Van Name, 1945). Puget Sound populations are likely among the youngest regional populations of this species. Despite early mentions of B. schlosseri in a naval shipyard in the Seattle area (US Navy, 1951), this species has not been recorded in civil marinas for several decades. Indeed, the partial absence in our 1999 survey may indicate a very recent establishment and coincides with the finding of B. schlosseri in Des Moines and Shelton (Cohen et al., 1998), but not in Edmonds. Thus, Puget Sound provides a unique opportunity to study relatively recently established (ca. three decades ago) B. schlosseri populations.
The results of the present study reveal significantly differentiated populations on both a spatial and temporal scale. Moderate genetic drift (ranging from 0.027 to 0.071) within populations and limited gene flow (up to 0.031) between populations from different locations were also observed. Population clustering showed that the Des Moines and Shilshole populations were assigned to different clusters, whereas the Edmonds and Shelton populations represented a mix of the clusters. The low number of microsatellite alleles found in the Puget Sound populations compared to that of other sites worldwide is remarkable. These results differed from a study on European harbor populations of the solitary ascidian Styela plicata, which revealed differences only on a spatial scale but not on a temporal scale (Pineda et al., 2016).
The absence of B. schlosseri in two surveys conducted in Edmonds marina in 1998 (Cohen et al., 1998) and 1999 (this study) made this marina a candidate for the most recently established site (out of the four studied) for the B. schlosseri Puget Sound populations. Four years later (2003), the allelic indices of this newly established population significantly differed from those of the populations from the subsequent sampling period (two years later). In contrast, the Shilshole populations showed fewer fluctuations in allelic indices between the sampling points, despite only a single colony being found in 1999. This indication of a stable population is further supported by genetic drift, which was much lower at Shilshole than at the other three sites. Genetic indices showed diverse behaviors in the marinas studied. While the allelic richness (AR) and gene diversity (He) were highly fluctuating in Des Moines, Edmonds, and Shelton, these two indices slowly declined over time in the Shilshole marina. The two northernmost marinas, Edmonds and Shilshole, showed a peak in allelic richness during the 2005 sampling period, whereas in Des Moines and Shelton, the highest numbers were observed in 2018 and 2003, respectively. Only in Shilshole did the peaks of He ( 2003) and private allelic richness ( 2007) not coincide with those of the AR. The differences in genetic indices suggest that even close populations exhibit distinct genetic patterns. However, Shelton had by far the lowest values in genetic indices, all indicating its remoteness.
Heterozygote deficiency and significant inbreeding coefficients (Fis) found in all Puget Sound populations were in line with those of previous studies on the US west coast (Karahan et al., 2016;Reem et al., 2013a;Stoner et al., 2002) and worldwide populations (Ben-Shlomo et al., 2010;Lacoursière-Roussel et al., 2012;Paz et al., 2003;Reem et al., 2013bReem et al., , 2017)). Temporal fluctuations in allelic richness seem to be characteristic of B. schlosseri populations (this study) (Karahan et al., 2016;Reem et al., 2013a) and differ from long-term trends in isolated populations of other organisms, where a decreasing (García-Navas et al., 2015) or increasing (Jason Kennington et al., 2012) allelic richness was observed.
On the US west coast scale, the four Puget Sound populations were clustered together, despite more than 120 km separating the two furthest populations. In contrast, the Moss Landing (Karahan et al., 2016) and Santa Cruz (Reem et al., 2013a) populations in California were always assigned to separate clusters, despite being only 20 km apart. This indicates that, despite significant Fst and Dest values and different allelic richness peak timings, the Puget Sound populations at the west coast level remained genetically close. It is also of interest to find that Moss Landing populations are genetically closer to the Puget Sound populations than to the Santa Cruz populations, as revealed by the Netstruct clusters and gene flow charts, a result supported by (Reem et al., 2013a) with the assumption that Santa Cruz populations are isolated. Additionally, null alleles were frequent in both Puget Sound and Moss Landing populations (Karahan et al., 2016) but not in Santa Cruz populations (Reem et al., 2013a). The year 2007 was marked by the highest gene flow between the Puget Sound and Moss Landing populations, a year following a severe flood occurring in Moss Landing, where B. schlosseri was temporarily eradicated from that marina (Karahan et al., 2016). Thus, the genetic flow originating from the Puget Sound might have a footprint in the recolonization of the new population in the Moss Landing marina.
In the global comparison, Puget Sound's B. schlosseri populations had one of the lowest allelic richness and gene diversity values. These low numbers, even when compared to those of remote populations, such as of New Zealand (Ben-Shlomo et al., 2001) and South America (Ben-Shlomo et al., 2010), further illustrate the late establishment of the Puget Sound populations and can be explained by the limited time for the accumulation of mutations and inflow of genetic material from distant populations. The low allelic richness in Puget Sound could also reflect a low number of founding individuals. Furthermore, B. schlosseri microsatellite loci exhibited disparate repertoires of allele sizes and richness in different populations worldwide. While BS-811 had numerous alleles in all locations, other loci, such as PBC-1, revealed fewer alleles in regions recently invaded by B. schlosseri colonies. In addition, the amplification success of B. schlosseri microsatellites showed local differences even when following the same protocols. In the Puget Sound populations, the amplification success of BS-8, BS-9, and PB-41 was relatively low (up to 30% missing data), whereas the amplification of these loci in samples from the Israeli coast resulted in no failures (Tamir et al., 2022). Moreover, the native European populations (Reem et al., 2017;Reem and Rinkevich, 2014) revealed more alleles per microsatellite locus, further supporting the Mediterranean/Eastern Atlantic origin of B. schlosseri.
Studying newly established B. schlosseri populations over a period of 19 years in the Puget Sound revealed highly fluctuating patterns in genetic indices, such as allelic richness (AR), gene diversity (He), inbreeding coefficients (Fis), and population structures (Fst), showing a significant deviation from the HWE in all populations. Owing to these fluctuations, no temporal trend could be observed, and it is suggested that despite remarkable variations between the different sites, the Puget Sound populations remain isolated and still closely related. A comparison at a worldwide level revealed a considerably lower number of alleles in the Puget Sound population, supporting the recent introduction hypothesis.