Marine Dinoflagellate Assemblage in the Galápagos Marine Reserve

It is likely that harmful algal blooms have increased in frequency, intensity and geographic distribution in the last decades in response to anthropogenic activities. The Galapagos Islands are renowned for their exceptional biological diversity; however, marine dinoflagellate communities have not been represented in biodiversity assessments. Therefore, this study aims to provide key information about dinoflagellate diversity and abundances, with special attention to harmful species, during a weak La Nina event in the Galapagos Marine Reserve (GMR). The study was performed during March-April 2017 and four transects were conducted at four Islands (Santa Cruz, Santa Fe, Seymour and Pinzon) representing the southern region of the GMR. Water net samples were collected at 2, 5 and 10 nautical miles (nm) from the coast, at a total of 48 sampling sites. The presence of toxic species, and their cell abundance was estimated in seven transects at 0, 15 and 30 meters of depth. A total of 152 taxa belonging to 7 orders, 22 families, and 38 genera were registered. The number of taxa found is almost three times higher than the maximum observed in previous studies. Dinoflagellate species richness among stations ranged between 53 and 23 taxa and was higher in northern sites. From the applied cluster analysis, five dinoflagellate assemblages were identified as a discrete community structure, one was found only in Santa Fe Island, which is probably related to the presence of the Equatorial Undercurrent (EUC). Regarding cell abundance estimations, low abundances were registered throughout the sampling sites and no blooms were detected. Higher abundances were registered in the northern transects coinciding with one of the most productive areas of the archipelago, situated north of Santa Cruz. Among the identified taxa, 19 of them were potentially toxic, including epiphytic species, allowing the possibility of blooms in benthic areas. This study presents the first record of several dinoflagellate species in the area (both non-toxic and harmful species) and thus, emphasizing the need for the implementation of phytoplankton monitoring programs by the government to prevent potential ecological, sanitary and economic impacts in the GMR.


INTRODUCTION
Marine microalgae are primary producers and constitute key components of marine food webs: they fix carbon and produce nearly 50% of the oxygen on the planet (Field et al., 1998). Most phytoplankton communities are directly affected by human activities, especially through the excess inputs of organic matter and nutrients to the system (Hallegraeff, 2010), usually linked to high phytoplankton biomass in coastal areas (Davidson et al., 2014). Moreover, phytoplankton dynamics are affected by climate change, due to altering environmental conditions such as timing of large-scale climate events (Boyce et al., 2010) or sustained warmer temperatures that increase stratification and reduce the resupply of recycled nutrients to the upper mixed layer (Doney, 2006). Consequently, physiological responses of phytoplankton species, such as phenology, abundance and calcification rate may be affected (Poloczanska et al., 2016), leading to a reduction in ocean productivity (Behrenfeld et al., 2006) and changes in bloom timing and community structure (Henson et al., 2018). Thus, the study of phytoplankton communities and its response to environmental drivers is crucial to gain a better understanding of the functioning of the marine ecosystem.
Particularly, some species can form toxic proliferations of cells, thereby causing harm to aquatic ecosystems, including the resident plants and animals, and substantial economic losses and disturbances in touristic areas (Berdalet et al., 2015). These, so called "harmful algal blooms" (HABs), can affect humans through direct exposure by skin contact or inhalation of aerosols, or by ingestion of seafood contaminated by toxins that have bio-accumulated along the food web. Examples of this are well known, for instance in the digestive tract of shellfish (mussels, clams, oysters, scallops) or finfishes (Berdalet et al., 2015). Dinoflagellates are of particular interest among the harmful phytoplankton species because of their high species richness, morphological diversity and strategies for adaptation to thrive in different ecological niches (Smayda and Reynolds, 2003). Most dinoflagellates have cyst-forming stages, which are in a dormant state until environmental conditions, such as temperature, salinity, light intensity or turbulence, favor their growth, resulting in abundant proliferations in a short period of time (Delwiche, 2007). In addition, there is a high number of mixotrophic dinoflagellate species, whose predatory behavior enables them to increase their nutrient uptake, allowing them to survive under conditions that are unfavorable for strict autotrophs (Stoecker, 1999). From an ecotoxicological point of view, dinoflagellates are of vital importance since they show the highest representation among toxic phytoplankton with 99 species, in contrast with the number of diatom species (29), Haptophytes (8), Raphidophyceans (4), Dictyochophyceans (3), Pelagophyceans (2), and Cyanobacteria (35) (Moestrup et al., 2009; IOC-UNESCO Taxonomic Reference List of Harmful Microalgae). These features, mentioned above, make dinoflagellates crucial to marine ecosystems.
Most dinoflagellates' toxins are neurotoxic, but others cause specific poisoning syndromes including gastrointestinal disturbances (diarrhea, nausea, vomiting), muscle paralysis, or amnesia (Hallegraeff, 2003). The most frequent and studied syndromes caused by dinoflagellate toxins are Paralytic Shellfish Poisoning (PSP), caused by Alexandrium spp., Gymnodinium catenatum, and Pyrodinium bahamense; Diarrhetic Shellfish Poisoning (DSP) produced by Dinophysis spp. and Prorocentrum spp. Poisoning syndromes not only come from planktonic species, but also from benthic species. For example, Gambierdiscus spp., can cause ciguatera fish poisoning, a disease with worldwide impact, though mainly in tropical zones (Parsons et al., 2012), and Ostreopsis cf. ovata blooms that have been seen to affect human health directly by inhalation of marine aerosols that cause respiratory irritation (Vila et al., 2016).
Nowadays, there is a big concern about proliferations of toxic species that may be increasing in frequency, and expanding their biogeographic distribution under future global warming scenarios (Hallegraeff, 2010;Glibert et al., 2014;Wells et al., 2015). Thus, it is important to survey ecological hotspots such as the Galápagos Marine Reserve (GMR), which is one of the largest and most biologically diverse marine regions in the world (Schaeffer et al., 2008;Liu et al., 2014). Its great diversity derives from the situation of the GMR lying in a complex transition zone between tropical and subtropical waters and intense equatorial and local upwelling (Palacios, 2004;Liu et al., 2014). The GMR area is largely affected by oceanic-atmospheric perturbations such as the seasonal migration of the Intertropical Convergence Zone (ITCZ) and the El Niño Southern Ocean Oscillation (ENSO) (Palacios, 2004;Liu et al., 2014). During the dry season (July-December), the ITCZ is north of the equator, south trade winds increase, and sea surface temperature (SST) diminishes. In turn, in the warm, wet season (December-June), the ITCZ migrates southward (toward the equator), the northeast trade winds become prevalent and precipitation and SST increases. In general, during El Niño events, warmer surface conditions and the deepening of the thermocline inhibit the upwelling of cooler, nutrient-rich subsurface waters to the surface. As a consequence, during an ENSO, distributions of phytoplankton communities undergo large-scale disruptions and chlorophylla (Chl a) levels diminish (McCulloch, 2011). For example, in the equatorial Pacific Ocean, during the strong ENSO event in 1997/1998, Chl a reached a very low concentration (Chávez et al., 1999), which may have influenced higher trophic levels. During La Niña events, the cool ENSO phase, strong winds and cooling of ocean temperatures cause an elevation of the thermocline that increases nutrient supply and phytoplankton productivity along the equator (Ryan et al., 2002(Ryan et al., , 2006. In addition, there is some evidence that climate change may influence ENSO to some extent (Sachs and Ladd, 2002;Cai et al., 2014), resulting in unpredictable changes in phytoplankton assemblages and HAB frequencies.
There is little data regarding taxonomic characterization of phytoplankton assemblages in the coastal area of the Eastern Tropical Pacific (ETP), thus, the diversity of phytoplankton species remains poorly described. Regarding HABs, isolated studies have documented the prevalence of Gymnodinium catenatum, Akashiwo sanguinea, and Alexandrium spp., among others (Colombia - Giraldo et al., 2014;Costa Rica -Vargas-Montero et al., 2008;Calvo-Vargas et al., 2016;South of Mexico -Maciel-Baltazar, 2015). The Oceanographic Institute of the Ecuadorian Navy (INOCAR) has conducted sporadic cruises in the GMR since 1968. Unfortunately, a reduced number of reports are available to the public and most are written in Spanish Tapia, 2000, 2002;Torres, 2002;Tapia and Naranjo, 2012;Torres and Andrade, 2014;Naranjo and Tapia, 2015). Toxic species such as Dinophysis spp., Prorocentrum mexicanum, and Ostreopsis cf. siamensis have been registered, however, the number of species indicated within the reports is greater than the list of taxa supplied, thus, there is an under-representation of the presence of dinoflagellates in the GMR. Moreover, Torres (2015) reviewed fish mortality and algae proliferation events from 1968 to 2009 along the Ecuadorian coast and in the GMR, reporting Prorocentrum gracile in some bays of the GMR, associated with fish mortalities in 1980.
The absence of a systematic monitoring program leads to a lack of scientific data to evaluate phytoplankton communities and the future risk of HABs in the GMR (Kislik et al., 2017). Considering the unique features of dinoflagellates and their relevance among HAB species, the aim of this study is to give insights into the dinoflagellate community of the GMR, specifically focusing on the island of Santa Cruz and nearby islands, with special attention to toxic species, during a weak La Niña phase in the wet season (March-April, 2017).

Study Area
The GMR is located in the ETP, ∼900 km west from the coast of Ecuador. The GMR area is influenced by cold and saltier Equatorial Surface Waters (ESW) with temperature (T) < 25 • C and salinity (S) > 34, and by the warm and fresher Tropical Surface Waters (TSW) (T > 25 • C, S < 34) (Fiedler and Talley, 2006;Sweet et al., 2007). The main currents affecting the GMR are the westward South Equatorial Current (SEC), which drives surface waters over the entire region around the Galápagos (Kessler, 2006;Schaeffer et al., 2008) and the upwelled waters of the EUC, developing a sub-surface (subthermocline) compensation against the westward SEC (Schaeffer et al., 2008). The collision of the subsurface iron and nutrientcarbon rich EUC with the Galápagos platform results in topographically induced local upwelling areas, which favor an enhanced production in many areas of the archipelago (Palacios, 2004;Schaeffer et al., 2008;Kislik et al., 2017).
The study took place in the southern area of the GMR (Figure 1) during the warm -wet season (March-April, 2017) coinciding with a weak La Niña phase. During this season, the EUC becomes stronger, shallower, and with higher salinity values than during the dry season (Sweet et al., 2007). Another source of nutrients in the GMR is partially attributed to the Island Mass Effect (IME) where iron provided by continental sediments, together with topographic upwellings, contribute to an increase in Chl a concentration (Kislik et al., 2017).

Sample Collection
Sampling was performed in four islands; Pinzón (27th March), Santa Fé (28th March), Santa Cruz (6th April), and Seymour (7th April). At each site, four transects were covered (A, B, C, D) and for each transect, three sampling stations were selected at 2, 5, and 10 nm from the shore (Figure 1). In every station, vertical net water samples (30 m) were collected in 200 mL plastic bottles using a 20 µm mesh plankton net, then samples were preserved in neutral Lugol's solution (final concentration FIGURE 1 | GMR sampling sites: Pinzón, Santa Fé, Santa Cruz, and Seymour Islands. Four transects were performed in each Island at 2, 5, and 10 nm from the coast. of 0.4%). These vertical net samples were used to identify the presence of dinoflagellates in each sampling site. In addition, seven transects (Santa Fé A, B; Pinzón D; Seymour A, B; Santa Cruz A, D) were selected based on the presence of toxic species found in net samples, for phytoplankton abundance estimation. To do this, at each sampling station (2, 5, and 10 nm) three water samples using a Van Dorn bottle were collected at surface, 15 and 30 m depth, then preserved in neutral Lugol's solution (final concentration of 0.4%).

Environmental Parameters
Vertical profiles were obtained using an EXO2 Multiparameter Sonde from Yellow Springs Instrument Company (YSI Inc.) equipped with sensors for temperature, conductivity, pressure, dissolved oxygen and pH, by deploying (at each station) from the surface down to 30 m depth. Only the downcast was selected for processing the data.

Phytoplankton Identification and Abundance Estimation
Species identification in vertical net samples was performed from 50 mL aliquots settled for 24 h and observed thoroughly under an inverted microscope (Nikon Eclipse TE2000-S).
Cell abundance estimation was performed from 50 mL of Van Dorn water samples settled in Utermöhl chambers (Utermöhl, 1958) during 24 h. Cell counting was performed under an inverted microscope (Nikon Eclipse TE2000-S), where the entire bottom of the chamber was examined at 200× magnification.
To enumerate the small and more abundant organisms, one or two transects were counted at 200× or five/ten fields randomly chosen at 400× until reaching 100 cells.
The Utermöhl method is widely used, it enables identification and quantification of phytoplankton samples (Intergovernmental Oceanographic Commission [IOC], 2010), and it is recommended for phytoplankton samples with low abundance (Hallegraeff, 2003). The method has been standardized (CEN, 2006) and it is frequently used in research and monitoring programs for quantitative phytoplankton analysis (e.g., Reguera et al., 2016;Wasmund, 2017).
Species identification was based on Steidinger and Jangen (1997), Hoppenrath et al. (2009Hoppenrath et al. ( , 2014, Omura et al. (2012), and Lassus et al. (2016). The validity of names of the different taxa was checked on the World Register of Marine Species (Horton et al., 2018), the list of toxic species was checked on the Taxonomic Reference List of Harmful Micro Algae (Moestrup et al., 2009). Those species identified only at the genus level, but with clear morphological differences, were listed as sp.1, sp.2, etc.

Data Analysis
Exploratory data analysis (EDA) was conducted to both, environmental and biological databases previous to any statistical analysis. The EDA allowed us to confirm general assumptions of the 'family' distribution of the variables (i.e., binomial distribution for presence-absence), the presence of homoscedasticity, the independence of the values (i.e., lack of spatial autocorrelation) and collinearity (i.e., correlation among explanatory variables), as suggested by Zuur et al. (2010). Before fitting regression models all the oceanographic variables which showed a Pearson correlation coefficient (rho) higher than 0.75 under the collinearity analysis, were dropped from the model (Supplementary Figure S1).
Dinoflagellate community analysis considered the calculation of the species richness over the net samples. From an ecological perspective we looked for discrete groups (clusters) of species composition over the sampling area. To achieve this, the Sorensen dissimilarity matrix for presence-absence (binary) data (Gower and Legendre, 1986) was obtained using the "dist.binary" function from the ade4 package (Dray and Dufour, 2007). Then, a hierarchical clustering was performed using the Ward's minimum variance method by the use of the "hclust" function and "method = ward.D2" argument in the stats package (R Core Team, 2018; version 3.4.4).
In order to evaluate the role of environmental variables over the community structure of the most common species observed in net samples (frequency of occurrence higher than 75%), a canonical correspondence analysis (CCA) (ter Braak and Verdonschot, 1995) through the "cca" function of the vegan package (Oksanen et al., 2018) was fitted. This multivariate constrained ordination technique allows us to include unimodal relationships inside the correspondence analysis since the response of species to abiotic variables typically follows an unimodal response (Greenacre and Primicerio, 2013). To evaluate the significance of the CCA model the "anova.cca" function (from the vegan package, Oksanen, op.cit) which performs an ANOVA like permutation test using 999 permutations (Borcard et al., 2018) was used.
From a species-specific perspective, the effects of environmental variables over the presence probability were assessed for three toxic species present in net samples. Dinophysis caudata was selected considering the magnitude of the damage caused by their blooms  and the high occurrence in water samples (see section "Harmful Species"). Karlodinium spp. was the potentially toxic genus observed in highest abundances in Van Dorn samples (see section "Harmful Species"). Finally, considering its benthic origin, Ostreopsis cf. ovata was studied as well.
Due to the binomial nature of the data in net samples (presence/absence) the Generalized Additive Logistic Models (GALoM) were fitted. The implementation of the "GALoM" models was done using the mgcv package through the "gam" function and the "family" argument defined as binomial (Wood, 2017). The model selection followed a backward selection method, based on the significance of each explanatory variable and using the Akaike's Information Criterion (AIC). This criterion negatively penalizes excess parameters, preventing from over-parameterization and allows for multimodel comparison where lower AIC values represent more parsimonious models.
All the statistical analysis and visualizations were conducted with the R software (R Core Team, 2018; version 3.4.4) and supervised by OneMind-DataScience 1 .

Hydrography and Environmental Variables
The spatial distribution of SST and sea surface salinity (SSS) (<5 m depth) was practically uniform in the area of study, with mean SST of 27.7 • C, and with mean SSS of 34.5. The vertical profiles were generally characterized by the presence of an upper mixed layer with a varying thickness delimited by a sharp thermocline, halocline and oxycline (Supplementary Figures S2-S4).
Below the mixed layer, temperature and dissolved oxygen declined with depth and salinity tended to increase (Supplementary Figures S2-S4). In this study, a subsurface salinity maximum was revealed in several profiles at different depths, showing peaks of salinity maxima (S > 35.5) in Santa Fé (transect A, B, and D) and Santa Cruz Island (transect, B), while the higher subsurface salinity value (S = 36.5) was registered in Santa Fé Island, transect A, 5 nm, at about 13 m depth (Supplementary Figure S2).
A slight increase in dissolved oxygen concentrations was observed at depths where such subsurface maximum of salinity was present (Supplementary Figure S4). The lowest salinity value (S < 34.0) belongs to Santa Cruz (B, 10 nm from the coast) at surface (0 m depth). The mean surface pH value was homogeneous for the study area with values about 7.9. The pH profiles (Supplementary Figure S5) followed similar trends as those of temperature and dissolved oxygen, even though showing higher variability (Supplementary Figure S5).
Regarding cell abundance estimation from Van Dorn bottles samples, no blooms were detected but higher abundances were registered in Seymour Island, especially in transect B (2 nm surface), with maximum values registered for Heterocapsa spp.
Dinoflagellate species richness varied among stations and ranged from a maximum of 53 to a minimum of 23 taxa per station (Figure 2). The lowest number of species was recorded at Santa Fé Island showing a maximum of 29 and a minimum of 23 taxa per station.
From the community structure analysis, a total of five discrete community composition groups were found. These assemblages showed a contrasting spatial distribution (Figure 3) with groups 4 and 5 being more representative of the coastal areas of Pinzón, Santa Cruz, and Seymour.
From the CCA analysis (based in the most frequent species), the total inertia (i.e., total variance in species distributions) was 0.49, representing the 100% of total inertia in the model. From this, the total variance explained by the environmental variables (i.e., constrained inertia) was low (8.4%) even though significant (F = 1.346, p-value < 0.05). The first two axes explained 7.2% of the constrained inertia with a total of 4.7% and 2.5% being explained by the first and second axes, respectively (Figure 4).
All the species represented in the CCA diagram, showed a homogenous distribution across the environmental variables, with some exceptions. A negative association with temperature was found for Protoperidinium conicum, Dinophysis caudata, and Prorocentrum compressum. In contrast, a positive relationship with temperature was found for Protoperidinium pellucidum, Karlodinium spp., and Tripos fusus; even though Karlodinium spp. and T. fusus are closer to a unimodal response with temperature (Figure 4). In relation with salinity, the species that showed a positive relationship were Podolampas bipes, Ceratocorys horrida, Podolampas palmipes, and Histioneis costata. However, a negative relationship was found for the following species: Tripos lineatus, Tripos candelabrum, Cucumeridinium coeruleum, and Ostreopsis spp. For Ostreopsis spp. this pattern is consistent with the GALoM model results (next section, Figure 5). No clear association between the species composition and pH were found.

Harmful Species
Nineteen potentially harmful taxa were identified in net samples ( Table 1). The most frequent were Gonyaulax spinifera, that was present in all samples, Phalacroma mitra, which was present in 63% and followed by Dinophysis caudata that was found in 46% of the samples. Benthic species were also present in the water   Table 1).
For Karlodinium spp. the model explained 17.3% of the deviance and the presence probability was significantly correlated with all the variables. The effects of temperature over presence probability showed a unimodal response, with the lower probability at about 26.5 • C, thereafter the probability started to increase (χ 2 = 238.2, p-value < 0.001). In the case of salinity, the response showed a negative linear decrease with higher values of salinity (χ 2 = 157.6, p-value < 0.001). The effects of pH showed a positive non-lineal relationship (χ 2 = 192.1, p-value < 0.001; Figure 5B).
In the case of Ostreopsis cf. ovata the model explained 13% of the total deviance and the presence probability diminished significantly with temperature and salinity (χ 2 = 103.1, p-value < 0.001 and χ 2 = 183.2, p-value < 0.001, respectively). Additionally, the presence probability showed a unimodal response with pH (χ 2 = 104.7, p-value < 0.001), with the highest values found at pH about 7.9 ( Figure 5C).

Dinoflagellates Assemblages
The number of taxa found in the present study is almost three times higher than the maximum observed in all the previous studies (Marshall, 1972;Tapia, 2000, 2002;Torres, 2002;Tapia and Naranjo, 2012;Torres and Andrade, 2014;Naranjo and Tapia, 2015), with 102 dinoflagellate taxa detected for the first time in the GMR. Some of the specimens, in this study, were only identified to genus level. Identification of armored dinoflagellates can be improved by the use of fluorescent stains such as calcofluor, which is an aid to examine the thecal plate pattern. Another method to improve identification consists in adding a drop of a diluted (5%) solution of sodium hypochlorite (bleach) applied to the edge of the coverslip which can be used for the separation of the dinoflagellate thecal plates (Tomas, 1997). These methodologies were not applied to the phytoplankton samples collected in the GMR in the present study. The use of these methods would probably have contributed to a deeper taxonomic identification of some of the specimens resulting in a longer list of species identified.
In addition, this study constitutes the first approach to the study of dinoflagellate assemblages in the GMR, since former studies have focused only on identification and estimation of cell abundance. The study was performed during a precise temporal sampling period corresponding to the wet season (March-April) during a weak La Niña event. Further discussion will be limited to the dinoflagellate community assemblages, its spatial variability and the response of specific species to the environmental conditions under this temporal span.
One of the first studies performed in the archipelago was done by Marshall (1972), covering 25 stations at seven depths in August 1968, during a moderate El Niño event. A total of 26 dinoflagellate taxa were reported at very low concentrations with a maximum of 326 cell L −1 in a station located north of Santa Cruz Island. Similarly, in August 2000, during the dry season and under the influence of a weak La Niña event, Torres and Tapia (2002) found 36 species of dinoflagellates during an expedition in the GMR and reported a high Chl a concentration north of Santa Cruz Island. This area has been previously defined as one of the most productive habitats in the GMR (Schaeffer et al., 2008), probably supported by topographic upwellings of the EUC and by the island-derived iron enrichment due to the IME (Palacios, 2002). This spatial abundance distribution coincides with the present study, where not only higher dinoflagellate abundances were observed in northern sites (North Seymour and Santa Cruz), but also the dinoflagellate richest stations (up to 53 taxa). Furthermore, Santa Cruz Island is considerably larger in size than the other islands. Biogeographically, there can be influenced differently either by the nutrient discharge or composition, depending on currents and winds. Those factors should be taken into account in future investigations in the area to study their contribution to the dynamics of dinoflagellate communities.
In addition, in Santa Fé Island, which is located of south of Santa Cruz Island, the lowest dinoflagellate richness was observed. Moreover, a particular community assemblage (group 1 in Figure 3B), not observed in the rest of islands, was present in this area. This particular dinoflagellate composition found at Santa Fé Island could be related with special features of this area, characterized by the presence of high-salinity waters likely associated with the cooler, carbon-nutrient rich EUC (Supplementary Figures S2, S3). Regarding specific species, Protoperidinium conicum was linked with Santa Fé sites, being isolated from the rest of the species in the CCA analysis and related to cold waters (Figure 4). However, further investigation should be addressed to this area in order to identify potential species associations to different environmental conditions, representative of the entire niche occupied by the species, since the present study includes only one seasonal period under particular climatic perturbations like the presence of a weak la Niña event.
Reports performed by the INOCAR highlight consistent diatom dominance over other groups in phytoplankton communities in the GMR. Chávez et al. (1996) reported low phytoplankton biomass in the ETP, where more precisely, open ocean waters were dominated by small and solitary organisms such as dinoflagellates, while diatoms were associated with upwelling and the Equatorial Front. Moreover, during La Niña in May 2007, McCulloch (2011 reported a dominance of pelagophytes, haptophytes, euglenophytes, and chrysophytes relative to Chl a in the phytoplankton community of the archipelago. Thus, it is probable that these groups were dominant over dinoflagellates in the present study, which was conducted during a weak La Niña, explaining in part, the low dinoflagellates abundances. There is scientific evidence that ENSO events explain, to a certain extent, Chl a concentration along the GMR (Kislik et al., 2017). During La Niña events, a shallower thermocline increases macronutrients supply (Chávez et al., 1996), however, it has not been linked with an increase in phytoplankton abundances in the ETP (Strutton et al., 2008). The study conducted by Kislik et al. (2017) pointed out that from 2003 to 2016, only one high Chl a concentration out of six, coincided with a La Niña event in the GMR. During this period of time, the variability in Chl a concentration found across the GMR occurred a few months later than the SST highest anomalies were registered. It is thus possible that Chl interannual variability could be partly explained by annual post-ENSO events. As the present study occurred during a weak La Niña event, the low dinoflagellate abundances registered could be expected based on this statement. In addition, during the wet season, warmer waters from the Panama Current flow southward, southeast trade winds weaken, ITCZ migrates to the south and rainfall increases (Palacios, 2002;Edgar et al., 2004;Kessler, 2006;McCulloch, 2011;Kislik et al., 2017). Under those conditions, a low Chl a level has been observed, more pronounced in the far northern Islands, Darwin and Wolf (Kislik et al., 2017).
Moreover, changes in the composition and distribution of phytoplankton are not only due to seasonal sub-surface variations driven by large marine currents, but also due to the influence of local winds and local currents, different tidal inflows, mixing and sediment resuspension the long-distance larval transport capabilities of species, among other factors (Chávez et al., 1990;Lafabrie et al., 2013;Venrick, 2015;Anabalón et al., 2016;Jacox et al., 2016;Zhao et al., 2018). This makes the characterization and evaluation of the phytoplankton distribution a more complex process.

Harmful Algae
In this study, we registered the highest number of harmful dinoflagellates ever recorded in the GMR, with a total of 19 species ( Table 1). The only other study focused on HABs in Ecuador reported proliferations of two species of diatoms, Mesodinium rubrum and Bellerochea malleus in the GMR (Torres, 2015). In addition, a bloom of Prorocentrum gracile (reported as P. gracilis) was associated with a fish mortality event in the archipelago in 1980, however, in the present study, this species was not observed.
Furthermore, Torres and Andrade (2014) in a study conducted during the dry season in September 2006, in Aeolián Bay (Baltra Island, south of Seymour Island), registered harmful species such as Prorocentrum cf. mexicanum and Ostreopsis cf. siamensis. The morphological identification of Ostreopsis species is very difficult since they overlap in size and have the same tear-drop shape (Penna et al., 2005). In the present study, identification of O. cf. lenticularis and O. cf. ovata is supported by molecular phylogeny and SEM images from epiphytic samples taken during the same period (Carnicer et al., in preparation). The presence of benthic species was also reported in a phylogenetic study on Prorocentrum lima (Nagahama et al., 2011) which included strains from Santa Cruz Island; unfortunately, the location was not specified. The present study also reported three other benthic taxa in the water column, O. cf. lenticularis, P. lima, and Coolia spp., with the dominant species being O. cf. ovata (46% occurrence). Their presence may be in response to blooms that occur on the substrate in the surrounding areas (Carnicer et al., 2015). However, the knowledge of epiphytic dinoflagellates assemblages and their blooming behavior along the GMR is not known.
Blooms of O. cf. ovata occur when a threshold temperature (25 • C) is achieved during summer in the Mediterranean Sea, probably linked to cyst germination (Accoroni et al., 2015). In the Galápagos, the SST variability over the year (Schaeffer et al., 2008) is similar to temperate areas. During 2017, in Charles Darwin Research Station located in Academy Bay on Santa Cruz Island, a difference of 8.2 • C was registered between April (27.3 • C) and August (19.1 • C). Thus, in the present study, this temperature variability may have contributed to O. cf. ovata bloom formations during the warm season. The environmental parameters used for the calibration of the GALoM model have to be taken with caution considering the benthic origin of this species. However, the statistical analysis showed that Ostreopsis spp. has a negative relationship with both, temperature and salinity. Other studies have found similar results, for instance those studies conducted in the Florida Keys and Hawaiian Islands (Parsons et al., 2012 and references therein). Few studies exist about the influence of pH over Ostreopsis spp. Di Cioccio et al. (2014) through their study about the possible effects of pH decrease on benthic HABs pointed out that O. cf. ovata seems to be tolerant to a wide range of pH values. In the present study, the pH values showed a narrow variability with the response of Ostreopsis spp. being unimodal with the highest occurrence probability at a pH about 7.9. Because the already cited work of Di Cioccio et al. (2014) was carried out in a much more dynamic scenario with respect to pH [stations located along marked pH gradients (6.8-8.1)], we can no longer compare and/or discuss their findings.
In the present study, the first record of Karlodinium spp. in the GMR is reported. There are at least 12 species described within the genus Karlodinium (Luo et al., 2018). Some species have been associated to fish kills: K. armiger (Garcés et al., 2006), K. australe (Lim et al., 2014), K. conicum (de Salas et al., 2008), K. corsicum (Paulmier et al., 1995), K. gentienii (Nézan et al., 2014), and K. veneficum (Place et al., 2012). The water column stratification has been suggested as the main factor that favors blooms of Karlodinium species such as K. veneficum (Place et al., 2012), as well as the abundance of preys (Lin et al., 2018). The stratification of the water column due to warm conditions and rising precipitation, a property of wet season, could benefit the presence of Karlodinium in the GMR area. Relative to the pH variable, the statistical analysis showed that the probability of presence of Karlodinium spp. rises with pH > 7.8. A study of a bloom of K. australe conducted in west Johor Strait, Malaysia, by Lim et al. (2014) showed that the pH values found during the period of bloom in the strait was also >7.8. Those authors theorize the existence of a relation between the outbreak of K. australe blooms and the influx of nutrients mainly from anthropogenic activities from nearby areas. Under this scenario, further studies are required to shed light on the positive response of Karlodinium spp. to some pH values in order to known if this response corresponds to anthropogenic or natural forces. The abundance of Karlodinium spp. detected in the present study (5 × 10 4 cell L −1 ) was low in comparison to the abundances that can be reached during bloom events in different areas of the world (Place et al., 2012;Lim et al., 2014). The fact that this study was conducted during a weak La Niña event could have contributed to a lower probability of blooms of Karlodinium spp. due to the stratification conditions being reduced and SST and precipitation being lower.
Dinophysis caudata is a neritic species associated to production of DSP toxins . In this study, D. caudata presented a 46% of occurrence (Table 1), being one of the most common species in the GMR. The presence of this species was associated with temperate and tropical waters, as previously cited by Kofoid and Skogsberg (1928). A negative association with higher pH and temperature values, together with the increase in the probability of presence of D. caudata to higher salinities (Figure 5A), could indicate local upwellings associated to the presence of D. caudata. The presence of D. caudata in the ETP is limited to Mexico, where a monitoring program has been running since the 1990s (Daguer et al., 2018). This species has been reported in the Gulf of California (Esqueda-Lara et al., 2013) and in Chiapas (Hernandez-Becerril et al., 2008;Maciel-Baltazar, 2015), but in very low abundances . Thus, its limited distribution and its general low abundances could be the reason why it was not observed in the quantitative sampling (Van Dorn bottles).
The evaluation of the dinoflagellates-climate link provides baseline information to forecast biological responses under future environmental changes in the oceans (McCulloch, 2011). Hence, the study of dinoflagellate assemblages presented in this research constitute an advance in the knowledge about their diversity, abundance and spatial variability during the wet season in a weak La Niña event around Santa Cruz Island and nearby islands. This study provides original information about phytoplankton communities that will be useful for future work as a reference, not only to identify follow-up studies on potentially toxic species, but also to identify possible new species introduced by human activities and that may contribute to damage of the marine ecosystems of the GMR (Tian et al., 2017).

CONCLUSION
The dinoflagellate assemblage found in this study represents the highest diversity observed until now in the GMR, with 152 taxa identified. However, abundances were low and no blooms were detected, in agreement with similar abundance ranges previously reported. This is probably related with the environmental conditions occurring at the time of sampling (during the wet season and under a weak La Niña event). Although this work is from a restricted temporal scale, this study proposed and demonstrated a systematic data analysis that was able to reveal discrete dinoflagellate assemblages and among-island variability of environmental conditions. Moreover, a link between those assemblages and some particular dinoflagellate species was established with the environmental conditions prevailing during the sampling period.
Dinoflagellate species richness found in the north of Santa Cruz Island, together with higher abundances could be associated to the IME. A different dinoflagellate assemblage was distinguished in the surrounding waters of Santa Fé Island correlated with cooler and high-salinity waters possibly connected with the EUC.
A long-term monitoring program needs to be set up following an internationally validated protocol to assure accurate results considering the high ecological value of the GMR. Furthermore, the government must invest, both in sophisticated equipment and taxonomy training, in order to carry out a more extensive identification up to species level to avoid possible data misinterpretation.
Because this methodology is mainly focused in the response of dinoflagellate community (and particular species) to different environmental conditions, this approach can be easily extrapolable to different areas. However, this extrapolation should be taken carefully and must include the whole analytical procedure, including the spatio-temporal (when available) variability of environmental conditions and species presence (or abundance) in the new area, as well as the whole model calibration and variable selection procedure.
Efforts should also be addressed to the characterization of potentially toxic strains collected in the region considering the intra-specific variation among strains with different geographic origin. A parallel analysis of toxin content from water samples would further contribute to estimations of the potential HAB risk in the GMR.

AUTHOR CONTRIBUTIONS
OC and MF-T made substantial contributions to conception and design of the work. OC and IK conducted the sampling. ER-M and MF-T contributed to data collection. OC, MF-T, PD, and AC made substantial contributions in data analysis and interpretation, and participated in drafting the article. JD revised the article critically. OC, PD, AC, IK, ER-M, JD, and MF-T gave the final approval of the version to be submitted. FUNDING a color scale (red and blue represents positive and negative Pearson's correlation coefficients, respectively). Additionally, the value of each coefficient is given inside each comparison box.
FIGURE S2 | Vertical profiles for temperature ( • C). The information for each island and transect is showed in columns and rows, respectively. Station is color code: blue color for 2 nm, green color for 5 nm, and yellow color for 10 nm.
FIGURE S3 | Vertical profiles for salinity. The information for each island and transect is showed in columns and rows, respectively. Station is color code: blue color for 2 nm, green color for 5 nm, and yellow color for 10 nm.
FIGURE S4 | Vertical profiles for dissolved oxygen (mg L −1 ). The information for each island and transect is showed in columns and rows, respectively. Station is color code: blue color for 2 nm, green color for 5 nm and yellow color for 10 nm.
FIGURE S5 | Vertical profiles for pH. The information for each island and transect is showed in columns and rows, respectively. Station is color code: blue color for 2 nm, green color for 5 nm, and yellow color for 10 nm.