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 (sub-thermocline) compensation against the westward SEC (Schaeffer et al., 2008). The collision of the subsurface iron and nutrient-carbon 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).

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 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%).

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.

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. (2009, 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.

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 (Reguera et al., 2014) 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 multi-model 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¹.

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).

A total of 152 dinoflagellate taxa were identified in net samples, belonging to 7 orders, 22 families, and 38 genera. The most diverse genus was Tripos (25 species), followed by Protoperidinium (22 species), Dinophysis (13 species), Oxytoxum (9 species), and Phalacroma (7 species) (Supplementary Table S1). The taxa described as Gonyaulax spp., Scrippsiella spp. and Heterocapsa spp. were found in all stations, as well as Gonyaulax spinifera. Only 12 taxa were found in 75% of the samples, whereas 100 species appeared in 25% of them. Species with more occurrence were Tripos furca (88%), Tripos muelleri (79%), Oxytoxum tesselatum (85%), Podolampas bipes (85%), and Protoperidinium steinii (83%).

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. (5,832 cell L^-1) and Protoperidinium pellucidum (4,000 cell L^-1) (Supplementary Table S2).

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.

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 column: Ostreopsis cf. ovata (46%), Ostreopsis cf. lenticularis (4%), Prorocentrum lima (17%), and Coolia spp. (17%) (Table 1).

No HABs were detected in Van Dorn bottle samples. The highest abundance corresponded to an ichthyotoxic genus, Karlodinium spp. (3,940 cell L^-1) that was registered in transect A (2 nm, 10 m depth) at Seymour Island. Other harmful species were recorded in very low abundances (<280 cell L^-1), such as Gonyaulax spinifera, Porocentrum micans, P. lima, P. cordatum, P. mexicanum, Karenia spp., Karenia papilionacea, and Alexandrium spp. (Supplementary Table S2).

For Dinophysis caudata the GALoM model explained 7.94% of the deviance and the presence probability was significantly correlated with all the variables, showing a clear negative association with temperature (χ² = 106.2, p-value < 0.001) and a non-lineal positive response with salinity and pH (χ² = 60.9, p-value < 0.001 and χ² = 122.2, p-value < 0.001, respectively; Figure 5A).

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 (χ² = 238.2, p-value < 0.001). In the case of salinity, the response showed a negative linear decrease with higher values of salinity (χ² = 157.6, p-value < 0.001). The effects of pH showed a positive non-lineal relationship (χ² = 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 (χ² = 103.1, p-value < 0.001 and χ² = 183.2, p-value < 0.001, respectively). Additionally, the presence probability showed a unimodal response with pH (χ² = 104.7, p-value < 0.001), with the highest values found at pH about 7.9 (Figure 5C).

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; Torres and 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.

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⁴ 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 (Reguera et al., 2014). 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 (Reguera et al., 2014). 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).