Spatio-Temporal Pattern of Dinoflagellates Along the Tropical Eastern Pacific Coast (Ecuador)

Among marine phytoplankton, dinoflagellates are a key component in marine ecosystems as primary producers. Some species synthesize toxins, associated with human seafood poisoning and mortality in marine organisms. Thus, there is a large necessity to understand the role of environmental variables in dinoflagellates spatial-temporal patterns in response to future climate scenarios. In that sense, a monthly four-year (2013-2017) monitoring was taken to evaluate dinoflagellates abundances and physical-chemical parameters in the water column at different depths. Sampling sites were established at 10 miles in four locations within the Ecuadorian coast. A total of 97 taxa were identified, corresponding to 7 orders, 19 families, 28 genera. Eight potentially harmful genera were registered but no massive blooms were detected. The most frequent dinoflagellates were Gymnodinium sp. and Gyrodinium sp. Environmental variables showed different mixing layer thickness and a conspicuous and deepening thermocline / oxycline / halocline and nutricline depending on annual and seasonal oceanographic fluctuations. This study confirms that seasonal and spatial distribution of the environmental variables are linked to the main current systems on the Eastern Tropical Pacific, thus the warm Panama current lead to a less dinoflagellates abundance in the north of Ecuador (Esmeraldas), while the Equatorial Upwelling and the cold nutrient-rich Humboldt Current influence dinoflagellates abundance at the central (Manta, La Libertad) and South of Ecuador (Puerto Bolivar) respectively. Inter-annual variability of dinoflagellates abundance is associated with ENSO and upwelling conditions. Climate change scenarios predict an increase in water surface temperature and extreme events frequency in tropical areas, so it is crucial to involve policy-makers and stakeholders in the implementation of future laws involving long-term monitoring and sanitary programs, not covered at present.


INTRODUCTION
The Equatorial Pacific region is characterized by a unique complex oceanographic variability with a well-defined Equatorial Front where El Niño Southern Oscillation (ENSO) events eventually occur (Santos, 2006;Gierach et al., 2012). The Ecuadorian coast is divided in two biogeographical regions, a northern "Bay of Panama" ecoregion from Azuero Peninsula to Bahía de Caráquez and a southern "Guayaquil" ecoregion from Bahía de Caráquez to the Illescas peninsula in Perú (Sullivan and Bustamante, 1999). The northern coast of Ecuador is influenced by the warm current of Panama and the current of El Niño, while the Humboldt current brings cold water to the southern Ecuadorian coast which contributes to outcrops formation in the area (Cucalon, 1989). The southern sector is the most productive area regarding fisheries and shrimp production (Pennington et al., 2006). Due to the variability found in oceanographic conditions over the Ecuadorian coast, there is a noticeable heterogeneity concerning physical, chemical, and biological parameters which leads to a high diversity in phytoplankton communities (Smayda, 2008) and the potential formation of harmful algae blooms (HABs).
The oceanographic conditions identified as drivers of HABs, are linked to those conditions associated to upwelling systems (cold nutrient-rich waters), high column water stratification (mostly associated to ENSO conditions), and long-term sea surface temperature (SST) increase due to climate change (Hallegraeff, 2010;McCabe et al., 2016;Kudela et al., 2017). Those studies that have related the particular oceanographic conditions with the occurrence of HABs in the Pacific Ocean (Brown et al., 2003;Franco-Gordo et al., 2004;Moore et al., 2010), have shown that those factors module the phytoplankton community and dynamics. In general, the phytoplankton community on the Eastern Tropical Pacific (ETP) region is characterized by relatively low levels of biomass, with the dominance of small species and rare large bloom-forming diatoms (Brown et al., 2003). Diatoms and dinoflagellates are two major groups of phytoplankton that flourish in the oceans (Abate et al., 2017), mostly registered in the south-eastern Equatorial Front (Jimenez-Bonilla, 1980;Torres and Tapia, 2002), and upwelling systems (Shulman et al., 2012). Particularly, dinoflagellates exhibit environmentally induced adaptation and survival to changing environmental conditions (Kudela et al., 2010).
Among the dinoflagellate community, there are some toxin producer dinoflagellates which may impact negatively marine ecosystems causing massive marine organism mortality, cause economic damage in coastal locations and affect human health by seafood consumption (Li et al., 2014). In addition, ballast waters contribute to the transfer of phytoplankton in commercial and oil tankers resulting in the dispersion of non-native phytoplankton species, increasing the risk of phytoplankton proliferation (Hallegraeff, 2003). Moreover, physiological changes may affect phytoplankton communities in response to climate change (Poloczanska et al., 2016) with consequences in their distribution and abundances worldwide. In that sense, it is important to perform studies related to phytoplankton biodiversity characterization and its spatio-temporal variability to guarantee an optimal integrated coastal management, contributing to future decisions from policy-makers to optimize waters vigilance.
Some countries of Latin America, such as Chile or Mexico, have implemented HABs monitoring regarding seafood production safety (Daguer et al., 2018;León-Muñoz et al., 2018). However, in Ecuador there is a lack of a national regulation plan, related to phytoplankton observations. Nevertheless, the Oceanographic Institute of Ecuador (INOCAR) has performed periodic monitoring since 1989 along the Ecuadorian coast. Torres (2015) observed the presence of algal proliferations along the Ecuadorian coast, which were more frequent in the Gulf of Guayaquil, where most aquaculture activities are located. Therefore, through the analysis of the INOCAR data bases, this study aims to provide information related to dinoflagellates abundances dynamic across a broad latitudinal area from 2013 to 2017 that constitute a base-line for future investigations.

Study Area and Samples Collection
The study area was located between 0 • 89 North and 03 • 036 South, in the ETP Ocean. Along the study area, four geographic points were settled as sampling stations: three of them located 10 miles from the coast: (1) Esmeraldas (Galera San Francisco), (2) Manta, (3) La Libertad (Santa Elena Peninsula), and (4) Puerto Bolívar (Figure 1). All stations have an approximate depth of 100 m. except for Puerto Bolívar station, which is located approximately 27 miles, inside the submarine platform in a pit that is 9 m deep, northeast of Santa Clara Island in the Gulf of Guayaquil.
At each sampling station, 200 mL water samples were collected in plastic bottles from surface, 10,20,30,40,50, and 75 m depth using a 3 Liter volume Vand Dorn bottle, then preserved in neutral Lugol's solution (final concentration 0.4%). Samplings were performed monthly from February 2013 to December 2017, between 09h00 to 12h00 am.

Phytoplankton Identification and Abundance Estimation
Water samples were settled in 25 mL volume Utermöhl chambers (Utermöhl, 1958), and observed in a Leica (ML) inverted microscope after 24 h. Cell abundance estimation was performed in horizontal transects of the chamber at 400× magnification. Results are reported in cell. L −1 . The validity of names of the different taxa was checked on the World Register of Marine Species (Horton et al., 2018). Dinoflagellate species identification was based on Steidinger and Jangen (1997);Hoppenrath et al. (2009)Hoppenrath et al. (2014; Omura et al. (2012), and Lassus et al. (2016). The list of toxic species was checked on the Taxonomic Reference List of Harmful Micro Algae (Moestrup et al., 2009).
In Esmeraldas and Puerto Bolívar no sampling was performed during February to December 2016 due to economic issues.

Environmental Variables
Temperature and salinity vertical profiles were obtained using a SBE 19plus V2 SeaCAT Profiler CTD, at each sampling site. Data was processed using the SEABIRD software. Water samples for dissolved oxygen and nutrient estimation were collected with Van Dorn bottles at same depths as for phytoplankton analysis. Determination of dissolved oxygen for each sample was carried out in situ following the protocol for Dissolved Oxygen Test PEE/LAB-DOQ/01-INOCAR based on Apha (2005). Water samples for nutrient analysis were filtered with 0.45 µm millipore filters and kept at −20 • C. Nitrate, silicate and total phosphate were analyzed by techniques described in Test PEE/LAB-DOQ/01-INOCAR based on Strickland and Parsons (1972). For phosphates, as values in Equatorial Pacific waters are below 3 µmol. L −1 , the analysis followed is detailed in Part II of the protocol cited above. In Esmeraldas and Puerto Bolívar there is no data for dissolved oxygen and nutrients from February to December 2016 due to economic issues. In this study N/P is the ratio between Nitrate and Phosphate.

Data Analysis
Before any formal statistical analysis, an exploratory data analysis (EDA) was conducted to both, environmental and biological databases in order to avoid further statistical problems (i.e., proper family distribution of the variable, homoscedasticity, independence of the values), as suggested by Zuur et al. (2010). Beside this, the environmental variables were inspected in order to avoid collinearity (i.e., correlation among explanatory variables), thus all variables who showed a Pearson correlation coefficient (rho) higher than 0.7 were dropped from the analysis (Supplementary Figure S1). Community analysis considered the calculation of the common Shannon-Weaver diversity index and its spatio-temporal variability. From an ecological perspective those species with a percentage of occurrence higher than 2% were related with the environmental variables using a canonical correspondence analysis (CCA) for each year (ter Braak and Verdonschot, 1995) through the "cca" function of the vegan package (Oksanen et al., 2018). The significance of the models were achieved through the function "anova.cca" (from the vegan package, op.cit) which performs an ANOVA like permutation test using 999 permutations (Borcard et al., 2018). All the statistical analysis were conducted under the R software (R Core Team, 2018; version 3.4.4) and supervised by OneMind-DataScience 1 .

Dinoflagellate Occurrence
During the studied period, a total of 102 dinoflagellate taxa were identified, corresponding to 8 orders, 22 families and 31 genera ( Table 1). The genera Tripos and Protoperidinium were the most diverse, comprising 22 and 16 species, respectively ( Table 1).
With respect to the temporal variability of the species mentioned above (most representative in terms of occurrence), Gymnodinium sp, Gyrodinium sp, Oxytoxum spp, and O. scolopax were present during the entire study period (Table 1). Gonyaulax spp was present also for both seasons (wet and dry seasons) 1 https://onemind-datascience.com/ and all years except 2016, which was not observed (Table 1). Gyrodinium acutum and Scripsiella spp was found in both wet and dry seasons, but only for the 2015 to 2017 period (Table 1).
Protoperidinium spp was found also in both wet and dry seasons but for the 2013-2015 period (Table 1). Instead, for the years 2016 and 2017, Protoperidinium spp was only observed during the dry season. P. micans was only found in dry season for 2013 and for wet season in 2016; instead during the period 2015-2017, P. micans was present in both seasons (Table 1).

Spatio-Temporal Variability
The two most dominant species, Gymnodinium sp and Gyrodinium sp, were more abundant in Manta, with bloom formation up to 1.2 × 10 6 cell. L −1 for Gymnodinium sp in April 2017 (Figure 2A) and 5.4 × 10 4 cell. L −1 for Gyrodinium sp in February 2016 ( Figure 2B).
Regarding cell abundance estimation, 2013 registered for the entire water column studied (0-75 m depth) the lowest abundances comparing with the rest of the years (Figure 3). For 2013, the maximum values were found at surface in La Libertad (mostly corresponding to March (2.8 × 10 5 cell. L −1 ) and at 10 m depth in Puerto Bolívar (mostly in February, 2.0 × 10 5 cell. L −1 ). In 2014, abundances were higher in the first 10 m of the water column except for Esmeraldas where abundances were similar in the first 30 m, being lower comparing with the other sampling sites for the same year (Figure 3). The maximum abundances in 2014, were observed at 10 m in Manta and at surface in Puerto Bolivar (with the maximum value corresponding to May, with 2.9 × 10 5 cell. L −1 and 3.8 × 10 5 cell. L −1 , respectively). In 2015, Esmeraldas registered again the lowest concentrations in the first layers, comparing with the other sampling sites for the same year. The same vertical pattern was observed, as in 2014, with similar abundances from surface until 30 m (Figure 3). Maximum of the year occurred in La Libertad at surface, in March (1.1 × 10 6 cell. L −1 ) and at 20 m depth in Puerto Bolívar in July (maximum value of 4.5 × 10 5 cell. L −1 ). No peaks were registered during the year 2016, but data is not available for Esmeraldas and Puerto Bolivar. Thus, the maximum concentrations registered were observed in the first 20 m for Manta and La Libertad sites (Figure 3). The maximum abundances were found at surface and at 75 m depth for Manta (in February with 1.7 × 10 5 cell. L −1 and 1.4 × 10 5 cell. L −1 , respectively) and at surface for La Libertad site (with maximum value for March with 1.0 × 10 5 cell. L −1 ). The highest concentrations of the sampling period were reported in 2017 (Figure 3), mostly between January and June. The maximum abundance was found in Manta, in April at surface, 1.5 × 10 6 cell. L −1 and at 10 m depth in Manta with an abundance of 4.0 × 10 5 cell. L −1 . In February, at the first 10 m layer, Manta and Esmeraldas sites registered also high abundance values (∼3.5 × 10 5 cell. L −1 and ∼ 2.5 × 10 5 cell. L −1 , respectively). The pattern of decreasing dinoflagellate concentration with depth, (very significant from 20 m onward) observed in the vertical abundances' distribution for the rest of the study period, in 2017 is not so marked (Figure 3). In 2017, high abundances are observed in more deep layers (30-40 m) respecting to the

Species Richness
Species richness values ranged from zero to a maximum of 21 species recorded in Puerto Bolivar during 2015 (August). Higher values were found in 0 m with a negative trend with depth (Figure 4) Figure S2).

Temperature and Salinity
In general, water temperature was higher in Esmeraldas compared to the rest of the stations during the studied period with a maximum in surface of 28.86 • C in May 2017 ( Figure 5A).  the lowest values occurred in September in all depths of the water column ( Figure 5A).
In 2016, lowest temperatures were reported in March in Manta for the first 20 m and in April in Puerto Bolívar from 10 m (Figure 5A).
During 2017, lowest temperatures were registered in February and March for Esmeraldas, Manta and La Libertad, while in Puerto Bolivar there were also registered in April ( Figure 5A).
Contrary, higher temperatures (>25 • C) registered in 2013 occurred from April to December in the first 40 m in Esmeraldas and Manta, however, in La Libertad and Puerto Bolívar, warmer waters were limited to the first layers, 20 and 10 m, respectively ( Figure 5A).
In 2014 and 2015, Esmeraldas and Manta registered warmer waters in the 50 m of the water column from April to December, evident in La Libertad and Puerto Bolivar (Figure 5A).
The thermocline was deeper in northern stations, Esmeraldas and Manta (between 30 and 50 m), while in La Libertad and Puerto Bolivar it was between 10 and 30 m ( Figure 5A). Annual vertical profiles of monthly mean temperature are available in Supplementary Figure S3.
Lower salinity values (<32) were observed mainly in

Nutrients and Dissolved Oxygen
Regarding nutrients, highest nitrate concentrations (5 a 27.6 µg-at·L −1 ) were registered from 20 to 50 m, mainly during the warm season (February and March) and dry season (From July to October) in all sampling sites ( Figure 6A). Lowest nitrate concentrations (<1 µg-at·L −1 ) were observed in surface, up to 30 m and there was not a seasonal trend among sampling sites. In general, Esmeraldas registered lowest values compared with the rest of sampling sites. Annual vertical profiles of monthly mean nitrate are available in Supplementary Figure S5.
Concerning phosphates, highest concentrations (>1 µg-at·L −1 ) were registered in 2013 in Manta and Esmeraldas in the first 50 m of the water column in February and March and between 30 and 50 m from May to September. For the rest of the years, higher values decreased between 0.3 and 1.9 ug.at/l and no seasonal pattern was observed among stations ( Figure 6B).
Lowest phosphate concentrations (<0.2 µg-at·L −1 ) were mainly observed in Esmeraldas without a particular pattern. Manta also registered low values but less frequency. Annual vertical profiles of monthly mean phosphate are available in Supplementary Figure S6.

Relationship With Environmental Variables
For those species with a percentage of occurrence higher than 2.8% the association between the environmental variables and the structure of dinoflagellates was inspected through CCA per year. Results showed for year 2013 a significant proportion of total constrained inertia (TCI) explained by the environmental variables (TCI = 8%, p-value < 0.05) with the first two axis explaining more than 80% of this constrained inertia ( Table 2). Permutation showed that salinity, phosphate and temperature significantly affected the community structure for this year ( Table 3). The community composition (CCA analysis) showed that the species, Protoperidinium simulum and Prorocentrum compressum, were positively correlated with salinity and with N/P (this last being non-significant). The taxa Oxytoxum sp and Protoperidinium sp. were positively associated with high values of nitrite. The group formed by Gymnodinium cf. catenatum, Oxytoxum turbo and P. micans were positively associated with high values of temperature and negatively with salinity and N/P values. N/P is the ratio between nitrate and phosphate in this study, low N/P ratio have been associated to HAB occurrence. In relation with salinity the species O. scolopax, Gyrodinium sp and Gonyaulax polygramma showed a negative association with this variable. Finally, the group formed by the species Gonyaulax sp, Gymnodinium sp, Tripos fusus and Tripos furca were negatively associated with nitrite and phosphate, but showed an unimodal response with the rest of the variables (a central position in the CCA-Triplot, Figure 7A).
For year 2014 a lower than previous year, but still significant proportion of TCI explained by the environmental variables (TCI = 5.8%, p-value < 0.05) with the first two axis explaining 78.8% of this constrained inertia ( Table 2). For this year (same as 2013) salinity, temperature and phosphate significantly affected the community structure (Table 3). Moreover, a clear positive correlation between phosphate with nitrite and N/P (both being non-significant) and those three variables showed a negatively association with temperature ( Figure 7B). The association of species with environmental variables showed that Dinophysis sp was positively associated with N/P and nitritephosphate variables. With respect to salinity Oxytoxum turbo was positively associated whereas T. fusus, O. scolopax and P. simulum were negatively associated. T. furca and C. catenatum (in a lesser extent) showed a positive association with temperature. The rest of the species showed a central position in the CCA-Triplot representing unimodal response to the environmental variables ( Figure 7B).
For year 2015 the proportion of TCI explained by the environmental variables was also significant (TCI = 5.8%, p-value < 0.05) with the first two axis explaining the highest amount of this constrained inertia along the study period (89.4%, Table 2). During 2015 the temperature, the N/P ratio and salinity were those variables which significantly explained the dinoflagellate's community structure (Table 3). A positive relationship was observed among salinity and nitrite and phosphate and a negative association between temperature and N/P (Figure 7C). From the species-specific relationship, the CCA analysis revealed a positive association between Dinophysis sp and P. simulum with N/P ratio and negatively with temperature. The species Gymnodinium cf. catenatum was associated with positive values of nitrite and phosphate and in a lesser extent, showing a non-linear response, with salinity. Finally, the species P. micans, P. compressum and T. furca were associated to low values of salinity (negative relationship) ( Figure 7C). Interestingly, this year none species was associated with high, but low temperature values.
For year 2016 and 2017 a decrease in the TCI from 7.8% to 5.5 was evident. However, for both years the model was able to significantly explain this constrained inertia with the first two axes explaining 82.6% and 65.1% of this constrained The total number of observations (n), the total model inertia (TMI), the proportion of total constrained inertia explained by the environmental variables [TCI (%)], the model significance at α = 0.05 (p-value ), the proportion of constrained inertia explained by axis 1 [CI-1 (%)] and by axis 2 [CI-2 (%)], and the total proportion of constrained inertia explained by the first two axis CI (%) is showed for each year.
inertia, respectively ( Table 2). Particularly during 2016 the environmental variables who were significantly associated with the structure of the dinoflagellate community were temperature, nitrite and phosphate ( Table 3). Salinity and phosphate were positively correlated and temperature and nitrite were negatively correlated. The CCA analysis for this year showed that Gymnodinium cf. catenatum was correlated positively with salinity and negatively with nitrite. The species Gyrodinium acutum was associated with high values of phosphate with a nonlineal response (center position of the species related with the environmental vector). The species C. catenatum, O. scolopax, Scripsiella sp, G. polygramma, T. furca, Protoperidinium sp and Gymnodinium sp. showed a negative association with phosphate and Oxytoxum sp, Gyrodinium spirale, G. polygramma, Scripsiella sp, O. scolopax, and Gymnodinium sp a negatively association with nitrite ( Figure 7D). During 2017 all the environmental variables were significant when related to the observed dinoflagellate community structure ( Table 3). Among the environmental variables, a negative association was evident for salinity and phosphate ( Figure 7E). From the community composition perspective, the species P. simulum, T. fusus, Gymnodinium cf. catenatum, and Prorocentrum dentatum showed a clear positive relationship with nitrite and negatively with salinity and with the N/P ratio. The species T. fusus, P. compressum, O. scolopax, and Gyrodinium sp were positively associated with phosphate and negatively with salinity. For this year, non-species were related positively with temperature, but some species showed a negative association with this environmental variable (e.g., Oxytoxum turbo, Gyrodinium spirale, Oxytoxum sp, and Gyrodinium sp).
Interestingly, the total model inertia (TMI), which represents the total variability n-dimensional inside the model, remained more or less constant for the first 3 years (2013-2015) ranging from 4.1 to 4.9. However, during years 2016 and 2017 this TMI was significant lower with values for 2016 of 1.95 and for 2017 a TMI of 1.82 ( Table 2).

Dinoflagellate Community
In this study, there were reported a total of 97 taxa, corresponding to 8 orders, 22 families, and 31 genera during the sampling period from 2013 to 2017 in 4 stations in the coast of Ecuador at different depths. Phytoplankton dynamics across the ETP coast is more characterized in other countries and there is a lack of information in the central area corresponding to the Ecuadorian coast. Interestingly, recent studies regarding phytoplankton distribution in response to ENSO events have been published , but, unfortunately, species identification is not indicated as a result. Another recent study explored the oceanography of red tides using remote sensing data from 1997 to 2017, confirmed that potential HABs have been dominated by dinoflagellates during wet season mostly at the Gulf of Guayaquil (Borbor-Cordova et al., 2019). Furthermore, dinoflagellates abundance reported in the present study was not high and the number of proliferations was low, only Gymnodinium sp exceeded 10 6 cell. L −1 in April 2017. In the coast of Ecuador, INOCAR and other researchers have reported blooms of dinoflagellates, for example Gymnodinium sp was reported in 2003,2004,2005 reaching high concentrations (1 × 10 5 cell. L −1 -1.1 × 10 7 cell. L −1 ) in the Gulf of Guayaquil both in wet and dry season (Torres and Tapia, 2002;Torres et al., 2017).
It is important to mention that species identification resulted difficult to perform without specific equipment and adequate formation. In that sense, some of the taxa found in the area were not identified at species level. Dinoflagellates identification under inverted microscope based on morphological features constitutes a challenge since some species from several genera share the same plate pattern and overlap in size (e.g., the genus Ostreopsis - Carnicer et al., 2016). In the last decades, molecular techniques have allowed to perform accurate identifications in dinoflagellates using rDNA sequences (Litaker et al., 2007). Unfortunately, there is no molecular characterization of marine dinoflagellates in the Ecuadorian coasts, with the exception of Ostreopsis cf. ovata (Carnicer et al., 2016). It is important for Ecuadorian institutions to invest in the implementation of molecular techniques for future investigations, especially important to identify potentially toxic species.
Higher species richness in the water column was observed in surface layers, with a maximum of 21 taxa. These results are in agreement with a study performed in 2015 at the southern coast of the country, in El Oro province, where they reported species richness average between 28 and 23 dinoflagellates in surface samples (Prado et al., 2015;. The time series length of the present study avoids us to conduct a formal time series Significant environmental variables (terms) are showed in bold and the value of the F statistic (an p-value) are also given. n.s., means non-significant term.
analysis to statistically prove this increase pattern observed. In addition, some oscillatory behavior is evident and probably associated with the interannual variability pattern .
In 2015, associated with El Niño ENSO event, higher dinoflagellate species richness values were registered during the wet season in La Libertad, Puerto Bolivar and Esmeraldas, reaching a maximum of 28 species in dry season in Puerto Bolivar, that could be linked with upwelling and higher nutrients concentrations (Borbor-Cordova et al., 2019). In 2017, higher values are also observed in each station with the highest found in Puerto Bolivar in agreement with Prado et al. (2015), associating the event with the presence of the Humboldt current. The higher amount of species during El Niño ENSO events can be linked with the transport of species from the advected waters due to the current, as well as upwellings occurring in the equatorial region (Pennington et al., 2006;Gierach et al., 2012).
Regarding potential harmful dinoflagellates, a total of 8 taxa were observed, Alexandrium sp, Cochlodinium sp., D. acuminata, D. ovum, Gymnodinium cf. catenatum, Karenia sp., P. lima and P. mexicanum but not HABs were observed. In Torres (2015), a total of 131 HABs were registered from 1968 to 2009 in the same sampling sites as the present study, where toxic dinoflagellate community was similar to the species listed above. Unfortunately, no statistics analysis was performed in Torres (2015), however, they concluded that during wet-warm season the highest number of HABs events occurred, being more abundant in the Gulf of Guayaquil, southern station. Due to the low abundances registered in the present study, CCA analysis was not possible to perform with toxic species, which do not allow to confirm this hypothesis. Furthermore, an historical reconstruction of blooms from 1997 to 2017 in the Ecuadorian coast, highlighted the presence of toxic species such as Gymnodinium cf. catenatum, P. micans, C. catenatum and Dinophysis caudata among others. As for the present study, lower abundances were registered for toxic species, slightly higher in the Gulf of Guayaquil than in La Libertad and Manta (Borbor-Cordova et al., 2019). Unfortunately, Esmeraldas was not considered in the study and data corresponded only to wet season.
During these studies performed by the INOCAR, no toxic analysis has been performed in water samples due to a lack of specialized laboratory. This fact highlights the necessity of specific equipment acquisition by the government in order to control the presence of toxins in the productive areas, not only limited to monitoring phytoplankton abundances.

Environmental Variables and Seasonality
Overall, in this study, SST and salinity anomalies, nutrients, as well as the N/P ratio were strongly associated to the spatiotemporal distribution of dinoflagellate community (Figure 7). Studies elsewhere have found that HAB have been related during positive anomalies of SST, stratification with a deeper thermocline, pulses of upwelling during summer or dry season and the biological interaction among phytoplankton community (Hallegraeff, 2010;Díaz et al., 2016), some of these conditions were found in this analysis.
The environmental variables presented in this study confirm the main patterns of current systems on the ETP, which are driven by cold Humboldt Current in the central and south coast of Ecuador (La Libertad and Puerto Bolivar), associated with the equatorial upwelling (EU) system and the warm Panama current toward the north of Ecuador (Esmeraldas), these oceanographic settings drives the biogeochemical conditions that may promote extensive phytoplankton growth and potential HABs (Cucalon, 1989;Pennington et al., 2006). The north of the Ecuadorian coast is under the influence of the inter-tropical convergence zone (ITCZ), receiving higher levels of precipitation as well as runoff from Panama Gulf generating a low salinity and warmer temperatures which lead to stratified conditions, shallow halocline-nutricline (30-40 m), and low nutrients concentration (Pennington et al., 2006). Those conditions are reflected in the northern station of Esmeraldas, where minimum dinoflagellate abundances were reported (Figure 2 and Supplementary Figure  S2) compared to other stations, highlighting the influence of the oceanographic characteristics and strong seasonality of the ETP (Chavez et al., 2002).
Oceanographic studies have established the seasonal pattern of the coastal EU at the center of Ecuador, and the upwelling from the Humboldt Current on the south of the Ecuadorian Coast (Pennington et al., 2006). Seasonality of those upwelling system are reflected in all stations (based on higher nutrient concentrations and lower temperatures) except at Esmeraldas, being deeper and high nutrient concentration during wet season (December to May) but also between August to September. It is expected that Manta is mostly influenced by the EU while at La Libertad and Puerto Bolivar it is the upwelling of Humboldt Current, considered one of the most productive of
the world (Chavez et al., 2002;Pennington et al., 2006;Oyarzún and Brierley, 2018). In this study, highest concentrations of dinoflagellates were found in Manta (April 2017), La Libertad (March 2015), and Puerto Bolivar (May 2014 and May 2017), during wet season. During this period, there is a large amount of nutrients coming with the intense runoff from the agricultural Guayas Basin (Borbor-Cordova et al., 2006) and interestingly, in the present study it is reported higher levels of nutrients, mainly nitrate and phosphate (Supplementary Figures S5, S6), coming from deeper water, suggesting advection from the coastal Humboldt Upwelling in dry season with peaks on August. Thus, data on this study suggests that a combination of extreme hydrodynamic conditions, climate variability associated to the warm phase of the Pacific Decadal Oscillation (PDO) and ENSO in the coast of Ecuador (2015-2016), influencing salinity and SST anomalies and nutrient limitation conditions, led to shifts in the phytoplankton community and maybe to HABs. In this region phytoplankton productivity and growth are influenced by the upwelling dynamics which can be associated to the pycnoline-nutricline-oxycline depth, and their variability along the coast is driven by local specific conditions of coastal waves, seasonality and inter-annual variability of ENSO cycle (Montecino and Lange, 2009). The thermo-nutricline determines the supply of limiting nutrients (related to N/P ratio) to the euphotic zone and influencing the phytoplankton community assemblage during each ENSO event. This is evidenced in the present study with the inter-annual variability of dinoflagellates associated to ENSO events which is characterized by high positive anomaly of SST, changes on the timing of seasonality and strengthening (cold phase) or weakening (warm phase) the coastal upwelling (Pizarro and Montecinos, 2004;Pennington et al., 2006;Montecino and Lange, 2009). Similar dynamics have been reported in North and ETP during ENSO events (Chavez et al., 2002;Jacox et al., 2016). Previous ENSO studies pointed out that southwest winds weakness lead to a coastal upwelling decline which diminish phytoplankton blooms occurrence (Tam et al., 2008;Wang et al., 2008;Montecino and Lange, 2009). However, during warm phase ENSO event in 2015 occurred anomalous conditions related with sustained wind stress, intermittent coastal upwelling, that could be the reason why dinoflagellates were able to persist even in a nutrient limited environment, and warm stratified conditions (Smayda and Reynolds, 2001;Du et al., 2015;McCabe et al., 2016;McKibben et al., 2017;. During the period of this study (2013)(2014)(2015)(2016)(2017), 2013 is considered a normal year with the lower abundances of dinoflagellates of the study, a weak ENSO warm phase in 2014, followed by a strong warm phase ENSO in 2015-2016 (with a different seasonality of the one in 1997-1998) beginning a transition to cold phase ENSO during the dry season 2016 is extended to 2017 characterized by abundance of dinoflagellates (see Figure 3). Each year of the study develop specific environmental conditions; however, SST anomalies, salinity and nutrients persistently drove the dinoflagellates community for all the study period (see Figure7).
In 2013 a normal year with a slightly cold SST anomaly (see Supplementary Figure S9), lower salinities, and nutrients were the drivers for the dinoflagellates. Gymnodinium cf catenatum, Oxytoxum turbo and P. micans were associated to less saline, low nutrients and warm waters typical conditions during wet season in Esmeraldas influenced by the Panama Current. While P. simulum, P. compressum, Oxytoxum sp, and Protoperidinium sp are positively related to nutrient and salinity which are characteristics of upwelled waters, suggesting a preference of nutrient-rich water along the coast of Ecuador ( Figure 7A).
In 2014, considered a weak ENSO year with ONI strong anomalies (>1.5) during wet season decreasing rapidly toward December 2014 (see Supplementary Figure S9). T. fusus, O. scolopax and P. simulum were associated to low salinity conditions mostly on Esmeraldas, Manta and La Libertad suggesting iN/Put of freshwater runoff and precipitation. While Dinophysis sp was associated to nutrients from upwelled waters from EU system at the Esmeraldas station, the upwelling process is related with the nitrate distribution at 30-50 m (Supplementary Figure S5). T. furca and C. catenatum were associated to positive anomalies of temperature, suggesting that have adapted to the warm conditions on the Esmeraldas site.
Although the SST anomalies of ENSO events in 2015 were of similar magnitude to those in 1997-1998, the responses were diverse in relation to the oceanographic characteristic of the equatorial front . By April 2015, the thermocline deepened to 40-75 m until December 2015 generating stable thermal stratification, but nutrients were sustained by a weak upwelling and a pycnocline-nutricline between 30-50 m. Even though the nutrients were reduced at the surface layer (see Figures 6, 7), surprisingly the abundance levels of dinoflagellates in La Libertad and Puerto Bolivar where high between 10-30 m (see Figure 3). By August 2015, nutrients (nitrate, and phosphate) were upwelled between 30 and 50 m at La Libertad, and Puerto Bolivar. These conditions of thermal stratification, deepened thermocline, slightly upwelling, light for growth and ciliate prey (Mesodinium rubrum) appear to establish the optimal combination for the development of potential HABs such as Dinophysis sp and P. simulum (Reguera et al., 2012;Díaz et al., 2013). G. catenatum have affinity for the nitrate and phosphate and less with salinity, bloom of this specie has been related to advected vegetative populations to the coasts during the relaxation of coastal upwelling (Bravo and Figueroa, 2014). P. micans, P. compressum, and T. furca were associated to less saline water characteristic of the Panama current or influenced by precipitation at the northern stations of Esmeraldas and Manta. Considered an extreme ENSO year and warm PDO, expecting a limited algal growth, surprisingly there was a relative high abundance of dinoflagellate and higher richness reflexed in the diversity of species (21). If this increase in richness is a pattern of the ENSO need to be verified in other years.
Year 2016 shows the transition between a warm ENSO to cold and normal conditions, with ONI decreasing from 2.5 in January to -0.7 in September in the dry season (see Supplementary Figure S9). In this year, Gymnodinium cf. catenatum appeared during dry season associated with salinity and inversely with nitrite at the Manta and La Libertad, unfortunately no data in Esmeraldas and Puerto Bolivar are available. Gymnodinium cf. catenatum has been recognized to cause paralytic shellfish poisoning across regions and as a survivor by using strategies based on its motile forms and cysts in plankton assemblages and in surface sediments (Yamamoto et al., 2004;Bravo and Figueroa, 2014). Gymnodinium cf. catenatum is also identified as one of ten most damaging potential domestic target species, considering the environmental and economic impact in aquaculture and also in human health (Hallegraeff, 2003). C. catenatum is another specie that seems to be associated to the upwelling at La Libertad in dry season.
During 2017, persistent negative SST anomalies associated to upwelling and water nutrients richness lead to an increase in the productivity of the region in all the ETP (Brown et al., 2015). Concentrations of dinoflagellates increased and reached their maximum in all the stations due to conditions of a persistent upwelling specially in Puerto Bolivar station. Some of the species that were associated with upwelled nutrients were P. simulum, T. furca, Gymnodinium cf. catenatum and P. dentatum.
In the present study, the most abundant species were Gymnodinium sp and Gyrodinium sp which were present every year in both seasons of the studied period and registered up to 40 m depth. Those species correspond to free living large cell dinoflagellates, ubiquitous in the Pacific Ocean, previously registered in México, California and Ecuador (Gomez, 2005;Meave-del Castillo et al., 2012;Torres, 2015). Both species appeared in all coastal stations during the seasonal transition in the weak ENSO in 2014 and strong ENSO in 2015 (Figure 2 and Table 1) indicating their capacity to adjust to warm and deeper thermo-nutricline.
Even though some microalgae exhibit a rapid adaptability to high temperatures in short-term experiments, the physiological plasticity and genetic response of many microalgae under future environmental conditions is unknown (Hallegraeff, 2010). Organism size and elemental composition of the phytoplankton community will influence processes at the level of individuals, populations, communities and ecosystems (Finkel et al., 2010). Some researchers have found that dinoflagellates are species with a smaller surface and have the ability to grow and sustain during ENSO, ocean heatwaves, stratified conditions, and nutrient poorer waters (Smayda and Reynolds, 2003;Glibert et al., 2005).
However, there is also evidence that nutrient enriched waters can stimulate dinoflagellates blooms (Glibert et al., 2005), therefore, dinoflagellates are considered strategist and survivors either in extreme conditions.
The contrasting behavior of Tripos species in divergent lake types (with different climatic, morphometric, geological, hydrological, and trophic features) explains the existence of ecotypes of these species adapted to diverse environmental conditions and exhibiting high intra-and inter-population and morphological variability (Cavalcante et al., 2016). T. furca is also known to perform active vertical migration, depending on nutrients and water temperature in both natural and laboratory conditions, which means that different species of dinoflagellates have different ecological abilities. These abilities are possibly linked with their ecological responses for surviving under local environmental conditions. T. furca has a competitive advantage because of physiological adaptations to low nutrient concentration waters (Baek et al., 2011). Therefore, dinoflagellates can alter their vertical position in the water column by swimming, which allows them to maintain an optimal depth in terms of light or nutrients. Consequently, the actual duration of high irradiance in the coastal bays is considered to regulate the maintenance of the bloom (Baek et al., 2011).
In general, under N or P limitation or some specific N/P ratio, some toxic dinoflagellates are probably outcompeted, and toxin production may be an adaption strategy to offset the ecological disadvantages of dinoflagellates with low nutrient affinity (Smayda, 1997). How these biochemical shifts shape phytoplankton community structure and the presence of potentially toxic species need to be explored in Ecuador. Montecino and Lange (2009) proposed that during warm ENSO as 2015, there is microbial trophic web path mostly dominated by small-sized phytoplankton, in this study taxa such as Gonyaulax polygrama, P. dentatum and Gynmodinum spp. were characterized to be trophic web dominated. Still it is not very clear which are the main ecological drivers to explain interannual variability of phytoplankton community structure, including toxin producer species. In addition, mixotrophy relationships need to be explored because many of the motile bloom formers are mixotrophic protists (Stoecker et al., 2017). Some argued that PDO is a good indicator to explain those shifts on the community structure and potential toxins generation (Du et al., 2015). Unfortunately, in this study only dinoflagellates community was considered, so it is not possible to confirm this ascertain. An effort in future works need to be done in order to encompass more taxa to better understand phytoplankton distribution in the area.

CONCLUSION
Dinoflagellates abundances reported in the present study were not elevated, and no HABs were observed in the studied period. The taxa presented in this study represents a baseline in which future work may referred in order to understand certain dinoflagellate dynamics. In synthesis, the environmental variables and the oceanographic ENSO index served to understand how the current system of the ETP, their thermocline dynamics, and the coastal upwelling affect euphotic zone nutrient supply and hence, dinoflagellates abundance and richness along the coast of Ecuador.
Moreover, toxic species were detected and, considering that in future decades it is predicted an increase in seawater temperature and extreme events frequency (IPCC, 2014), phytoplankton dynamics may be influenced. In order to be able to predict future scenarios of HABs, it is important to study their distribution in relation to environmental variables in areas where extreme events may affect their proliferations. Thus, it is crucial to involve policymakers and stakeholders in the implementation of environmental policies and climate change adaptation strategies involving longterm monitoring and sanitary programs, not covered at present.

AUTHOR CONTRIBUTIONS
GT, SR, RN, and EP made substantial contributions to conception and design of the work and conducted the sampling, nutrient, and dinoflagellates analysis. OC, GT, MB-C, PDLF, and AC made substantial contributions in data analysis and interpretation and participated in drafting the article. MB-C gave the final approval of the version to be submitted.

FUNDING
This research has been developed as an international and interinstitutional collaboration among different research institutes and universities. We would like to acknowledge the following projects and institutional support: Project "Oceanic Early Warning System, " funded by National Secretary of Science and Technology (SENESCYT) and INOCAR, Project T2-DI-2014 "Climate Variability and Harmful Algae Bloom Interaction and their impact on human health on the coast of Ecuador, " funded by the Escuela Superior Politecnica del Litoral (ESPOL).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2019.00145/full#supplementary-material FIGURE S1 | Collinearity analysis of the environmental variables measured in the study. The Pearson's correlation coefficients for each pair of variables are shown in the upper-right panel with the size of the text proportional to its value. Significance code represents: * * * p-value < 0.001. Additionally, the density histogram of each variable is represented in the diagonal and a scatter plot with a smoothing spline (in red) is shown in the left-lower panel.