Rapid expansion of the invasive oyster Crassostrea gigas at its northern distribution limit in Europe: Naturally dispersed or introduced?

The Pacific oyster, Crassostrea gigas, was introduced to Europe for aquaculture purposes, and has had a rapid and unforeseen northward expansion in northern Europe. The recent dramatic increase in number of C. gigas populations along the species’ northern distribution limit has questioned the efficiency of Skagerrak as a dispersal barrier for transport and survival of larvae. We investigated the genetic connectivity and possible spreading patterns between Pacific oyster populations on the southern Norwegian coast (4 localities) and Swedish and Danish populations by means of DNA microsatellite analysis of adult oysters, and by simulating larvae drift. In the simulations we used a 3D oceanographic model to explore the influence of recent climate change (1990–2010) on development, survival, and successful spreading of Danish and Swedish Pacific oyster larvae to Norwegian coastal waters. The simulations indicated adequate temperature conditions for development, survival, and settlement of larvae across the Skagerrak in warm years since 2000. However, microsatellite genotyping revealed genetic differences between the Norwegian populations, and between the Norwegian populations and the Swedish and Danish populations, the latter two populations being more similar. This patchwork pattern of genetic dissimilarity among the Norwegian populations points towards multiple local introduction routes rather than the commonly assumed unidirectional entry of larvae drifted from Denmark and Sweden. Alternative origins of introduction and implications for management, such as forecasting and possible mitigation actions, are discussed.


Introduction
The Pacific oyster, Crassostrea gigas, was repeatedly introduced to Europe for aquaculture purposes in the second half of the 20 th Century (see [1] for a review), and has established wild  [23], B [17], C [3], D [24], E [9] & F [7]) and the genetic differenciation boundary between a documented southern and northern genetic group delineated by a dotted line ( B [17], G [15] & H [1]). The six C. gigas collection sites used in this study are indicated by the oyster symbole (See Table 1 for details). For Norway, valid and withdrawn aquaculture licenses for Ostrea edulis (http://www.fiskeridir.no/register/akvareg/?m=utl_lok&s=1; 20. May 2014) and C. gigas (Directorate of Fisheries) are indicated by open circles and stars, respectively. The map is produced using ESRIs GIS software ArcMap v species in northern Europe has raised a concern for further uncontrolled northwards expansion through massive larvae supply across Skagerrak from southern countries. This would cause severe problems for any mitigation actions against further northward spread of the species. In this study we used genetic analysis to investigate the origin of 4 established C. gigas populations along the Norwegian coast. We expect that if the main origin of the Norwegian populations is larvae dispersal from Swedish and Danish populations, then these populations would be genetically similar. Alternatively, if the origin is from post-introduction dispersal from local populations founded through other origins (e.g. aquaculture, shipping, or live trade), we expect these populations to be genetically different. We also examined what influence recent climate change and temperature conditions might have on dispersal of oyster larvae from Swedish and Danish populations, using a 3D oceanographic model, modelled sea water temperature for the region for selected years, and known temperature thresholds for larval development, spawning, and survival.

Sampling and DNA preparation
A total of 262 individuals of Crassostrea gigas oysters sampled from six Scandinavian populations in 2010 (Table 1 and Fig 1), were analyzed. As the Pacific oyster is an invasive species considered to be a threat to marine ecosystems in all three countries, permission was not required before sampling. The field studies did not involve endangered or protected species. Shell length varied between 6.5 and 18 cm with an average length of 10.9 cm among sampling locations, implying that all individuals were adults and probably from multiple generations. Oyster mantle samples were collected in 15 ml cap tubes and preserved in ethanol 96% (the ethanol was changed once). We used a new simple DNA preparation protocol without any purification step [25]. Briefly, individual tissue samples were washed in deionized water and about 5 mg mantle tissue was transferred to 100 μl 0.3% SDS with 2 μl proteinase K, incubated at 65˚C for 10 min followed by 98˚C inactivation for 2 min. The lysates were further diluted 10 −2 in Tris EDTA buffer (Fluka, Chemie GmbH, Switzerland) prior to performing PCR.

Genetic diversity and structure
Genotyping results were analyzed with Micro-Checker v2.3.3 [26] to identify potential inconsistencies and errors (e.g., null alleles and large allele drop-out). All incidences identified by Micro-Checker were chromatographically inspected before proceeding with further analyses. GenAlEx software v6.5 [27] was used to report overall observed (H O ) and expected (H E ) heterozygosity. Genepop v4.2 [28,29] was used to report observed (H O ) and expected (H E ) heterozygosity for each locus within each sampling location. The number of observed alleles (N A ) and calculated allelic richness (A R ), compensating for less individuals in N B (Rarefaction option), at each locus within each location was assessed using HP-RARE [30]. Independence among loci was tested by linkage disequilibrium (LD) and Hardy-Weinberg equilibrium (HWE) was calculated to identify loci and populations departing from theoretical equilibrium of allele frequencies, using the Arlequin software v3.5.1.3 [31]. Calculations of statistical significance were corrected for multiple tests according to the B-Y FDR method [32]. A multivariate Discriminant Analysis of Principal Components (DAPC) was used to resolve genetic connectivity between populations through sequential clustering and model selection. The DAPC was performed in R v3.3.2 [33] using the adegenet package [34] and the sampling locations as prior groups. The genetic relationship between samples, at the population level, was evaluated according to Chords distances (D CE ) [35], calculated and bootstrapped 2000 times with MSA v4.05 [36], and presented in a neighbour joining (NJ) tree [37] using the PHYLIP v3.68 software package [38] and SPLITSTREE v4.0 [39]. The significance of splits in the NJ tree were evaluated according to Hillis & Bull [40], i.e., all splits > 70% are considered statistically significant. Calculated pairwise D CE 's were also used to generate a Principal Coordinates Analysis (PCoA) in GenAlEx v6.5 [27]. Pairwise genetic differentiation was estimated by calculating the fixation index, F ST [41], and the statistical significance of the differences between populations was tested by 10,000 permutations of individuals between samples using MSA v4.05.

Larvae dispersal and survival simulations
Simulation of larvae dispersal and survival to settlement was performed with an open source, numerical 3D oceanographic model (ROMS) [42] with a spatial resolution of 800 m [43]. ROMS has shown accurate results when compared with field observations [44][45][46] and is a widely used model at both local and global scale (myroms.org). We focused on the influence of recent climate change by performing the simulations for 6 years (1990, 1998, 2002, 2006, 2007 and 2010) representing the climate since the 90s, one cold (2007), two warm (2002 and 2006) and one moderately warm year (2010). For survival, the simulated larvae have to had experienced 225 recruitment degree days [47], and the temperature at the landing site had to be ! 18˚C (according to Mann, Burreson [48]). Due to the high reproduction capacity (several million larvae per spawning individual), true individual-based modelling was impossible.
Hence we used the super-individual approach suggested by Scheffer, Baveco [49], where each modelled individual represents a large number of actual individuals. From each of 44 locations equally distributed along the Danish and Swedish coastline, 7 simulated larvae were released between 1 and 14 August (1 larvae every second day, i.e. 7 Ã 44 = 308 larvae per year [50]), and their floating path, experienced degree days and temperature at the landing sites were recorded. This could represent one viable super-individual from each of seven individuals in a small colony on each location. We did not include any behavior or random walk approach, but chose to distribute the time of release within the two first weeks of August, known to be the most relevant period for oyster spawning [12]. The number of landed larvae in Swedish and Norwegian coastal areas was counted within coastal grid cells of 50x50 km resolution for each simulated year, by summing the landed larvae within all the 800 m cells that fall within each of the coarser grid cells.

Genetic diversity
A total of eight microsatellite markers [15] were initially tested using a Norwegian oyster sample (N O ). Two of these microsatellites, CG108 and Cgsili29, failed to amplify and were not further used. The remaining six markers were first tested in simplex PCR to determine optimal annealing temperatures, and thereafter tested in combinations using various primer concentrations and cycling conditions for multiplex PCR testing. Both MgCl 2 and BSA were also tested as PCR helpers for multiplex optimization. Successful conditions were found for a fourplex and duplex PCR used in this study as described in the methods and Table 2.
Among the six microsatellite loci analysed, four amplified for all samples while CG49 and Cgsili44 each amplified in 259 out of 262 samples (Table 2 and S1 Appendix). All loci were polymorphic and the total number of alleles detected per locus varied from eight to 33 (Table 2 and S1 Appendix). Populations did not differ markedly in allelic richness, with average number of alleles ranging between 7.19, for N B which also had the lowest analysed sample size (n = 12), to10.35 for N G (Tables 1 and 3). The number of private alleles ranged from none at sampling location N B to eight at S S (Table 3). When the four Norwegian sampling locations are compared to a compiled Swedish and Danish sample, private alleles are 25 and 26 respectively. Genetic diversity was homogenous among the populations with the highest expected heterozygosity found in population N G (0.89), the lowest in N B (0.83), and highest and lowest Only one pair of loci (CGE009 & CG49) was identified as significantly linked by the LD analysis (B-Y FDR adjusted α = 0.00984). However, since this non-random association among alleles was not consistent among sampling locations and only identified in N O , loci CGE009 and CG49 were retained for further analysis. Significant departure from HWE was identified in 13 of the 36 tests (B-Y FDR adjusted α = 0.01198, Table 3). However, no locus showed significant departure from HWE in all sampled locations.

Population genetic structure
The PCoA analysis of all samples (Fig 2A) separated location N B (at the Norwegian west coast) from the remaining samples along the first dimension representing 36.5% of the total variance (72.2%). The remaining samples clustered into two groups along the second dimension, representing 18.6% of the total variance in the data set. The two clusters were made up of the Swedish (S S ) and the Danish (D A ) population in one group and the remaining Norwegian  Table 1.  (Fig 2A). The third dimension separated outer Oslo fjord (N O ) from the N I and N G samples, representing 17.1% of the variance. The same clustering pattern was also shown by the NJ tree (Fig 2B). In the NJ tree, sampling locations D A and S S were significantly split from the remaining sample locations. The N B sample occurred at an intermediate position within the NJ tree, between the other Norwegian localities (i.e. N O , N I and N G ), and the foreign countries localities (i.e. D A and S S , Fig 2B). For both the PCoA and the NJ tree, removal of location N B caused clustering of the remaining samples into three groups: 1) N O , 2) N I and N G , and 3) S S and D A (Fig 2c and 2d). The DAPC analyses, with all sampling locations included, showed a clear separation of N B from the remaining samples, primarily separated by the first principal component (Fig 3A). The second and third principal components did not vary in information value based on DA eigenvalues. With the removal of N B a stronger tendency for structuring between D A and S S versus the remaining Norwegian samples occurred (N I , N O , and N G , Fig 3B). Despite overlapping of individuals, the population ellipses for the Norwegian samples did not cross the centers of the D A and S S samples. The separation of sampling location N B was supported by high and significant F ST values when compared with the remaining sampling locations (Table 4). Among the remaining sampling locations, the pairwise F ST values showed non-significant genetic differences, except for the outer Oslo fjord location (N O ) versus the Danish and Swedish populations (D A and S S , respectively, Table 4). Hence all the statistical analyses indicate genetic differentiation.

Changes in dispersal and survival of Danish and Swedish oyster larvae
The simulation results showed that the water temperatures were too cold for the oyster larvae to develop and settle in Norwegian coastal waters in the 1990s (1990 and 1998), but warm summers since 2000 had adequate temperatures for development and survival of transported larvae across the Skagerrak (Table 5 and Fig 4). In the warmest year, 2002, a high fraction (36%) of the released larvae landed in Norwegian coastal waters, whereas in the following cold and moderately warm years (2007 and 2010) only 2-6% of the released larvae experienced sufficient water temperatures to successfully develop and settle on the Norwegian coast. In the two warm years (2002 and 2006), the simulated larvae could reach beyond the limit of the Skagerrak region and into the North Sea region. The hot spot for receiving the highest supply of oyster larvae along the Norwegian coast (i.e. the 50x50 km grid cell with the largest number of landed oyster larvae, Fig 4) included the sampling site at Hui in the Outer Oslo fjord (N O ).

Discussion
The analyses of the Pacific oyster (Crassostrea gigas) samples showed that the studied Norwegian populations were genetically different from the Danish and the Swedish populations. This contradicts the hypothesis that the Norwegian populations mainly origin from natural dispersal of larvae drifting from established populations in Denmark and Sweden [11].
Among the studied Norwegian Skagerrak C. gigas populations, sampling location N O differed genetically from the studied Swedish and Danish populations although it is located within the hot spot area for larvae supply as indicated by the larvae drift simulations. However, the identification of this area as a hotspot area for supply of larvae is otherwise supported by being the only area in Norway with oyster reef formation in 2015 [58]. The genetic difference identified between N O and the Swedish (S S ) and Danish populations (D A ) concurs with unpublished genetic studies [10] showing high similarity between Swedish and Danish populations, whereas the population at Hui (N O ) differed from the studied Swedish and Danish populations. The remaining two Norwegian Skagerrak populations (N I and N G ), located on each side of N O , form a cluster differentiated from both S S /D A and N O . This patchwork of dissimilar  [7,9], when all licenses for Pacific oyster aquaculture in Norway were revoked (Pers. comm. Directorate of Fisheries) it seems more likely that the identified genetic differences separating N I /N G from N O and from S S /D A are a consequence of multiple introduction events such as from aquaculture or shipping activities. On the other hand, this does not preclude the possible existence of other populations originating from larval drift across the Skagerrak.
The differentiation and independency of the studied Norwegian samples towards the Danish and Swedish samples is furthermore supported by the presence of private alleles in both groups in high and almost equal numbers (25 and 26, respectively). In order for the Norwegian populations to be a result of a frontier/range expansion (drift) scenario, as known from terrestrial [59] and marine organisms [60,61], the Norwegian Skagerrak group should have had lower genetic diversity, i.e. less private alleles, than the Swedish and Danish populations. Even among the Norwegian populations private alleles occur, providing further evidence of a lack of a uniform population structure, and pointing towards multiple introductions from separate sources of origin.

Origins of invasive Crassostrea gigas in northern Europe
The larvae drift simulations indicated that the water temperatures in the Skagerrak were too cold for larvae development, survival, and settlement to support successful natural dispersal of Danish and Swedish larva until 2000. Since 2000, summer temperatures have several years been sufficiently high for such natural dispersal and successful crossing of the Skagerrak barrier into Norwegian coastal areas. However, the lack of genetic similarity between the population at Hui (N O ), which is situated within the hot spot area for landing of foreign C. gigas larvae in the simulation study, and the studied Danish and Swedish oyster populations, indicates that so far there has been low success rate of this pathway. However, other genetic studies [10] have revealed some similarities between one Norwegian population (approximately 20 km north of N G ) that was not included in this study, and another Swedish population, which indicates the possible occurrence of successful recruitment of Swedish oyster larvae in Norwegian waters. Future climate change with rising summer temperatures is likely to increase the risk of C. gigas larvae dispersal [62]. Analysis of sea surface temperature data along the Swedish  (1990,1998,2002,2006,2007,2010), summed per coastal grid cell (50x50 km). Number of landed larvae (super-individuals) per grid cell is shown (see legend). The location and names of the sampled DNA stations in this study are indicated (black circles, cf. Table 1). For simulation details see [50]. Origins of invasive Crassostrea gigas in northern Europe Skagerrak coast [50] suggests a 125 km northwards displacement of the 19˚C temperature isocline in August. This implies a northward shift of the summer temperature needed to enable Pacific oyster spawning in wild populations. This will further push the distribution range of the Pacific oyster northward into previously unfavourable areas/ecosystems, as previously documented for other species [63]. Our findings of a theoretical possible increased supply of foreign C. gigas larvae in recent and future years clearly indicate a need to monitor and investigate the newly established populations of C. gigas along the Norwegian coast to assess the connectivity link across the Skagerrak area. There are many factors that may cause high pre-or post-settlement mortality of drifting C. gigas larva (e.g. predation, starvation, etc.) and that could counteract successful dispersal and colonization across the Skagerrak. In addition, selection imposed by strong environmental gradients, such as the temperature gradient in the studied region, promotes adaptive differentiation [64]. Local adaptation of earlier introduced C. gigas in Norwegian waters would imply that the local genotypes would have higher fitness than genotypes from foreign habitats [65]. Accordingly, recently landed foreign C. gigas larvae, would have lower chances of survival than locally adapted oysters.
The simulation model does not include any other mortality factors than the influence of temperature on the larvae's possibility to develop successfully during the planktonic phase, and sufficient temperature for the larvae to survive at the landing site. This implies that the predicted rate of success of transported Pacific oyster larvae, is likely to be higher than the real success rate. Other mortality factors in the planktonic phase (e.g. starving and predation), when settling (finding suitable substrate), and post settlement (including spatial competition with other species), will all reduce the larvae's chances of successful spreading. Hence the predicted rate of successful spreading in the two warm years, are likely to be higher than the actual rate because of these limitations, further reducing the possibility of connectivity between the populations.
The genetic differences found among the Norwegian Pacific oyster populations suggests that multiple introductions may have occurred along the coast. This could involve previous aquaculture activities or other introduction pathways such as shipping activities and live trade. Unfortunately, no aquaculture sources were included in the present analysis, so its role as a potential source of introduction cannot be established. Introduction of oysters to ports by shipping is possible since C. gigas larvae and adults have been found in ballast water and on ship hulls, respectively [66]. The Norwegian Skagerrak coast houses some large ports with high shipping activities and hence the potential for this introduction pathway exists [67]. Spreading of C. gigas from aquaculture sites has occurred in several countries [13,68]. In Norway aquaculture licenses for both native (Ostrea edulis) and the invasive oyster species (C. gigas) (Directorate of Fisheries, http://www.fiskeridir.no/register/akvareg/?m=utl_lok&s=1) have been given, indicating that introduction for aquaculture purposes is another plausible introduction pathway. Indeed, the Norwegian population from Bergen (N B , western Norway) was collected in the vicinity of a former aquaculture site (Espevik, Tysnes) for C. gigas, for which the origin of the imported larvae was reported to be Scotland [7]. This population, N B , showed strong genetic differentiation from the other studied populations. The possible aquaculture origins for the remaining sampled populations in Norway are difficult to establish from literature. Despite strong restrictions on the import of molluscs for cultivation purposes in Norway in 1986 [9], these aquaculture licenses were still assigned until 2001. The restrictions may have reduced the likelihood of recent repetitive aquaculture introductions.
Although being genetically different, the studied Norwegian oyster populations had low genetic diversity. This agrees with other studies [1,15,17] indicating a general pattern of low genetic diversity in the north. Among these studies [15] used the same six microsatellites as this study. Few differences in mean allelic richness were shown for all the analyzed sampling locations (Table 3), except for the westernmost sampling location, Bergen (N B ). Despite the relatively low number of analyzed individuals, the low allelic richness of the Bergen population could be due to a founder effect or subsequent bottleneck effects [69][70][71]. However, a bottleneck analysis [72,73] using the two-phase mutation model with default settings did not identify any bottleneck events in the analyzed samples (data not shown). Moreover, the simulation study indicates low chances for natural dispersal of Danish and Swedish larvae so far along the Norwegian coastline as to Bergen given recent year's climate. It therefore seems plausible that the Bergen population has adapted to local environmental conditions within the area since the 1970s when the species was introduced for aquaculture purposes. The clustering software STRUCTURE v2.3.4 [74] was run with and without N B (data not shown), to detect any population structure among the sampled populations. The program failed to show any structuring pattern. This concurs with previous C. gigas studies showing no population structuring except between the northern and the southern European populations [1] or between aquaculture and feral populations [75].
Initial introduction of C. gigas to Europe entirely originates from Japan, USA, and BC in Canada. Several studies [1,15,17] have demonstrated that the southern European populations genetically cluster with the introduction source populations. This indicates that the northern group has developed locally in Europe. As the history of introduction to the north (i.e. to Denmark, Norway and Sweden) mainly originate from the UK following the 1970s, it seems that genetic differentiation between the northern and southern group may have begun in the UK. Hybridization is known to be an evolution mechanism by which genetic variations may be swiftly introduced and fixed in populations, and in particular when species colonize new environments [76]. Considering that Crassostrea angulata has been introduced both prior to C. gigas, and in parallel into the UK [24], hybridization between the two species within the UK may have been possible during this spatial and temporal overlap. Hybridization is documented for C. gigas with C. angulata [77,78], supporting that such an event may have occurred. Although considered as conspecific to C. gigas by some authors [79], C. angulata has been shown to be sufficiently genetically distinct to be a separate species [78,80] and listed under the World Register of Marine Species (http://www.marinespecies.org/aphia.php?p= taxdetails&id=146900, July 2016).
This study contributes to the understanding of the genetic pattern of C. gigas in northern Europe and shows so far, low connectivity across the Skagerrak. Furthermore, it also demonstrates a likely future increase in successful dispersal of C. gigas larvae across the Skagerrak, as sea surface temperatures keep rising [14]. The current expansion might therefore temporarily be mitigated by reducing the density of the species in locations with suitable conditions for oyster growth and spawning, e.g. semi-enclosed bays, traditionally used for oyster aquaculture in Norway [7]. The suggested importance of aquaculture, shipping and import for live food as likely introduction pathways is highly relevant for nature management and implies the need to inform aquaculture industries and the public about the risk of introducing invasive species. The future risk of successful dispersal of Pacific oyster larvae from Danish and Swedish populations, due to climate change, emphasize the need for monitoring to detect any massive expansion as basis for targeted management of affected ecosystems. These conclusions may be extended to other invasive species with pelagic larvae stages also affected by climatic change.