Temporal Changes in Population Structure of a Marine Planktonic Diatom

A prevailing question in phytoplankton research addresses changes of genetic diversity in the face of huge population sizes and apparently unlimited dispersal capabilities. We investigated population genetic structure of the pennate planktonic marine diatom Pseudo-nitzschia multistriata at the LTER station MareChiara in the Gulf of Naples (Italy) over four consecutive years and explored possible changes over seasons and from year to year. A total of 525 strains were genotyped using seven microsatellite markers, for a genotypic diversity of 75.05%, comparable to that found in other Pseudo-nitzschia species. Evidence from Bayesian clustering analysis (BA) identified two genetically distinct clusters, here interpreted as populations, and several strains that could not be assigned with ≥90% probability to either population, here interpreted as putative hybrids. Principal Component Analysis (PCA) recovered these two clusters in distinct clouds with most of the putative hybrids located in-between. Relative proportions of the two populations and the putative hybrids remained similar within years, but changed radically between 2008 and 2009 and between 2010 and 2011, when the 2008-population apparently became the dominant one again. Strains from the two populations are inter-fertile, and so is their offspring. Inclusion of genotypes of parental strains and their offspring shows that the majority of the latter could not be assigned to any of the two parental populations. Therefore, field strains classified by BA as the putative hybrids could be biological hybrids. We hypothesize that P. multistriata population dynamics in the Gulf of Naples follows a meta-population-like model, including establishment of populations by cell inocula at the beginning of each growth season and remixing and dispersal governed by moving and mildly turbulent water masses.


Introduction
Marine planktonic organisms can grow extremely fast. Such fast growth, sustained by plentiful resources and temporally relaxed predation pressure, can lead to episodic, rapid and vast increases in their population sizes. The huge numbers of individuals and the moving and mixing water masses they inhabit are expected to foster large-scale population genetic homogeneity. Yet, a series of recent studies demonstrated that geographic structuring can occur in marine planktonic organisms [1][2][3]. In the case of the jellyfish Aurelia, trans-oceanic populations exist genetically in isolation-by-distance because the restricted lifespan of its planktonic medusa-stage prohibits gene flow across such extensive tracts of ocean [1].
Unicellular phytoplankton species usually show high genotypic diversity and in cases where genetically distinct populations are observed, they are often correlated with hydrographic or geographic features [1,[4][5][6][7][8][9][10]. Marine eukaryotic microalgae grow by means of mitotic division, but in contrast to daughter cells in macrophytes and animals, microalgal daughter cells disconnect and drift apart in their mildly turbulent environment, thus forming widely distributed clones. Episodic sexual reproduction in a population composed of large numbers of clones generates huge numbers of F1 cells with distinct genotypes, each of which in its turn can form a clone [11]. Therefore, the likelihood of sampling multiple individuals belonging to the same clone in a large phytoplankton population is very small, given the sample sizes normally deployed in population genetic studies [12].
Although the emergence of genetic differentiation without geographic barriers remains highly controversial, speciation can occur in sympatry [13] that is, if populations reproduce in distinct temporal windows, and/or have distinct ecological niches [14][15]. In phytoplankton, Casteleyn et al. [3] demonstrated that strains of the diatom Pseudo-nitzschia pungens established from cells collected in Belgian, Danish and Irish waters grouped into two genetically distinct, but apparently sympatric populations. Such genetic distinctness could merely be temporal, for instance resulting from contemporary establishment of founder populations from distinct sources, to be homogenized if sexual reproduction can still occur amongst them. In fact, marine habitats are among the most heavily invaded systems on Earth [16] and this is not necessarily restricted to invasions of alien species, but also to alien populations of resident species (e.g., [17]). Alternatively, mate preference and/or slightly offset bloom windows may keep these sympatric populations genetically segregated. If this is the case, then subtly different performance, e.g., different growth rates and environmentally governed differential mating success, could explain radical shifts in their proportions from one growth season to the next.
Few studies have addressed the structure of planktonic microalgal species over a temporal scale. A considerable genetic differentiation was detected over two consecutive years for the dinoflagellate Alexandrium fundyense in a coastal pond, where different populations were detected even amongst samples collected after 7 days. These highly diverse and dynamic patterns contrast with the constant genetic structure of the diatom Skeletonema marinoi, where samples composed of strains resulting from the germination of up to 100-years old resting stages collected from laminated sediment cores were found to belong to the same population as samples composed of strains obtained from the present-day plankton from the same site [18].
In the present study we explore population genetic structure of the planktonic, chain-forming diatom Pseudo-nitzschia multistriata at the Long Term Ecological Research station MareChiara (LTER-MC) in the Gulf of Naples (GoN, Tyrrhenian Sea, Mediterranean Sea, Fig. 1). The species has been recorded here from summer until late autumn of every year since its first appearance in the GoN, in 1996 [19][20][21].
Pseudo-nitzschia multistriata, as almost all diatoms, possesses rigid siliceous cell wall elements. During cell division each daughter cell inherits a parental cell wall element and makes a new one that fits precisely inside the parental one. Hence average cell size diminishes with ongoing mitotic division and the variance around the average increases. Populations of clonal cell lines escape from this miniaturization trap through sexual reproduction, upon which large cells (82-72 mm long) are formed within the zygote, thus reestablishing the initial large cell size. Yet, cells can be sexualized only within a restricted window of cell sizes [22]. A life cycle model of P. multistriata in the GoN by D'Alelio et al. [21] predicts that sexual reproduction in this species occurs once every two years. If the model is correct, then only six or seven sexual reproduction events have occurred between its first appearance in the GoN in 1996 and the onset of our population genetic sampling, in 2008.
To investigate intraspecific genetic diversity and population genetic structure of P. multistriata, we obtained seven polymorphic microsatellite markers [23] and used these to genotype strains generated from single cells, or clonal chains of cells, isolated from 22 net samples obtained during the species' seasonal appearance at the LTER-MC station over four consecutive years (September 2008 to July 2011). The genotype data were analyzed utilizing Bayesian clustering software to assess the most probable number of genetically distinct populations. Evidence emerged for the occurrence of multiple populations in the GoN, and therefore we assessed if these were separated in time, i.e., by sampling date, season or year, and if there were signs of them merging, e.g., by the increasing proportion of hybrids. In order to support the idea that sampled strains assigned as putative hybrids resulted from crossings between strains assigned to different populations, we included the genotypes of F1 offspring strains and their parental lineages, generated in Tesson et al. [11] in our analysis.

Ethics Statement
The present study on the planktonic diatom Pseudo-nitzschia multistriata was performed in accordance with the Italian laws and did not investigate endangered or protected species. No specific permissions were required to collect surface phytoplankton samples at the LTER-MC (40˚48.59N, 14˚159E, Fig. 1) in the Gulf of Naples.

Sampling site
The Gulf of Naples (GoN; Fig. 1) is an embayment open to the Tyrrhenian Sea (Western Mediterranean) and connected to the nearby Gulfs of Gaeta and Salerno (North-west and South-east, respectively) by slow seasonal cyclonic and anticyclonic eddies and offshore currents [24]. The sampling site, the LTER-MC, is located two miles offshore in the Gulf and is characterized by a surface temperature between 14˚C in February-March and 26˚C in August, a salinity of 37.5-38.0 (¡0.1) and the presence of a seasonal thermocline from April to September [19].

Cell collection, identification and culture maintenance
Surface phytoplankton samples were collected weekly using a rosette sampler equipped with Niskin bottles. Water samples for cell counts were fixed with neutralized formaldehyde at a final concentration of 0.8%. Depending on the total phytoplankton cell abundance of the samples, between 1 and 50 ml were allowed to settle in combined sedimentation chambers and cell concentration was estimated at 400x magnification [25].
For the isolation of living cells and chains, phytoplankton samples were collected with a standard phytoplankton net (20 mm mesh size). Cells and chains of Pseudo-nitzschia multistriata were isolated from 21 net samples collected between September 2008 and September 2010, and from a single net sample taken July 2011 (Table 1). Isolated cells were incubated each in 2 ml f/2 filtered medium at 22˚C, ca. 80 mmol photons?m 22 ?s 21 , and a photoperiod 12L:12D hour (Light : Dark). After a few days, isolation success and purity of the cultured strains were checked using an inverted light microscope; average cell size was determined by measuring cell length over the apical axis of at least 5 cells per strain at 200x magnification. After one week of growth, the obtained strains were transferred into culture bottles (Corning Flask, Corning Inc., NY, USA) containing 25 mL of f/2 medium to reach a concentration of about 10 3 cell?ml 21 . A total of about 1100 single cells or short chains were isolated, and of these, 735 (66.8%) were grown successfully in culture. Strain isolation success was comparable among sampling dates.

DNA extraction and strain genotyping using microsatellite markers
Genomic DNA of P. multistriata was extracted as described in Tesson et al. [23]. DNA-purity and quantity were assessed using the NANODROP Spectrophotometer (ND-1000, NANODROP, Wilmington, DE, USA) and agarose gel electrophoresis (1.5%). Pure DNA was diluted to 25 ng?ml 21 and 2 ml of the dilution was added to a well in a 96-well-plate containing the polymerase chain reaction (PCR) mix [23]. Seven microsatellite markers (PNm1, PNm2, PNm3, PNm5, PNm6, PNm7, and PNm16) were amplified in multiplexes by PCR, purified and combined for microsatellite identification as described in Tesson et al. [23] using forward labeled primers. Microsatellite reactions were prepared in automation with a robotic station BIOMEK FX (Beckman Coulter, Fullerton, CA) and analyzed on an Automated Capillary Electrophoresis Sequencer 3730 DNA Analyzer (Life technologies, 5791 Van Allen Way Carlsbad, CA 92008 USA). Electropherogram analysis was performed using PEAK SCANNER (version 1.0, Applied Biosystems). Peak assignment was refined manually in order to avoid scoring mistakes. In the few cases where stutter bands were present and interpretations were not univocal, samples were re-run up to three times to confirm scoring.

Microsatellite data analysis
MICRO-CHECKER (version 2.2.3 [26]) was applied to calculate the frequency of null-alleles (i.e., non-amplification of alleles) for each of the loci under the assumption that the loci are in HWE (S1 Table). A Brookfield-1 estimation was applied following recommendation by the program manual. Tests of presence of stutter bands due to slight changes in microsatellite length and low efficiency of large alleles amplification (large allele dropout) in PCR were performed as described in the manual (S1 Table).
Average number of alleles per locus, expected and observed heterozygosity, the inbreeding coefficient F IS , the fixation index F ST , and the G ST index of diversity across population were calculated for groups of strains using GENALEX (version 6 [27]). F IS is defined here as the proportion of the genetic variance in a group of strains contained in a sample. A significant positive F IS implies a significant inbreeding/homozygote over-representation, whereas a significant negative value, significant outbreeding/heterozygote over-representation; F ST , is defined here as the proportion of the total genetic variance contained in a group of strains relative to the total genetic variance among all of the strains; G ST , is a relative measures of F ST that takes into the account differences in sample size and the amount of genetic variation among populations versus all populations.
Allelic richness was re-calculated for each sample using GENCLONE [28]. In order to compare groups (i.e., samples, sampling years) with different numbers of strains, allelic richness was estimated in each group, with a resampling procedure, performing 1000 permutations for all possible numbers of sampling units (from 1  to N). Values obtained for the same number of samples were taken for comparison. Genotypic diversity (G/N) within samples was estimated using GIMLET (version 1.3.3 [29]). The number of genetically distinct populations and their occurrence through time were inferred using a Bayesian cluster analysis performed with STRUCTURE (version 2.1 [30]) The number of clusters (K; populations) was estimated following Evanno et al. [31] using the web-based program STRUCTURE HARVESTER [32]. Population structure was obtained applying a burn-in of 15,000 iterations, with 75,000 MCMC repetitions from K51 to K523. The parameters used to infer genetic structure were an ancestry model with admixture, along with an independent frequency model (infer l estimated to 0.1414). The LOCPRIOR assumption was not applied. The Assignment Probability (AP) threshold for a strain to a population assigned by STRUCTURE was set at 90%. A separate STRUCTURE analysis with the same settings was performed also including the genotypes of a set of parental strains and their F1 progeny obtained from Tesson et al. [11]. The existence of a bottleneck effect in the two populations and in consecutive years was tested using the Wilcoxon test implemented in the software BOTTLENECK [33], after 1000 iterations and hypothesizing the two-phased model (TPM) of mutation.
GENETIX (version 4.05.2 [34]) was used to identify linkage disequilibrium (10,000 permutations and a significance level of 0.05) within the whole data set. Factorial Correspondence Analysis was also performed with GENETIX, using genotypes assigned to populations as inferred from STRUCTURE and those not assigned to any of them. Principal Component Analysis (PCA) based on single strains as well as on groups of strains was performed using GENALEX and the pairwise matrix of Nei genetic distance. Statistical analysis and graphical representation of P. multistriata abundance and frustule size distribution -of cells observed in formalin fixed field samples as well as of cells from which the genotyped strains were raised -were performed using R [35].  Table 1). Twelve of these strains belonged to the frustule size class $55 mm, in which sexual reproduction cannot be induced, whereas the remainder belonged to the 55-30 mm size-interval, i.e. the size window in which sex can be induced (Fig. 3). We assessed cell size distribution of P. multistriata in a limited number of samples collected between 2008 and 2010 (Fig. 3). The cell size ranges observed in these natural samples matched those observed among the cells from these samples (Fig. 3).

Genetic diversity within samples and sample years
The 525 strains generated a clearly readable microsatellite fingerprint without any apparently missing values for any of the seven loci (S2 Table). All loci were polymorphic in all of the 22 samples, except locus PNm3, which was monomorphic in three samples (data not shown). The observed heterozygosity (Ho) was higher than the expected (He) -and often markedly so -in all samples except 8, 20 and 21. The inbreeding coefficient F IS was significantly negative for almost all of the samples from 2008 (samples 2-7) and for samples 15 and 16, implying significant outbreeding. The average inbreeding coefficient (F IS ) over all 22 sampling dates was negative (20.131¡0.039) indicating outbreeding and overall heterozygosity excess ( Table 1).
The number of alleles per locus ranged from 7 to 21. All loci, except PNm1 and PNm5, showed heterozygosity excess. No allele-drop out was detected in any of the loci (S1 Table). Only two loci (PNm2-PNm5) were in linkage disequilibrium within the whole data set (data not shown).
The strains exhibited 394 distinct genotypes (G). The genotypic diversity (G/N) ranged from 46.2% to 100% across the individual samples ( Table 1), from 66.2% for 2008, to 87.0% for 2010, with an average of 75.5% over all 22 samples ( Table 2). No clear relationship was detected between sample date and genotypic diversity, also because sample size varied markedly among sample dates.  Table 1).

Population genetic patterning
Results of the independent runs of K51 to 22 for assessing Bayesian population clustering as implemented in STRUCTURE, showed that the existence of two groups (K52; mean LnP(K)528683.8; SD LnP(K)51.78; DK5206.7) was by far the most likely (S1 Figure). The inferred groups are denoted from here on as populations POP_A and POP_B (green and red, respectively in Fig. 4, S2 Table). Strains assigned to POP_A dominated in 2008 and in the single sample of 2011, whereas the ones assigned to POP_B dominated in 2009 and 2010 ( Table 3). Since STRUCTURE results were obtained using the ancestry model with admixture, strains that could not be assigned to any of the two populations with a probability above the 90%-threshold were considered to be putative hybrids. These assigned hybrids constituted 15%, 27% and 24% of the strains sampled in 2008, 2009 and 2010, respectively, whereas no putative hybrids were assigned among the 2011samples. The proportion of assigned hybrids differed significantly among years (p50.02861 x 2 -test (x 2 57.11, df52)). The proportion of assigned hybrids did not differ significantly between 2008 and 2010 (p-value50.012), but both of these differed from that in 2009 (p-value.0.05).
Factorial Correspondence Analysis (FCA) of the strains' genotypes resulted in a plot in which all the strains were distributed along a stretched-out swarm without any apparent subdivision (Fig. 5). Strains assigned to POP_A grouped towards the right side, those assigned to POP_B towards the left side, and those classified as putative hybrids, lay in-between. PCA mapping of the 22 samples resulted in a plot (Fig. 6) in which the location of the samples corroborated the findings of STRUCTURE (Fig. 4). In addition, patterns not captured by STRUCTURE were uncovered. The 2008-samples, composed of strains mostly assigned to POP_A, clustered together on the lower right of  Fig. 6) were still dominated by POP_B strains, but, contained a minority of strains assigned to POP_A and others classified as For sampling dates see Table 1. putative hybrids, and were recovered on the right side of the POP_B-cluster, towards the POP_A-cluster.

Genetic diversity of the two populations
If all the strains analyzed by STRUCTURE are taken together, then a total of 103 distinct genotypes were found among the 164 strains assigned to POP_A (G/ N562.8%, Table 3), whereas 182 distinct genotypes were detected among the 246 strains assigned to POP_B (G/N574.0%). A total of 110 genotypes was identified (G/N595.7%) among the 115 strains that did not obtain an assignment probability $90% to POP_A or POP_B ( Table 3). We classified these strains as putative hybrids. The genotypic diversity among POP_A strains was comparable  Table 3). The genotypic diversity among POP_B strains was    Pseudo-nitzschia multistriata showed a marked seasonal cycle with high cell abundances restricted to a relatively short period. Results of the Wilcoxon test of each of the two populations sampled in consecutive years were negative, revealing no significant likelihood (p,0.05) for the existence of bottlenecks. Results were likewise after pooling all strains within single years.

Genetic structure of the F1 samples resulting from crosses
A Bayesian clustering analysis was carried out on the genotypes of the 525 strains from the field and the F1 strains generated from crossing experiments conducted in the study of Tesson et al. [11]. A cross between parental strains SY017 (POP_A) and SY278 (hybrid; cross 1 in Fig. 7, S2 Table) generated an F1 in which most strains were assigned to POP_A and a minority were classified as putative hybrids. A cross between parents SY017 (POP_A) and SY138 (POP_B; cross 2 in Fig. 7, S2 Table) generated an F1 in which most strains were classified as putative hybrids (i.e. 82.7% of F1). Crosses between parents SY138 (POP_B) and SY378 (POP_B) and parents SY138 (POP_B) and SY379 (POP_B; crosses 3 and 4 in Fig. 7, S2 Table) generated an F1 of which most strains were classified as POP_B and only a few as hybrids (i.e. 0 out of 25 F1 and 7 out of 62 F1 in crosses 3 and 4, respectively).

Discussion
Results of the Bayesian analysis and the FCA of the microsatellite genotypes obtained from the Pseudo-nitzschia multistriata strains uncovered two genetically distinct populations (denoted POP_A and POP_B) in the Gulf of Naples. A substantial proportion of strains could not be assigned with a $90%-probability to either of these populations and may represent hybrids between them. POP_A dominated in 2008, was largely replaced by POP_B in 2009 and 2010, and apparently returned to dominance in 2011. This latter observation is based on only a single sample with 13 strains, but all of these are assigned to POP_A. Samples dominated by strains assigned to one population usually included a minority of strains assigned to the other, as well as strains classified as putative hybrids, suggesting that the two populations persisted in sympatry throughout the four years of our study. Crosses between parents assigned to the two populations produced viable progeny with a genetic fingerprint by and large comparable to that of field strains classified by STRUCTURE as putative hybrids.

Genotypically distinct populations occur in sympatry
Genetically distinct populations in the same or overlapping areas, i.e., in sympatry or parapatry, have been reported before in planktonic species of diatoms (Pseudonitzschia pungens, Ditylum brightwellii, Skeletonema marinoi) [3,6,36,37,8] and dinoflagellates (Alexandrium fundyense and A. minutum) [9,10,38]. Our observation of two apparently sympatric populations among the Neapolitan strains of P. multistriata resembles results by Casteleyn et al. [3], who observed genotypically distinct populations within one of the three clades of Pseudonitzschia pungens co-occurring at their western European sample sites. Their sympatric populations exhibited an F ST value comparable to that between the two Neapolitan P. multistriata populations. Remarkably high genotypic difference was detected between populations of P. pungens along the Pacific Northwest (F ST 50.418; [6]), but these populations were attributed subsequently to different clades and are likely to belong to different cryptic species [36]. Similarly, four genetically distinct populations detected in Ditylum brightwellii [37] probably also represented closely related species [39]. In our case, results demonstrate that strains of the two P. multistriata populations can interbreed, and therefore, belong to a single biological species. When genetically distinct populations are detected within a species, environmental or biological factors can usually be identified that keep these populations apart. For instance, a population of the diatom Skeletonema marinoi in a Swedish fjord differed genetically from one in the nearby open sea (F ST ranging from 0.200 to 0.267) [8]. In this case, a shallow sill limits water exchange between the fjord and the open sea, maintaining environmental differentiation [8]. In addition, the species produces resting spores that can survive in the sediments for decades [18], thus permitting these populations to persist in their distribution areas [8,18]. Similarly, populations of the cyst-forming dinoflagellate Alexandrium fundyense from the open Gulf of Maine differed genetically from those in adjacent coastal ponds [10], although genetic differentiation was also recorded along the development of a single bloom [38]. In the case of P. multistriata, however, the two populations co-existed in a coastal area without physical barriers. Resting stages are not known for this species [40], nor do we have any evidence for gradual seasonal succession of POP_A and POP_B, thus excluding such explanations for the co-occurrence of genetically distinct populations in the Gulf of Naples.
Both P. multistriata populations appeared to persist over the entire growth season and throughout the four years of our study. Their relative proportions changed between consecutive years, but remained, by and large, similar across the growth season within years. The Gulf of Naples is an open system, exhibiting complex surface hydrodynamics and marked exchange with adjacent gulfs along the southwestern Italian coast and with the open Tyrrhenian Sea [41]. So, the detection of similar proportions between POP_A and POP_B throughout a growth season, or even between two consecutive growth seasons, leads us to conclude that these populations co-occur in the same proportions far beyond our sampling point.

Population structure may change in consecutive years
This study is the first in Pseudo-nitzschia -or in any phytoplankton species as far as we are aware of -revealing a turnover of populations between consecutive years, i.e., the dominant population becoming a minority and vice versa. Apparently, the turnover happens when the species is not observed in the plankton. Since benthic resting stages are unknown for Pseudo-nitzschia [40], the species must persist in the plankton outside the growth season either below the detection threshold for routine LM-observation or beyond the depth range at which we sample the plankton. Periods of extremely low densities could render populations prone to bottleneck effects. However, no such a bottleneck was evident from the results of the Wilcoxon tests. The genetic difference between POP_A in 2008 and POP_B in 2009 cannot be explained assuming a bottleneck since the alleles present in POP_B (2009 and 2010) do not represent a sub-set of the alleles present in POP_A (2008); POP_B exhibits several alleles absent in POP_A (2008). Therefore, other factors must be responsible for the turnover, possibly ones related to the life cycle or to mortality, affecting the populations differently.

Possible effects of the life cycle on population structure
The virtual disappearance of POP_A after 2008 could be a consequence of the episodic nature of sexual reproduction in P. multistriata and a differential success of sexual reproduction in the two populations. We do not have direct observations of the sexual phase of this species in the field. Yet, a model developed by D'Alelio et al. [21] from patterns of cell size distribution in field samples over time predicts that sexual reproduction occurs in the autumn of every second year. This peculiar timing is made possible because cell size in P. multistriata, as in most other diatom species, is a function of time due to a gradual reduction of average cell size (length) with on-going mitotic division. Once the size-range for sexual maturity has been reached (ca. 55 mm, after two years in P. multistriata), cells can reproduce sexually, giving rise to large F1 cells (82-72 mm) [22].
The model predicts that towards the end of the growth season of every odd calendar year, all cells in the populations belong to a single, two-years old cohort of short cells (55-40 mm) that have reached the size-range for sexual maturity and can give rise to an F1 cohort of very long, juvenile cells towards the end of the growth season. The following growth season (i.e., in an even calendar year) shows a bimodal cell size distribution with short cells (ca. 50-30 mm) belonging to a now three years old remainder of the parental cohort, and a juvenile, one-year old cohort of long cells (70-55 mm) that are not yet sexually mature. The latter will reach sexual maturity only at the end of the next growth season (i.e., an odd year; see figure 3 in [21]).
The model [21] was inferred from cell size distributions observed in the field until 2006, whereas our observations span the years 2008-2011. If the model is extrapolated to our sample period, then only one cohort is expected during the growth seasons of 2009 and 2011, with sex occurring towards the end of these seasons, and two cohorts are expected in the years 2008 and 2010 without sex. The distribution of cell lengths in the natural samples (grey dots in Fig. 3) was indeed bi-modal in 2008. In 2009, the distribution was narrow and unimodal, representing the now sexually mature, two years old cohort. In 2010, the distribution was continuous but very broad, representing in the lower part the short cells belonging to the remainder of the 2007 cohort and in the upper part the long cells of the juvenile 2009 cohort.
In 2008, the cells sampled for genetic analyses (triangles in Fig. 3) belonged basically all to the remainder of the cohort that underwent sex in 2007 and was destined to perish between the growth seasons of 2008 and 2009. The vast majority of this cohort belonged to POP_A (green triangles, Fig. 3). In 2009, the sampled cells belonged to a single, now sexually mature cohort, with the vast majority of cells being assigned to POP_B (red triangles, Fig. 3) or classified as putative hybrids (black triangles, Fig. 3). In 2010, POP_B and putative hybrids dominated again, both among the three-years old remainder of the 2007-cohort of short cells as well as among its offspring -the one-year old 2009-cohort of long cells.
The virtual disappearance of POP_A between 2008 and 2009 could be due to a low success of sexual reproduction within POP_A versus that within POP_B at the end of the growth season of 2007. If this is correct, then we expect the larger cells in 2008 to belong mainly to POP_B. Unfortunately, we did not sample this cohort well; only three cells sampled in 2008 were long. The dominance of POP_B in both 2009 and 2010 is in accordance with expectations, because in 2009 a sexually mature cohort dominated by POP_B and hybrid cells must have generated a new cohort of large cells detected in the following year, belonging mainly to POP_B or classified as putative hybrids.
Thus, different success of sexual reproduction of the two populations in 2007 and a comparable success in 2009 could explain the observed patterns until the end of 2010. However, this explanation requires that the two populations use different cues or triggers to commence sexual reproduction, or that they have a mating preference for partners belonging to the same population, for neither of which we have any evidence. Moreover, a different reproductive success cannot explain the dominance of POP_A in the only sample available to us in 2011 and neither does it explain the marked rise in the proportion of POP_A at the end of the 2010 growth season. Therefore other factors must be at work to explain these transitions.

Persistence of genetic structure in the face of hybridization
The persistence of POP_A and POP_B in large proportions over the four years is unexpected in the face of the large proportion of strains classified as putative hybrids and the apparent lack of reproductive barriers between POP_A and POP_B strains. At least under controlled laboratory conditions, cells of the opposite mating type undergo sexual reproduction and their F1 exhibit microsatellite genotypes largely in accordance with Mendelian inheritance rules [11]. Most F1 strains resulting from crossing POP_A and POP_B cells (cross 2) are not assigned to either one of these populations above the 90% probability threshold by STRUCTURE and are therefore classified as putative hybrids, suggesting that field strains not assigned to either one of these populations above the 90% probability threshold are hybrids in a biological sense as well. Such hybrids are fertile because a cross of a strain classified as hybrid with a POP_A-strain (cross 1) generated a perfectly viable F1, which in its turn is fertile as well (unpublished results). Results of the mating experiments in the lab show no evidence for mating barriers between the populations, though in these mating experiments partners were not offered a choice.
Our results indicate that hybridization is common between POP_A and POP_B also in the field. First, the results of the FCA reveal a continuum from POP_A via the assigned hybrids to POP_B, and even if the assigned hybrids are ignored, there exists marked gene flow between POP_A and POP_B. Second, many strains identified as putative hybrids by STRUCTURE shared alleles seen otherwise only among POP_A-strains with those observed otherwise exclusively among POP_Bstrains. Third, samples composed of strains of two reproductively isolated populations should reveal a Wahlund effect, i.e., a lower number of heterozygotes than expected if the samples were composed of a single population in Hardy-Weinberg equilibrium [42]. Instead, almost all samples show heterozygote overdominance. Such a pattern is typical for populations connected by substantial gene flow [42].
Interestingly, some F1-strains of cross 1 (hybrid6POP_A) were classified as putative hybrids and others as belonging to POP_A, indicating that a strain assigned to POP_A does not necessarily need to be derived from POP_A parents. Likewise, several F1-strains resulting from a cross between two POP_B parents were classified as putative hybrids. Two explanations can be given for these observations. First, assignment probability of a strain to, e.g., POP_B depends on the presence of private alleles for POP_B and on alleles found with higher frequency in POP_B strains. If an F1 strain happens to inherit from its POP_B parents alleles not very specific for POP_B, then its assignment probability to POP_B can fall below the 90% and hence, it is classified as a hybrid. Second, the precision of assignment probability improves with the number of loci, the number of individuals, and the F ST between the two populations. The smaller these values, the larger the imprecision around an assignment. In the case of our data, the total number of strains is high but the number of loci and the F ST are modest, suggesting that some of the assignments could be inaccurate. This argument does not explain why in the face of hybridization the majority of the strains remain assignable to POP_A or POP_B throughout the years of our sample campaign.
A reason why the two populations remain to be encountered in sympatry as genetically distinct entities in the face of hybridization might be that they both are relatively new arrivals. However, the species appeared in the Gulf of Naples in 1996 and must have gone since then through eight periods of sexual reproduction until 2011, according to the model by D'Alelio et al. [21]. If the two populations arrived simultaneously and remained confined in the Gulf of Naples ever since, then these eight phases of sexual reproduction sufficed to merge them entirely. However, this has not happened.

Towards a meta-population explanation
Periodical appearance and disappearance of populations have been observed in meta-population structures of organisms living in fragmented habitats, where single local populations crash and the site is re-populated by other, distinct, but not completely disjoint, populations [43]. Apparently, P. multistriata could follow such a meta-population-like structure, with distinct but connected populations blooming in different regions. This scenario could explain away the sudden rise in the proportion of POP_A strains sampled at the end of the 2010 growth season, or any other change in the proportions. We lack information about the population genetic structure of P. multistriata in other coastal regions in the Mediterranean Sea and in other basins, though the species is known to occur elsewhere along the Tyrrhenian coastline [44] and is, in fact, distributed globally (see references in [45]). Casteleyn et al. [3] demonstrated the existence of multiple genetically distinct, geographically distant populations in Pseudo-nitzschia pungens clade I, with samples taken at each of the geographically distant locations containing one or a few immigrants, or members of minor resident populations, exhibiting genotypes typical for populations dominating elsewhere. Our results show that within the PCA plot in Fig. 6, the 2011-sample dominated by POP_A strains was recovered distantly from the cluster of 2008 samples dominated by POP_A, indicating genetic changes over time, possibly resulting from exchange with yet unknown populations elsewhere. Extensive genotyping of P. multistriata samples from geographically distant places will reveal if such a patchwork of populations exists also in P. multistriata.
In conclusion, we hypothesize that P. multistriata population dynamics in the Gulf of Naples follows a meta-population-like model, which includes establishment of populations by cell inocula from previous populations at the beginning of each growth season as well as remixing and dispersal governed by water masses that host these populations. Further multiannual studies of population genetic diversity and structure will help clarifying the environmental and internal factors that govern the evolution and fate of the populations of unicellular microalgae.
Supporting Information S1 Figure. The number of clusters (K; populations) as estimated following Evanno et al. [29] using the web-based program Structure Harvester [30]. Independent runs were performed for K51 to 22; results for K.5 have been pruned from the figure. doi:10.1371/journal.pone.0114984.s001 (TIF) S1