Population Genetic Structure at the Northern Edge of the Distribution of Alexandrium catenella in the Patagonian Fjords and Its Expansion Along the Open Pacific Ocean Coast

In southern Chile, Alexandrium catenella is the main species generating HABs and over time it has expanded its range since it was first recorded in the Magallanes region in 1972. In 2016 a severe bloom of an Alexandrium species occurred, which was notable for its intensity and geographical extent, extending into new areas to the north of the Patagonian fjords including areas along the open Pacific Ocean coast. Given the exceptional nature of this event, we verified the taxonomic classification of the species that generated the bloom and evaluated the influence of the range expansion process on its genetic structure and population diversity. This was achieved by isolating clones collected in 2014 from cyst beds located at the northern limit of its then known distribution, along with clones isolated from the water samples taken during the 2016 bloom. These clones were characterized genetically with LSU rDNA and AFLPs molecular markers. Phylogenetic analyses showed that all clones were aggregated in the Group I of the A. tamarense species complex, which confirmed that A. catenella was the species responsible for the 2016 bloom. High genetic diversity was observed within populations though there were no significant differences between populations. Furthermore, genetic structure showed an isolation by distance relationship among populations, and several analyses consistently indicated a high divergence among the population groups derived from both cysts and vegetative cells. Despite this study not finding the patterns expected for species range expansion (i.e. diversity gradient and/or high population divergence), the genetic diversity and structure indicated that these were influenced by the geographic distance effect and the physical restrictions to gene flow, along with the demographic processes that occurred during the distinct phases of the life cycle of A. catenella.

In southern Chile, Alexandrium catenella is the main species generating harmful algal blooms (HABs) and over time it has expanded its range since it was first recorded in the Magallanes region in 1972. In 2016 a severe bloom of an Alexandrium species occurred, which was notable for its intensity and geographical extent, extending into new areas to the north of the Patagonian fjords including areas along the open Pacific Ocean coast. Given the exceptional nature of this event, we verified the taxonomic classification of the species that generated the bloom and evaluated the influence of the range expansion process on its genetic structure and population diversity. This was achieved by isolating clones collected in 2014 from cyst beds located at the northern limit of its then known distribution, along with clones isolated from the water samples taken during the 2016 bloom. These clones were characterized genetically with LSU rDNA and AFLPs molecular markers. Phylogenetic analyses showed that all clones were aggregated in the Group I of the A. tamarense species complex, which confirmed that A. catenella was the species responsible for the 2016 bloom. High genetic diversity was observed within populations though there were no significant differences between populations. Furthermore, genetic structure showed an isolation by distance relationship among populations, and several analyses consistently indicated a high divergence among the population groups derived from both cysts and vegetative cells. Despite this study not finding the patterns expected for species range expansion (i.e., diversity gradient and/or high population divergence), the genetic diversity and structure indicated that these were influenced by the geographic distance effect and the physical restrictions to gene flow, along with the demographic processes that occurred during the distinct phases of the life cycle of A. catenella.

INTRODUCTION
Around the world, factors such as climate change and eutrophication of coastal waters have favored the recurrence and severity of harmful algal blooms (HABs), increasing the probability that harmful species will expand their ranges to new areas (Wells et al., 2015). Species of the A. tamarense species complex are among the most prominent HAB species (Anderson et al., 2012b;Murray et al., 2015). These species, following a trend shown by other harmful algae (Anderson et al., 2012a), have over recent decades expanded their geographic ranges in both the northern and southern hemispheres (Penna et al., 2005;Persich et al., 2006).
Underlying the process of species range expansion are the colonization-extinction dynamics of the populations, whether a proportion of the migrating population will survive or not and successfully reproduce in the new area (e.g., Sexton et al., 2009). Theoretically, it is expected that the range expansion process will generate a gradient in genetic diversity and marked genetic differences between populations, due to a fraction of the migrants establishing in the new area changing due to the effects of genetic drift (Excoffier et al., 2009;Sexton et al., 2009). The genetic characterization of populations using AFLPs (Amplified fragment length polymorphisms) has permitted the study of genetic diversity and differentiation in several harmful species, with an important advantage over microsatellite markers due to the high amount of genome-wide loci observed (Rengefors et al., 2017) and its power in identifying weak genetic population differentiation (Alpermann et al., 2009). Although, AFLPs disadvantages include variation in the precision of sizing of fragments, and the assumption that bands of the same size are homologous, might not always be right, especially at higher levels of divergence (John et al., 2004;Fry et al., 2009). For instance, AFLP population characterizations have permitted the examination of the range expansion process of species such as Gonyostomum semen and Alexandrium ostenfeldii, which show no evidence of diversity gradients, but do show levels of population divergence and isolation by distance (Tahvanainen et al., 2012;Sassenhagen et al., 2015). However, in A. ostenfeldii, Brandenburg et al. (2018) recently described a contrasting genetic and phenotypic variability between a highly diverse population established in the Baltic Sea approximately at 3000 -8000 BP and the low variability of the populations that colonized the Netherlands localities in 2012, that represent examples of established and founder population respectively. Furthermore, genetic characterization using AFLPs and microsatellites of Alexandrium populations have shown the presence of high levels of genetic diversity, which is related in turn to high phenotypic variability, factors which are relevant to their adaptive potential (Alpermann et al., 2010;John et al., 2015;Kremp et al., 2016;Brandenburg et al., 2018). In the case of A. tamarense species complex, the population genetics characterization using microsatellites and AFLPs, have revealed that the population divergence could have been generated by demographic changes that have occurred over time, in the populations of cysts in the sediments, and/or in the population of vegetative cells in the water column, and as a result of reduced gene flow. For example, it has been hypothesized that the vegetative cells generated in a bloom would come from diverse banks of cysts, formed by different cohorts, that have accumulated in the sediment over time (Alpermann et al., 2009). In addition, during the bloom the environmental conditions could favor changes in allele frequency, due to selective pressures on certain genotypes, which could then be amplified by asexual reproduction (Alpermann et al., 2009;Dia et al., 2014). Finally, the potential migration of vegetative cells could be limited by oceanographic barriers (Casabianca et al., 2012) or restricted by geographic distance (Nagai et al., 2007), generating high levels of differentiation between populations.
The taxonomic identification of species comprising the A. tamarense species complex is fundamental not just for monitoring programs and implementing management strategies , but also, in order to verify that the genetic analyses of populations are made with conspecific individuals. Over time various authors have discussed the taxonomic classification of the species of the A. tamarense species complex which occurs in southern Chile, proposing a range of nomenclatural alternatives, such as A. tamarense (Jedlicki et al., 2012), A. fundyense Mcneill et al., 2014;Wang et al., 2014) or A. catenella (Cordova and Müller, 2002). However, this discussion was recently concluded with priority being given to the name A. catenella (Willem, 2017). This classification was based in large part on the homology of ribosomal genetic characters (rDNA) rather than morphology, as morphology can be highly plastic and variable within the species . For example, historically the capacity to form chains has been used as a character for distinguishing between the species of the A. tamarense species complex. However, evidence indicates that the formation of chains is influenced by environmental variables such as a high pH and low pCO 2 (Mardones et al., 2016b). On the other hand, the phylogenetic analysis with rDNA molecular markers consistently produces five clades, where clones are separated, including those with a sympatric distribution .
In southern Chile Alexandrium catenella is the principal species causing HABs, resulting in serious ecological, social and economic impacts throughout the Patagonian fjords, where this species has expanded its range northward (Varela et al., 2012). Since 1972, when it was first identified in the Magallanes region (ca. 56 • S), monitoring programs have registered the northward expansion of the species range distribution, arriving in the Aysén region in 1992, and in the waters of the southern end of Chiloé Island (ca. 43 • S), Los Lagos region, in 2002 (Molinet et al., 2003;Varela et al., 2012). Until 2015 the blooms of A. catenella were restricted to the fjords and channels between Magallanes and southern Chiloé Island (Molinet et al., 2003;Mardones et al., 2010;Varela et al., 2012). However, a severe bloom of A. catenella occurred during 2016 that extended from the central-northern region of Aysén (45 • 27 S) to the Chiloé inland sea, and for the first time, the bloom occurred in the open Pacific Ocean along the west coast of Chiloé Island (Figure 1; Buschmann et al., 2016).
The exceptional bloom of 2016 was surprising due to its wide geographical distribution, extending into new areas along the open coast of the Pacific Ocean and Chiloé inland sea ( Figure 1). The event chronology indicated that the first evidence of the bloom was recorded around November-December of 2015 in the central-northern area of the Aysén region (45 • 23 S). Later, during January-February of 2016 the highest abundances of vegetative cells were observed in the southern area of the Chiloé inland sea (Buschmann et al., 2016). Until then, the bloom had followed the same pattern of geographical distribution observed during the severe blooms that occurred in 2002 and 2009 (Molinet et al., 2003;Mardones et al., 2010). However, from the first week of March 2016 the bloom extend northward, through the Chiloé inland sea up to Desertores Island (42 • 41 S), from where high levels of paralytic shellfish poison (PSP) were detected in mollusks (Buschmann et al., 2016). By April, the bloom declined in the Chiloé inland sea, but surprisingly the highest concentrations of PSP in mollusks were found along the exposed Pacific coast of Chiloé Island (42 • 37 S). This was the first evidence of the spread of the bloom into the open Pacific Ocean. Subsequently, PSP increased rapidly in mollusks living further north along the exposed coast, vegetative cells abundance also increased in areas beyond the Patagonian fjords (up to 39 • 42 S). Finally, by mid-May, both the cellular abundance and PSP concentrations began to decline rapidly in most areas (Buschmann et al., 2016;Hernández et al., 2016).
Given the exceptional nature of the 2016 bloom, several hypotheses have been proposed to explain the origin and dynamics of the range expansion, including doubt as to whether the species that caused the bloom were the same species that have historically affected the Patagonian fjords. Thus, the objective of the current investigation was to verify the taxonomic classification of the species that generated this exceptional bloom and evaluate the influence of the range expansion process on its genetic structure and population diversity. In order to achieve this, rDNA large subunit (LSU) and AFLPs molecular markers were used for the genetic characterization of clones. This was done with two groups of samples; the first group were vegetative cells isolated from cysts collected in 2014 from sediments located at the northern limit of its known distribution; and the second group were vegetative cells isolated from samples taken from the bloom expansion areas observed during the 2016 bloom in waters of the open Pacific Ocean coast and from the waters of the Chiloé inland sea (Figure 1).

Sampling and Cultivation of Clones
The samples used for the genetic characterization of the exceptional bloom of 2016 were collected in April, during the peak of the event, from central-northern area of its distribution, and they were compared with samples collected previously, in 2014, from sites near the origin of the same bloom event.
In 2014 sediment samples containing A. catenella cysts were collected by SCUBA divers from the Aysén 1, Aysén 2 and Quellón localities (Populations from the northern limit of the then known species distribution (Figure 1 and Table 1). The sediment samples were transported to the laboratory, Centro i∼mar Universidad de Los Lagos, where the cysts were isolated using the methodology described by Varela et al. (2012). The cysts were then placed in multi-welled cultivation plates with L1 medium and exposed to conditions that facilitate germination (i.e., 12 • C, 30 psu, 35 ± 5 µmol photons m −2 s −1 , and a photoperiod of 16:8 h light/dark). From each well where cysts germinated, a single cell was selected and transferred to a new well (48-well plate) with the objective of generating clonal cultures. The single cell was isolated using an inverted microscope (Olympus CKX 42) by an extended Pasteur pipette (Anderson and Kawachi, 2005) and then deposited in a well with an area of 0.64 cm 2 (48-well plate) with 500 µl of L1. Subsequently, the well was carefully examined using an inverted microscope to verify the successful transfer. Additionally, in April of 2016, seawater samples were collected during the A. catenella bloom at Quemchi, Hueihue, Canal Chacao, and Bahía Mansa localities (Figure 1 and Table 1). These samples were part of an opportunistic sampling that attempted to determine the extent of the 2016 bloom beyond its historic distribution. Seawater samples were also transported to the laboratory (Centro i∼mar) where vegetative cells were isolated and maintained following the same methodology and environmental conditions described above. Moreover, after 24 days of culturing, cell characteristic were determined (N = 50 per locality), i.e., viability in culture (number of well with cells alive), cell density, growth rate and the number of the cells forming chains (c.f. Paredes et al., 2016). Finally, cell isolates, both from cysts and vegetative cells, once they began to multiply, were transferred first to culture tubes, then to 250 ml conical flasks, from which cells were harvested for DNA extraction. Thus, all the strains utilized in this study were established as clonal cultures and maintained under the same laboratory conditions (i.e., 12 • C, 30 psu, 35 ± 5 µmol photons m −2 s −1 , photoperiod of 16:8 h light/dark, and L1 cultivation medium).

Taxonomic Identification of the Clones
To verify the taxonomic classification of the clones a total of 96 clones isolated from the 2016 bloom and the cysts were characterized using rDNA LSU (D1D2) molecular marker as described by Scholin et al. (1994). DNA extractions were performed with the DNeasy plant mini Kit (Qiagen, Hilden, Germany) using 250 ml of culture harvested in the exponential phase of growth, and concentrated by centrifugation. The amplification of the DNA was performed using a thermocycler (Px2 Thermal Cycler, Thermo Electron Corporation, Waltham, MA, United States) with the following settings: one cycle of 5 min at 95 • C; 30 cycles of 45 s at 95 • C, 45 s at 56 • C and 1 min at 72 • C; and a final elongation for 10 min at 72 • C. The amplifications were verified using an agarose gel (1.4%), and then, the PCR products were sent to biotechnological company Macrogen R for sequencing. The obtained sequences were evaluated visually and the homology with the sequences available in GenBank was explored. Finally, a phylogenetic reconstruction was performed using the sequences generated by John et al. (2014) which were obtained from the Supplementary Material (i.e., Supplementary  Figure 1). Alignments of the SSU, ITS and LSU sequences for the phylogenetic analyses presented in this study. 1-s2.0-S1434461014001011-mmc1.doc). The reconstruction was performed with the maximum likelihood (ML) method, with a bootstrap analysis (1000 replicates) which evaluated the robustness of the nodes, and using an alignment of 826 base pair (bp). These analyses were performed with the software PAUP (Swofford, 2002). To perform the ML reconstruction, previous DNA sequences were analyzed with the MODELTEST 3.7 software which determined the most appropriate model of nucleotide substitution for the data (Possada and Crandall, 1998). Thus, the GTR + G + I model was selected with the following

Genotyping Using AFLP Analyses
The genetic characterization was performed with 150 clonal strains, using the AFLPs technique and 100 ng of DNA, as suggested by Vos et al. (1995) and Mardones et al. (2016a). DNA extractions were performed with the DNeasy Plant Mini Kit (Qiagen), using 250 ml of culture harvested in the exponential phase of growth and concentrated by centrifugation. The concentration and quality of the DNA (proportion 280/260 nm) were estimated using a NanoQuant Plate (Infinite M200 PRO). The primers EcoRI-5 GACTGCGTACCAATTCXXX 3 and MseI-5 GATGAGTCCTGAGTAAXXX 3 were used for the selective amplifications with the following combinations of primers marked with Fam dyes: EcoRI AAG x MseI CTA , EcoRI AAG x MseI CTT , EcoRI ACC x MseI CTA and, EcoRI ACC x MseI CTT . Selective amplifications were performed with the settings "touchdown" which included: 13 initial cycles with 30 s at 95 • C, 45 s at 56 • C (gradually reduced by 0.7 • C per cycle), and 30 s at 72 • C. This was followed by 25 cycles of 1 min at 95 • C; 1 min at 56 • C; 1 min at 72 • C, and finally an elongation of 10 min at 72 • C. The amplifications were verified in an agarose gel (1.4%) and then these were analyzed using an ABI PRISM 3500 sequencer (Applied Biosystems). Consequently, the AFLP fragments were visualized and edited with the software GeneMarker. In order to create a presence/absence matrix of comigrating AFLPs fragment, the sequence profiles were examined visually considering a fragment detection range between 100 and 15000 relative fluorescence units (RFUs) and between 50 and 400 bp. In addition, the reproducibility of the amplifications was determined visually by repeating 20% of the samples selected at random.

Intraspecific Distance, Genetic Diversity, and Population Structure
The intraspecific distance was determined using Nei's genetic distance (Saitou and Nei, 1987) while the genetic diversity of the populations (H) was calculated according to Nei (1973). The genetic diversity differentiation between populations was evaluated by an ANOVA, performed in a generalized linear mixed model (GLMM) framework. The genetic diversity for each locus was estimated, and the populations were considered to be a random effect to account for the lack of independence between the samples. In addition, a negative binomial residuals model distribution was used, and the null hypothesis (H 0 ) was evaluated using a likelihood ratio test with a significance level (α) of 0.05 (Venables and Ripley, 2002). The estimation of genetic diversity per locus and population was performed with the ARLEQUIN 3.5.1.2 software (Excoffier et al., 2005) while the GLMM was performed using R software (R Core Team, 2017) and the package "lme4" (Bates et al., 2015).
The genetic structure was estimated using different analyses. Firstly, analysis of molecular variance (AMOVA) and the significance of the F ST , F SC, and F CT were tested using a nonparametric analysis with 10,000 permutations. For this analysis the following sources of variation were considered: (1) between groups, composed of isolated populations in the northern limit of distribution (derived from cysts) and the populations isolated from the bloom; (2) variation between populations within groups, and (3) variation within the population. Furthermore, the level of differentiation between populations was estimated using the fixation index F ST . Both analyses were performed using the software ARLEQUIN 3.5.1.2 (Excoffier et al., 2005). Subsequently, in order to visualize the similarities/dissimilarities between populations, a Principal Component Analysis (PCA) was performed using a matrix of genetic distances between populations which was estimated using Nei (1978) method. This was done using the software GENALEX 6.5 with a binary model for haploid organisms (Peakall and Smouse, 2012). The identification of genetic populations and the level of admixture of each genotype, or proportion of genotype assigned to each population, were estimated using a Bayesian cluster analysis implemented in the software STRUCTURE 2.3.4. The simulation was performed using a correlated allele frequency model and an ancestry mixed model without using location as prior information (Falush et al., 2003). The simulations were conducted using values of K populations from 1 to 10; a burning length of 100,000 and 100,000 repetitions of Monte Carlo Markov Chains (MCMC), and the results were evaluated using likelihood of K populations [L(K)]. The evaluation of K populations and the visualization of the genotype assignments were made using the web application www.pophelper.com (Francis, 2016). Finally, the software ARLEQUIN was used to conduct a Mantel analysis to evaluate isolation by distance, which included the paired values of F ST (Y matrix) and the geographic distance in Km (X matrix).

Disequilibrium of the Multilocus Linkage
The index of multilocus disequilibrium − r d (Agapow and Burt, 2001) was calculated for each population and significance was determined using 10,000 permutations. The values of the index vary between 0 and 1 indicating equilibrium and disequilibrium respectively. The analyses were performed with R software (R Core Team, 2017) and the package "Poppr" (Kamvar et al., 2014).

Taxonomic Identification
Only 68 sequences of molecular marker rDNA LSU (D1D2) were used ( Table 1). The remaining sequences were eliminated after visual evaluations of the chromatograms which suggested ambiguity in the assignment of some nucleotides i.e., no clear nucleotide peak or nucleotide peak overlaying. The phylogenetic analysis indicated that sequences of the clones isolated from both the 2016 bloom and the northern limit formed part of the same clade, corresponding to the Group I of the A. tamarense species complex (Lilly et al., 2007;John et al., 2014) (Figure 2). Moreover, Chi-square test did not show  differences in the base frequencies among sequences used for the phylogenetic reconstruction (χ 2 = 62.243, p > 0.05). The maximum genetic distance, according to maximum-likelihood distance, was observed between A._pacificum_ver v/s A._leei_ver strains (0.924). Meanwhile, among clones isolated from southern Chile the distance ranged between zero, for several pairwise clones isolated from closer and separate populations, and 0.069 between clones isolated from Bahía Mansa v/s Canal Chacao. Furthermore, in the clones isolated from the 2016 bloom, the number of cells in chains ranged from 2 to 8, cell viability from 0.396 to 0.875 probability, and growth rate from 0.140 to 0.331 div/day (c.f. Paredes et al., 2016). These responses were concordant with the characteristics of A. catenella strains previously isolated from southern Chile.

Intraspecific Genetic Distance and Population Diversity
Based on the evaluation of the banding patterns, following the RFU criteria and amplification reproducibility, 264 loci were used for the genetic evaluation, from a total of 108 clones distributed across seven populations ( Table 1). Similar genotypes were not observed between clones. The genetic distance between clones varied between 0.100 and 0.708, with smallest distances between Hueihue and Canal Chacao; and the greatest distances observed between the clones sampled in Quemchi and Hueihue. On the other hand, the mean genetic diversity (H) had values that varied between 0.289 and 0.415, for Aysén 2 and Bahía Mansa populations respectively (Figure 3). Although, significant differences in the genetic diversity (p > 0.05, GLMM) were not observed among populations.

Genetic Structure
The genetic variability of A. catenella was significantly structured (AMOVA, p < 0.05, for each of the fixation indices; Table 2). Although most of the genetic variability was explained by variation within populations (90.450%), the genetic differences were significant between groups, explaining 6.730% of the genetic variation, and between populations within groups, which explained 2.820% of the variation ( Table 2). The paired index of fixation F ST was significant for the majority of comparisons between populations, with the highest values being observed for the comparisons made between the populations derived from cysts and those derived from vegetative cells collected in the 2016 bloom (Table 3). Of the comparisons, the population in Quemchi had the greatest differences compared with each of the three populations derived from cysts (F ST = 0.120 -0.122), these were followed by Bahía Mansa (F ST = 0.104 -0.109), Hueihue (F ST = 0.101 -0.103) and the population in the Canal Chacao (F ST = 0.076 -0.093). Among the populations derived from cysts the genetic differences were relatively similar, with the largest value observed between populations of Quellón and Aysén 2 (F ST = 0.032). The paired values of F ST between populations derived from the bloom were more heterogeneous. There were no significant differences in the values of F ST between Bahía Mansa and the populations of the Canal Chacao and Hueihue, the same was true for the comparison between the Canal Chacao and Quemchi. However, there were significant differences in the values of F ST for the comparisons between Quemchi and Hueihue (F ST = 0.048), Quemchi and Bahía Mansa (F ST = 0.054), and between Hueihue and Canal Chacao (F ST = 0.031). These differences were relatively higher than those found for the comparisons between populations derived from cysts. PCA indicated segregation of the populations in at least three groups. Thus, the first coordinate explained 65.76% of the variation and showed evident genetic similarities between the populations isolated from cysts beds in 2014 (Aysén 1, Aysén 2 and Quellon) and populations isolated from 2016 bloom. Moreover, among these last populations the second coordinates, which explained 18.35% of the variation, showed genetic similarities between Bahía Mansa and Hueihue; and between Canal Chacao and Quemchi (Figure 4).
A similar pattern described with PCA was observed when K = 3 (Supplementary Figure 1) or K = 5 Bayesian cluster analysis was plotted, distinguishing clearly between the strains from cysts and those from the bloom population, and recognizing higher variability within the group of isolated from the bloom. However, the Bayesian cluster analyses suggested a K = 5 as the more probable number of clusters [L(K) = −13490, Supplementary Figure 2]. The strains derived from cysts, in the main, were assigned more or less in the same proportions in two subgroups (pink and light blue, Figure 5). In the case of the strains derived from the bloom, despite the higher degree of mixing in their assignment, it was also possible to distinguish between geographic populations. Thus, the majority of the strains from Quemchi and Canal Chacao presented a high proportion of assignment to a third subgroup (yellow, Figure 5). Meanwhile, the strains from Hueihue and Bahía Mansa presented a higher degree of heterogeneity, with individuals assigned to different subgroups (blue, red or pink, Figure 5) and others assigned more or less equally to the same subpopulations.

Isolation by Distance
There was a significant positive relationship between the genetic and geographic distances based on the results of the Mantel

Disequilibrium of the Multilocus Linkage
The index of multilocus disequilibrium − r d gave significant values (p < 0.05) for all the populations ( Table 4)

Phylogeny and Genetic Diversity
The phylogenetic reconstruction based on rDNA LSU sequences indicated that all the clones analyzed from the populations associated with the 2016 bloom, as well as those isolated from the northern limit of distribution, were aggregated in the Group I of the A. tamarense species complex, verifying that the species being studied was A. catenella (Lilly et al., 2007;John et al., 2014;Willem, 2017). The current taxonomic revision of A. tamarense species complex, supported more by ribosomal phylogeny than morphological characters , has merged into A. catenella species group some toxic strains previously referred to as A. fundyense or A. tamarense Willem, 2017). As a result, the studies made by John et al. (2004), Nagai et al. (2007, Alpermann et al. (2009Alpermann et al. ( , 2010, Erdner et al. (2011), andRichlen et al. (2012) with the toxic strains of A. tamarense and A. fundyense now correspond to A. catenella, with A. tamarense and A. fundyense considered as morphotypes (c.f. John et al., 2014). In turn, the A. tamarense species is the only known nontoxic species among the species comprising the A. tamarense species complex, which segregate in Group 3 . Similar to our results, clones isolated from the Patagonian fjords between 1994 and 2016 have been placed in the same Group I (Cordova and Müller, 2002;Aguilera-Belmonte et al., 2011;Varela et al., 2012;Mardones et al., 2016a). Moreover, several studies have consistently found PSP toxins in the clones characterized and in the shellfish after blooms where A. catenella was the trigger (Krock et al., 2007;Molinet et al., 2010;Aguilera-Belmonte et al., 2011;Varela et al., 2012;Buschmann et al., 2016). These facts reduce the possibility that other relevant species comprising the A. tamarense species complex are present in the Patagonian fjords, as has been suggested by other authors (e.g., Jedlicki et al., 2012), and they confirm that the diversity and variation in genetic populations observed in the present study are intraspecific. However, taking account the huge geographical extent of the 2016 bloom and the impossibility of sampling many localities in different times, we cannot discount the presence of the A. ostenfeldii. This species has been identified before in the Patagonian fjords in low densities in the same zones where the 2016 bloom initiated, and on some localities within the northward expansion of the bloom (e.g., Salgado et al., 2015;Pizarro et al., 2018). In this study, the high genetic diversity observed among A. catenella populations, isolated from both cysts and vegetative cells, are consistent with the high genetic diversity reported in A. catenella populations in the northern hemisphere (Alpermann et al., 2010;Erdner et al., 2011). Mutation is the fundamental mechanism generating genetic variation in both the asexual and sexual phase of all organisms (Rengefors et al., 2017).
In Alexandrium species the mutation can be produced during asexual division of vegetative cells in the water column, and also during genetic recombination of distinct cells (i.e., heterothallic mating; Mardones et al., 2016a) that generate cysts that sink to the seafloor (Brandenburg et al., 2018). In the growth season, these cysts hatch continuously and thereby supply new genotypes to the water column, a mechanism by which they contribute to the genetic diversity (Brandenburg et al., 2018). Possibly, the high genetic diversity observed within the A. catenella populations is an attribute that confers to the species the capacity to adapt to the variable and heterogeneous environmental conditions observed across the Patagonian fjords. There, latitudinal gradients in the physical-chemical variables explain the presence of ecoregion, bioclimatic region, and biogeographic patterns (Camus, 2001;Spalding et al., 2007). However, further studies are needed to understand if the genetic diversity of A. catenella is coupled with phenotypic variability, as some authors have observed in the northern hemisphere (e.g., Alpermann et al., 2010;Brandenburg et al., 2018).
The genetic diversity observed between the populations did not indicate a process of recent colonization. Even though the reduction in diversity and segregation are expected results for a species during range expansion (Excoffier et al., 2009;Sexton et al., 2009), in this study the high diversity detected using AFLP does not differ significantly among populations (p < 0.05, GLMM) and it was distributed mainly within them (90.450%, AMOVA). Similarly, in some HAB species that have expanded their distribution, gradients of genetic diversity between established and founder populations have not been found (Tahvanainen et al., 2012;Lebret et al., 2013). Indeed, the genetic diversity of the A. catenella populations off the coasts of Japan, Korea, and the United States of America, estimated from samples taken in different months, years, localities, and life cycle phases, consistently show high values, even when population genetic divergence was found (Nagai et al., 2007;Erdner et al., 2011;Richlen et al., 2012). Although, when genetic differences in diversity and/or variability levels have been observed, they have occurred among populations separated by large geographic distances (Masseret et al., 2009) and when the colonization events have been separated by considerable periods of time (Brandenburg et al., 2018). In these latter cases, the gene flow limitation (Masseret et al., 2009) and the extended time to population diversification (Brandenburg et al., 2018) may have favored the maintenance of the diversity/variability gradient among populations. In a rapidly growing organism with a short generation time, it has been suggested that the huge population sizes could be considered as virtually infinite and genotypic diversity could thus be maintained during the bloom events, even if the organism primarily reproduces clonally (Dia et al., 2014). Thus, the lack of a genetic diversity gradient among A. catenella populations may be related to the potential observed in the genus Alexandrium for maintaining high genetic diversity even in a bloom, and the considerable gene flow among the populations, as is suggesting by the low and moderate pairwise F ST values.

Population Structure
The observations that described the progressive expansion of A catenella indicate that the presence of this species in the Aysén region dates from 1992, but it was only in 2002 when an intense bloom that started in this region, reached further north up to Quellon on the south coast of Chiloé Island (Molinet et al., 2003). Considering the local oceanographic currents between the populations of Aysén and those around Quellón the dispersion among populations should be restricted (e.g., Sievers and Silva, 2003). Two strong surface (0-30 m) counter currents meet in the Golfo Corvocado and from here flow out through the Boca Guafo toward the Pacific Ocean (Sievers and Silva, 2003): one current flows south from Chiloé inland sea and, the other, flows north from the Canal Moraleda (Figure 1). Taking into account these currents and that the abundances of vegetative cells are mainly found in the surface waters, down to about 15 m (Molinet et al., 2003), it is thought that the gene flow between populations from both sides of the Boca Guafo should be limited, generating high levels of genetic differentiation. However, the lowest levels of differentiation (F ST = 0.020-0.032) and similar patterns of genotype assignment to genetic population (Bayesian aggregation analysis) were observed between these populations, indicating that the current systems are not a barrier to gene flow. Indeed, the meteorological conditions that prevail during spring and summer, with the predominance of southerly winds, could be enough to push the vegetative cells northward, particularly during intense bloom events, such as those that occurred in 1998, 2002(Molinet et al., 2003Mardones et al., 2010). The descriptions of these events indicated that they started in the central-northern areas of the Aysén region and then expanded to the north (Molinet et al., 2003;Mardones et al., 2010). Thus, successive bloom events could have favored the dispersal of the vegetative cells from Aysén toward the north, increasing gene flow and obscuring the evidence of an initial founder effect between the populations. AMOVA, PCA, and Bayesian cluster analysis (K3 genetic population) indicated that the segregation of the genetic variance was important at least in two levels: between the populations coming from both the previous northern limit of the distribution and from the bloom of 2016, and also between the populations inside of each of these two groups. The differences between the former two groups of populations explained 6.730% of the variance in the AMOVA or 65.76% of the variance explained by the first axes (PCA), it also accounted for the greatest genetic differences (F ST = 0.076-0.122). Within each of the two groups, smaller differences were observed between populations collected during the 2016 bloom (F ST = 0.031-0.054) and the smallest differences were observed between the populations derived from cysts (F ST = 0.020-0.032). In A. catenella populations significant values of F ST have been observed in different geographical contexts: higher values comparing different geographical regions (F ST : 0.281-0.308; John et al., 2004); or from moderate to high when comparing isolated subpopulations within the same region (F ST : 0.07-0.36; Alpermann et al., 2009). Thus, the differentiation observed both between populations in the north of the distribution and those from the 2016 bloom, and within each of these groups, could be considered as subpopulations rather than clearly divergent populations.
The significant correlation between genetic and geographic distances of the populations (p < 0.05, Mantel Test) seems to support the assertion that differences observed between the populations from the northern limit of the distribution and those from the bloom of 2016 can be attributed, in part, to restrictions of the gene flow due to the limited dispersion of vegetative cells through the marine currents. Similar observations have been made for A. catenella in Japan, where the genetic structuring was correlated with geographic distance (Nagai et al., 2007), and also for A. minutum, where differentiation among samples from different locations in the Mediterranean Sea were attributed to the complex regional geography and hydrodynamics (Casabianca et al., 2012). Along over 300 km of exposed coast over which the great bloom of 2016 occurred (Buschmann et al., 2016), only some of the vegetative cells could have been advected by the unusual southern winds that were recorded throughout the month of March (Buschmann et al., 2016), limiting the dispersion to the other fraction of the population, and as a consequence, generating the observed differentiation.
Additionally, the differentiation and the changes in the genotype proportions assigned to ancestral populations, observed in the Bayesian aggregation analysis (K5 genetic population), between populations derived from cysts and those isolated from the bloom, could be the result of distinct population dynamics at each stage in the life cycle from which the samples came. The life cycle of A. catenella is characterized by alternating between a planktonic population of haploid cells that reproduce asexually and a benthic population of resting cysts that results from the occasional sexual reproduction of planktonic populations (Anderson et al., 2012b). The heterothallic behavior of the gametes (Mardones et al., 2016a), and the implicit recombination during sexual reproduction (Anderson et al., 2012b), along with the capacity of cysts to remain viable in the sediments over long periods of time (years in some cases), means that banks of cysts in the sediments act as reservoirs of genetic diversity (Alpermann et al., 2009;Anderson et al., 2012b). This could explain, in part, the relatively homogeneous assignment of the genotypes within the subpopulations derived from cysts, containing individuals with allelic proportions of the five clusters from detected ancestral subpopulations. This level of genotypic admixture and the dominance of two clusters could be the result of the allelic recombination of ancestral populations, probably of different geographic origins. However, in these subpopulations there is no genetic evidence of successive cohorts of temporally differentiated cysts, as has been observed by Alpermann et al. (2009). Despite the few studies undertaken in this area, low concentrations of cysts are observed in the Patagonian fjords (Mardones et al., 2016a). Furthermore, rapid reductions in cysts density have been observed (Diaz et al., 2014), probably due to the relatively short period of mandatory dominance (Mardones et al., 2016a) and the germination potential over the year (Varela, personal observation). This all indicates that the banks of cysts are a dynamic reservoir, without the accumulation of many cohorts over time.
Within the populations from the 2016 bloom, the heterogeneity of assignment proportions of the genotypes in genetic subpopulations (Bayesian aggregation analysis) and the high levels of the linkage disequilibrium index provided evidence of the incidence of demographic factors on the population genetic divergence. The Bayesian aggregation analysis revealed a relative dominance of individuals with a greater proportion of alleles from one of the ancestral subpopulations, with the insertion of some individuals that could be considered migrants between populations. At one extreme of the differentiation was the population from Quemchi, with a greater number of individuals with a dominant allelic proportion (yellow cluster), than the other populations which appear to be in a state of transition. This apparent process of rapid genetic spatial differentiation over the development of the bloom has also been observed in A. minutum (Dia et al., 2014) and in A. catenella (Erdner et al., 2011); differentiation which could also be temporal (Richlen et al., 2012;Dia et al., 2014). In these cases, it has been argued that these rapid genetic changes may reflect selective effects that are nevertheless not strong enough to reduce the genetic diversity (Dia et al., 2014). On the other hand, if encystment occurs during the bloom process (Brosnahan et al., 2016), and it is induced by adverse environmental conditions, this could magnify the selective effects on the proportions of different genotypes, reducing or eliminating lineages which are less successful under the prevailing conditions (Erdner et al., 2011). Both combined effects, the rate of reproduction and the differential encystment, could potentially reduce gene flow between subpopulations, generating differentiation between populations whilst maintaining the diversity of the regional population (Erdner et al., 2011).

CONCLUSION
In this study we did not observe the patterns generated by a range expansion process i.e., diversity gradients and/or high inter-population divergence. However, the genetic diversity and structure indicated that these were influenced by the geographic distance effects and the physical restrictions to gene flow, along with the demographic processes that occur during the distinct phases of the life cycle. Thus, the characteristics of asexual and sexual reproduction (e.g., heterothallism), the participation of ancestral populations in the allele proportions of individuals, as well as the selective value of some clonal lines amplified by asexual reproduction or by differential encystment, could be account of the genetic diversity and genetic structure observed in A. catenella populations in southern Chile.

AUTHOR CONTRIBUTIONS
JP was responsible for the laboratory and statistical analyses, and generation of the manuscript. DV was responsible for the generation of the manuscript. JP, CM, AZ, KC, AV, and BO were responsible for the sampling collection, isolation of strains, and clones maintaining.

ACKNOWLEDGMENTS
Special thanks to Dr. Federico Winkler and reviewer for significant comments to manuscript improvement, and to Dr. Matthiew Lee for English review. Moreover, thanks to Ma. Andrea Pérez Ríos from CTI Team (Competencias Transversales en Inglés) and Convenio Marco FID 1758 (Formación de Profesores 2017-2020, Universidad de Los Lagos) for the written English course imparted, which helped to substantially to the manuscript redaction.