Genetic population structure of the toxigenic dinoflagellate Alexandrium catenella in the Patagonian Fjords System, southern Chile

Alexandrium catenella is the main species that form harmful algae blooms (HABs) in southern Chile. Since its first record in 1972 in the Magallanes region this species apparently has increased its range distribution from south to north. In this study, we investigate the influence of the range expansions distribution on the A. catenella populations genetic diversity and structure. This was achieved by isolating 33 clones from different localities along the Magallanes and Aysén region which were genetically characterized with Amplified Fragment Length Polymorphism (AFLPs) molecular markers. Results showed a latitudinal genetic diversity gradient from the south to north populations. Inter-populations genetic divergences were low but significant between both geographically close and distant populations. Results indicated that the genetic diversity differentiation could be generated by a founder effect, which is expected in populations that have expanded their distribution range. On the other hand, low levels of genetic divergences between distant populations seems point out that high gene flow occurs along coast of the Pacific Ocean, but also, seems hints the connectivity route between the Magallanes and Aysén populations i.e., the vegetative cells dispersion among populations, would occur through coastal Pacific coast and Boca del Guafo. Finally, the high values of multilocus linkage disequilibrium found between closer population of Aysén indicates that divergence could be influenced along with the reproductive dynamic of the vegetative cells. Key word: Alexandrium catenella, Chilean Patagonia fjords system, genetic diversity, genetic structure, dispersion


Introduction
In southern Chile, Alexandrium catenella (Whedon and Kofoid) Balech is the main species forming harmful algae blooms (HABs), affecting ecosystems and human health (Varela et al. 2012, Paredes et al. 2019. This species has been monitored for the last four decades throughout the Patagonian fjords (Guzmán et al. 2002, Seguel et al. 2005, Alves-de-Souza et al. 2008, Mardones et al. 2010, showing an apparent northward geographic expansion (Guzmán et al. 2002, Molinet et al. 2003, Mardones et al. 2010, Varela et al. 2012). The first report of A. catenella vegetative cells was in 1972, at the southern end of Chile (ca., 53°S), Magallanes administrative region (Muñoz & Avaria 1997). Since then, A. catenella progressively increased its distribution range, reaching the northern boundary of the Magallanes region (ca., 50°S) in 1981 Vol. 55, N°1, 2020 Revista de Biología Marina y Oceanografía favor changes in vegetative cells allele frequency due to selective pressures on certain genotypes. 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.
In the southern Chile, the genetic characterization of A. catenella clones (AFLPs; LSU-rDNA and rDNA-ITS) have confirmed that they belong to the group I of the A. tamarense species complex, and have shown haploid variability and diversity (Aguilera-Belmonte et al. 2011, Varela et al. 2012, Mardones et al. 2016, Paredes et al. 2019. Moreover, Paredes et al. (2019) studied recently A. catenella population genetic structure and did not observed genetic diversity differentiation among populations of the Aysén, Quellón, and populations isolated from the 2016 bloom. The aim of this study was to analyze the influence of the geographic expansion process on the A. catenella genetic structure and diversity. To achieve this goal, we used AFLPs molecular markers for the genetic characterization of clones isolated from different locations across the Magallanes and Aysén regions.

Sampling and clone maintenance
In 2009, water samples were taken in 8 localities along the southern end of Magallanes region to the northern of the Aysén region (Fig. 1). From these samples, Alexandrium catenella single cells were isolated to establish clonal cultures. Thus, 12 strains were isolated from the southern boundary of the Magallanes region (south Magallanes population; 54.145°S, 70.686°W), 8 from the northern boundary of the Magallanes region (north Magallanes population; 49.521°S, 74.348°W), and 13 strains from the northern boundary of the Aysén region. The strains isolated from the northern Aysén were coming from close locations (Fig. 1A). Thus, 3 strains were isolated from Ester (45.099°S, 73.415°W), 3 strains were isolated from Canal Davis (44.465°S, 73.844°W), 2 strains were isolated from Valverde (44.322°S, 73.79°W), and 5 strains were isolated from Ceres (43.974°S, 73.782°W). These strains were established as clonal cultures and maintained under standard environmental conditions (i.e., L1 medium, 12 °C, 30 of salinity, 30-40 µmol photon m -2 s -1 and 16:8h lightdark cycles). In addition, the taxonomic status of the clones used in this study was supported by the results of previous researches which had evaluated both morphological and genetic taxonomic classification (c.f., Varela et al. 2012, Mardones et al. 2016. These works verified that the clones belong to the group I, and as discussed by Paredes et al. (2019), the proper name of the species comprising the A. tamarense species complex in southern Chile is A. catenella. (Guzmán et al. 2002). In 1992, an extensive monitoring revealed the presence of A. catenella in the Aysén region (ca., 45°S), but the first bloom event in this area was registered in 1996 (Guzmán et al. 2002). Between 1996 and2002, at least four blooms were recorded in the Aysén region which favored its northward expansion to the southern part of Los Lagos region (ca., 42.16°S). Since then, the southern part of the Los Lagos region has experienced frequent bloom events (Guzmán et al. 2002, Molinet et al. 2003, Varela et al. 2012. In 2016, a severe bloom of A. catenella was widespread from northern Aysén region to the Inner Sea of Chiloé (ca., 41.8°S), and surprisingly this bloom spread for the first time through open waters along the Pacific coast, extending the distribution to ca., 39.7°S (Paredes et al. 2019). Similarly, an apparent expansion in the distribution range of Alexandrium species has been described along the Atlantic coast of South America (Persich et al. 2006) and in some zones of the northern hemisphere (Vila et al. 2001, Penna et al. 2005, Masseret et al. 2009, Anderson et al. 2012. It is predictable that as a species expands its geographic range it could experience the founder effect, a process that can affect the population genetic diversity and promotes genetic differentiation between populations (Excoffier et al. 2009, Sexton et al. 2009). This is generated by the establishment of only a fraction of the migrants in the new area during an expansion event, moreover they may also suffer the effects of genetic drift (Excoffier et al. 2009, Sexton et al. 2009). Thus, if A. catenella has expanded its distribution through the Patagonian fjords system a northward genetic diversity gradient and population genetic divergence are expected.
Using molecular markers several studies investigated the influence of gene flow and genetic drift in the population genetic diversity and structure of HAB species, at geographical and temporal scales (e.g., Nagai et al. 2007, Erdner et al. 2011. This approximation has been useful to know the potential dispersion of HAB species, the range expansions and colonization process (Tahvanainen et al. 2012, Lebret et al. 2013, intra-population structures (Alpermann et al. 2009), and populations dynamics within and between blooms (Lebret et al. 2012). In the case of A. tamarense species complex, the genetics characterization of populations using microsatellites and AFLPs, revealed that the population divergence may be generated by demographic changes that occurred over time, in the populations of cysts in the sediments, and/or in the population of vegetative cells in the water column, as a result of reduced gene flow (Nagai et al. 2007, Alpermann et al. 2009, Paredes et al. 2019. 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, accumulated in the sediment over time (Alpermann et al. 2009). In addition, the environmental conditions during the bloom could

DNA extraction and genetic characterization
DNA extraction was performed with the phenol-chloroform reaction following the Sambrook et al. (1989) protocol with modifications. Cells were harvested in exponential phases and concentrated by centrifugation (15 min at 5500 rpm in a Universal 32 R Hettich Zentrifugen). DNA extracted was measured with a spectrophotometer (NanoDrop 1000 Spectrophotometer, Thermo scientific) and its purity was determined using both 260/280 ratios and visual evaluations through an agarose gel (1.5%). Clones genetic characterization was performed with amplified fragment length polymorphism (AFLPs) techniques following Vos et al. (1995) and John et al. (2004) procedures, and using 1,000 ng extracted DNA. The EcoRI-5´ GACTGCGTACCAATTCXXX 3´ and MseI-5´ GATGAGTCCTGAGTAAXXX 3´ primers were used for selective amplification. The following four primer combinations, labeled with Fam dyes were chosen: EcoRIAAG x MseICTA, EcoRIAAG x MseICTT, EcoRIACC x MseICTA and EcoRIACC x MseICTT. The Touchdown program was used for DNA amplification in the thermocycler (P x 2 Thermal Cycler, Thermo Electron Corporation, MA, USA) which included 13 initial cycles with 30 s at 95 °C; 45 s at 56 °C (gradually reduced by 0.7 °C every second 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 10 min of elongation at 72 °C. Amplifications were verified in the agarose gel (1.5%) and then analyzed in an ABI PRISM 3100 automatic sequencer (16 capillaries Applied Biosystems). The results were edited with the Peak Scanner 1.0 (Applied Biosystems) software. For this end, fragments between 50 and 1,000 bp, and fluorescence peaks higher than 51 relative units were considered. Moreover, in order to minimize the error of the loci score all the clones were analyzed twice (PCR amplification and fragment analysis). Finally, a presence/absence matrix of co-migrating AFLPs fragments was generated.

Figure 1. A) Localities where clones of Aysén region were isolated (Ceres= north Aysén population; Valverde, Canal Davis and Ester= south Aysén population). Grey circles indicate the locations where the samples were collected for the isolation of Alexandrium catenella clones. B) Map of southern
Chile outlining the Patagonian fjords system / A) Localidades donde se aislaron clones de la Región de Aysén (Ceres= población del norte de Aysén; Valverde, Canal Davis y Ester= población del sur de Aysén). Los círculos grises indican las localidades donde fueron aislados los clones de Alexandrium catenella. B) Mapa del sur de Chile que describe el sistema de fiordos Patagónicos

Intraspecific genetic distance
Genetic distance among A. catenella clonal strains was estimated using Nei's (Nei 1978) corrected genetic distance and were performed using the R software (R Core Team 2017) with the "APE" package.

Population genetic diversity and structure
In order to aggregate the few clones isolated from closer localities of the Aysén region as one population, we performed a pair-wise F ST analysis using ARLEQUIN 3.5.1.2 software (Excoffier et al. 2005). Results revealed that Ceres was the only differentiated population with all other (Table 1). Considering this result, we separated the clones isolated from Aysén region in two populations: one constituted by Ceres clones, hereafter denominated as north Aysén population; and the other composed for clones isolated from Valverde, C. Davis and Ester, hereafter called south Aysén population. Thus, for the estimation of genetic population diversity and structure we considered four populations, distinguishing between south and north Magallanes populations, and those two from Aysén. Population´s genetic diversity was estimated as the proportion of polymorphic loci (S), the observed number of alleles (na) and was performed using POPGEN 1.32 software (Yeh et al. 1997). Meanwhile, genetic diversity (h) (Nei's) (Nei 1987) estimation was performed using ARLEQUIN 3.5.1.2 software. In addition, population genetic diversity (h) differentiation was determined by analyses of variance (ANOVA, Type II model) and Tukey´s HSD multiple-comparisons. This analysis was performed within a General Linear Mixed Models framework (GLMM) considering population as a factor with a random effect to account for the lack of independence (pseudoreplication) among populations, and using a negative binomial model for residual distributions. The null hypothesis (H0) was assessed using a likelihood ratio test (Venables & Ripley 2002) and rejected at 0.05 significance level (α). To do this, genetic diversity (h) was estimated separately by each locus and the GLMM were performed using R software (R Core Team 2017) with the package "lme4" (Bates et al. 2015). Furthermore, the population genetic structure of A. catenella was estimated by several analyses. Thus, first we performed a molecular analyses of variance (AMOVA) locus-by-locus, and the significance of the global F ST value was tested using a non-parametric approach with 10,000 permutations. Moreover, to determine the level of differentiation among populations pair-wise F ST values were estimated. Both analyses were performed with ARLEQUIN 3.5.1.2 software (Excoffier et al. 2005). Identification of the genetic population and genotypic admixture were determined using Bayesian cluster analyses implemented in STRUCTURE 2.3.4 software (Pritchard et al. 2000, Falush et al. 2007). This was done with two simulations performed with different allele frequency models (independent and correlated), but both simulations were performed with an admixture ancestry model and without location as previous information (local prior option). Correlated allelic models assume that the frequencies among populations are to likely be similar due to migration or shared ancestry; meanwhile independent allelic model says that frequencies in different population to be reasonably different from each other (Pritchard et al. 2010). To perform the analyses both burning length and MCMC were setting with 500,000 repetitions and 1 to 10 K populations was used. Results were assessed following the mean of estimated log-likelihood over K method and the analyses and visualization were performed using the web app Pophelper 1 (Francis 2016). Besides, a principal component analysis (PCA) based in the Nei's genetic distance among populations was performed with the GENALEX 6.5 software, using a binary model for a haploid organism (Peakall & Smouse 2012). Finally, ARLEQUIN software was used to perform a Mantel test and evaluate isolation by distance (with 10,000 iterations). Thus, the Y matrix used corresponded to F ST values and the X matrix to the geographical distance (km).  correlated and independent allele frequency models. The mean of estimated log-likelihood over K supported K= 2 populations using both models (mean L(K): -11044 and -11234, respectively. Fig. S1). In the correlated allelic frequencies model all populations exhibited a degree of admixture in every clone membership with a major

Genetic distance
After evaluation and filtering of the AFLP raw data 770 loci were retained for genetic analysis. No identical genotypes were observed among the 33 clones characterized, and these showed high levels of genetic distance ranging between 0.118 (Canal Davis vs Ester) and 0.530 (south Magallanes vs north Magallanes).

Population diversity and structure
The genetic diversity parameters showed a gradient from the Magallanes populations to the Aysén populations (Table 2). Lower genetic diversity (h) was observed in the south Aysén population (mean: 0.166) and higher genetic diversity was observed in both south and north Magallanes populations with mean values of 0.311 and 0.312, respectively. Moreover, significant differences in the genetic diversity were found among populations (GLMM, ANOVA Chi square= 45.945, P < 0.01), being the south Aysén population the only differentiated (P < 0.05, Tukey post hoc. Fig. 2). Genetic variability of Alexandrium catenella populations was significantly structured (AMOVA, global F ST = 0.027, P < 0.05). Genetic variation was explained in a 2.73% between populations, meanwhile a 97.29% of variation was observed within populations (Table 3). Furthermore, pairwise F ST showed that south Aysén population was the only significantly (F ST , P < 0.05) different population (Table 4). Population genetic identification and genotypic admixture were detected with Bayesian cluster analyses using both the  proportion of the genetic black population in the south Aysén population (Fig. 3). Similarly, using the independent alleles frequencies model a high degree of admixture was observed across populations with a major proportion of the black genetic population in the south Aysén population (Fig. 3). Concordantly, the PCA analysis separated the population of the south Aysén from the other populations. This population distribution was sustained by the first axes which explained the 98.1% of the variability (Fig. 4). Finally, the Mantel test did not shown isolation by distance (P = 0.504), obtaining a regression coefficient of 0.0001; correlation coefficient of -0.236; and determination of Y by X of 0.056.  Figure 3. Bayesian clustering analysis of Alexandrium catenella clones isolated from Patagonian fjords. a) Results of simulation using correlated alleles frequencies models (K= 2). b) Results of simulation using independent alleles frequencies models (K= 2). In both simulations no sampling location is used as previous information. Proportion of the clone membership to the genetic populations is illustrated by a bar plot with different colors. Below are the names of populations from which the clones were isolated / Análisis de agregación Bayesiana de los clones de Alexandrium catenella aislados desde los fiordos Patagónicos. a) Resultado de la simulación utilizando un modelo de frecuencias alélicas correlacionadas (K= 2). b) Resultados de la simulación utilizando un modelo de frecuencias alélicas independiente (K= 2). Para ambas simulaciones no se utilizó la localidad de muestreo como información previa. La proporción de pertenencia de cada clon a las poblaciones genéticas se ilustra por un gráfico de barra con diferentes colores. Debajo están los nombres de las poblaciones desde donde fueron aislados los clones

Multilocus linkage disequilibrium
Significant statistical associations between alleles at different loci were observed in all populations (P < 0.05).

The lowest values of d r index was observed in the north
Magallanes with values of 0.008, whereas the highest was observed in the south Aysén population with values of 0.020 (Table 5).

Discussion
Geographical distribution of Alexandrium catenella population genetic diversity parameters showed a gradient from south to north. Considering the temporal record of A. catenella in the Patagonian fjords (Guzmán et al. 2002, Varela et al. 2012 the Magallanes south population could be the origin point of the expansion, where higher genetic diversity should be expected. The decrease in genetic diversity towards the North could be associated with founder effect, a process expected to be especially strong in small populations or in the case of a small number of colonizers (Hallatschek et al. 2007, Excoffier et al. 2009). Similar explanation had been discussed for A. catenella populations in the Thau Lagoon, France (Masseret et al. 2009), where low genetic diversity was attributed to the introduction of exogenous population from Japan in a new habitat. Besides, in Alexandrium 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 Netherlands localities in 2012, which represent an example of an established and founder population respectively. Furthermore, this genetic diversity and variability differentiation was observed despite the low number of clones characterized (e.g., ranged from 4 to 6 clones per Japanese locality; and 1 to 46 clones per French locality; Masseret et al. 2009. Ranged from 15 to 28 clones from Netherlands, and 5 clones from Baltic Sea; Brandenburg et al. 2018). Thus, based on the genetic diversity pattern herein observed, the northward expansion hypothesis through the Chilean Patagonian fjords system, originating with the most southern populations, seems to be supported.
In this study, low levels of genetic divergences were observed among south Magallanes, north Magallanes and north Aysén populations, as indicated by several genetic analyses including, non-significant isolation by distance (Mantel P > 0.05), lowest pairwise F ST , similar genotype assignation to genetic populations (Bayesian clustering analyses), and the populations distribution derived from PCA. In the context of the historical record, A. catenella was first reported in 1972 in the extreme of southern Chile, nine years later the species was observed throughout the Magallanes region, and in 1992 was recorded in the Aysén region (Guzmán et al. 2002). Since then, recurrent A. catenella HABs events have been triggered (Guzmán et al. 2002) which could promote progressive gene flow from south to north Magallanes populations. Despite the distance (approximately 880 km of coastline) and the intricate coastal geography between the north and south Magallanes populations, high gene flow seems to dominate the population dynamics in this region. As has been previously described in a bloom event, A. catenella cells can disperse northward by advection mechanism produced by superficial water drift, the current systems of the channels and wind forces (Molinet et al. 2003, Sievers & Silva 2003. Thus, under the right environmental and meteorological conditions, these blooms could spread, for example ca., 250 km in 4 months (Molinet et al. 2003). Indeed, the last severe bloom of A. catenella trigged in 2016 showed that vegetative cells can spread over approximately 562 km in 4 months (Paredes et al. 2019). This bloom occurred between December 2015 and May of 2016 and it was extended from northern of Aysén region to the Los Ríos region (i.e., from 44°S to ca., 43° by the inner sea coast of Chiloé; and from 44°S to ca., 39.7°S by the exposed Pacific coast) (Paredes et al. 2019). In a bloom scenario, the cells dispersion even would overwhelm geographical barriers such as the Peninsula of Taitao, which had been invoked to explain the biogeographic discontinuity or genetic lineages divergences in other species such as bulk kelp and cnidarian throughout the Patagonian fjords, (e.g., Häussermann & Försterra 2005, Fraser et al. 2010. Despite the scarce monitoring in the Pacific Ocean coast of the presence of A. catenella vegetative cells, our results suggest that the coastal populations dynamic and blooms of A. catenella population can be relevant to generate high levels of gene flow among the distant populations. The pattern of genetic divergence among A. catenella populations indicates that the northernmost individuals from the north Aysén population have similarly mixed ancestry to those from the Magallanes populations, while the most protected population of south Aysén diverged from them. Considering this pattern and the genetic diversity differentiation, we can propose a geographical route for the species range expansion distribution and/or interpopulation connectivity. Thus, vegetative cells generated from the Magallanes population perhaps migrated along the open coast of the Pacific arriving in the northern part of the Aysén region around the area close to Ceres. From this area, A. catenella could colonize southward following the Moraleda channel, the main channel of the Aysén region. According to the available oceanographic information, in the Moraleda channel (from 43.5 up to 46°S) the surface layer of the water mass (between 0 and 20-30 m) is characterized by a water structured in different proportions by Estuarine Waters (EW), which are waters coming out from fjords, and an intermediate mass water that corresponds to the Subantarctic Waters (SAAW) (Sievers & Silva 2003). The SAAW enter through the Boca del Guafo and move southward through the Moraleda channel (Silva et al. 1995, Rignot et al. 2003, Sievers & Silva 2003. As indicated by Molinet et al. (2003), the A. catenella blooms could be originated between the intermediate SAAW and the surface, and its dispersion would occur in the surface water mass due to wind advection. This general current pattern in Aysén suggests a complex connectivity between the Magallanes and Aysén populations, mainly associated with the SAAW.
In addition, the high and significant levels of linkage disequilibrium observed in the south Aysén population indicate that species life-cycle characteristics and/or selection could influence the inter-population divergence. Several works have shown that the spatial and temporal differentiation can occur within a bloom coupled with linkage disequilibrium (Erdner et al. 2011, Richlen et al. 2012. Within a heterogeneous environment and with fitness differences among strains, asexual reproduction may cause rapid changes in genotypic frequency (Erdner et al. 2011, Richlen et al. 2012, Dia et al. 2014. Although, considering the high population size involved in severe bloom, the genetic diversity can be maintained during blooms even if organisms primarily reproduce clonally (Dia et al. 2014). Furthermore, if adverse environmental conditions induce encystment during bloom (Brosnahan et al. 2016) the selective effects on the proportions of different genotypes could be magnified, reducing or eliminating lineages which are less successful under the prevailing conditions (Erdner et al. 2011). Thus, differentiation could result from the selective value of some clonal lines amplified by asexual reproduction or by differential encystment of some clonal lines (Paredes et al. 2019). The Aysén clones used in this work were isolated in 2009, when one of the most important outbreaks in terms of abundance and distribution occurred, covering a wide geographical area from the Aysén region to southern part of Los Lagos region, from ca., 43º45' to 46ºS (Mardones et al. 2010). Throughout the Chilean Patagonian fjords A. catenella inhabits a heterogeneous environment, with temporal and latitudinal patterns defined by biotic and abiotic factors (Camus 2001, Sievers 2008, Niklitschek et al. 2013, Iriarte et al. 2014, which would exert significant selective regimes. For example, the salinity gradient in the surface water through Moraleda channel in Aysén could be a selective pressure resulting in the relative dominance of only a few clonal lineages. Thus, life-cycle characteristics along with selection could have influenced the reproductive dynamic of species causing linkage disequilibrium. However, since the AFLP technique displays neutral markers, only an indirect inference about the selection process can be made in this study, (i.e., that selective regimes may sustain a high difference in multilocus linkage disequilibrium) (c.f., Agapow & Burt 2001, Nosil et al. 2007. Alternatively, the differences between the Aysén populations could be due to geographic separation, where distinct sources of cells contribute to each population. But this hypothesis contrasts with the broad dispersal ability shown by A. catenella, so severe environmental constraint should have existed between both Aysén populations, which explain the contrasting pattern. In summary, genetic diversity gradient is according to the genetic pattern expected in a species undergoing range expansion, suggesting that this species has expanded its distribution from southern Magallanes to Aysén region. The absence of genetic divergence among A. catenella populations separated by large distances indicate high gene flow in the Pacific coast, and also seems to indicate the connection route between Magallanes and Aysén populations and/or how the species expansion would have occurred (i.e., from the Magallanes population toward the inner Aysén populations through coastal Pacific coast and Boca del Guafo). Moreover, genetic divergence between close Aysén populations indicates variation in the reproductive dynamic during a bloom scenario rather than gene flow constraints, where some linages may have reproduced predominantly asexually under certain environmental condition or/and they have experienced differential encystment. In this study genetic population inferences should be viewed with caution, due to the limited set of clones sampled that could have skewed the genetic diversity and structure (Nei 1978, Sinclair & Hobbs 2009, however the impact of northward expansion, and the vegetative cells dynamics, over the genetic patterns could be estimated with relative confidence. When working with small samples there is a potential risk of underestimate the genetic diversity and or overestimate the population subdivision (Sinclair & Hobbs 2009), but these patterns were not observed in this study.