Abundance and Structure of the Zooplankton Community During a Post-eruptive Process: The Case of the Submarine Volcano Tagoro (El Hierro; Canary Islands), 2013-2018

The mesozooplankton community was analyzed over a 6-year period (2013-2018) during the post-eruptive stage of the submarine volcano Tagoro, located south of the island of El Hierro (Canary Archipelago, Spain). Nine cruises from March 2013 to March 2018 were carried out in two different seasons, spring (March-April) and autumn (October). A high-resolution study was carried out across the main cones of Tagoro volcano, as well as a large number of reference stations surrounding El Hierro (unaffected by the volcano). The zooplankton community at the reference stations showed a high similarity with more than 85% of the variation in abundance and composition attributable to seasonal differences. Moreover, our data showed an increase in zooplankton abundance in waters affected by the volcano with a higher presence of non-calanoid copepods and a decline in the diversity of the copepod community, indicating that volcanic inputs have a significant effect on these organisms. Fourteen different zooplankton groups were found but copepods were dominant (79%) with 59 genera and 170 species identified. Despite the high species number, less than 30 presented a larger abundance than 1%. Oncaea and Clausocalanus were the most abundant genera followed by Oithona and Paracalanus (60%). Nine species dominated (>2%): O. media, O. plumifera, and O. setigera among the non-calanoids and M. clausi, P. nanus, P. parvus, C. furcatus, C. arcuicornis, and N. minor among the calanoids. After the initial low abundance of the copepods as a consequence of the eruption, an increase was observed in the last years of the study, where besides the small Paracalanus and Clausocalanus, the Cyclopoids seem to have a good adaptive strategy to the new water conditions. The increase in zooplankton abundance and the decline in the copepod diversity in the area affected by the volcano indicate that important changes in the composition of the zooplankton community have occurred. The effect of the volcanic emissions on the different copepods was more evident in spring when the water was cooler and the mixing layer was deeper. Further and longer research is recommended to monitor the zooplankton community in the natural laboratory of the Tagoro submarine volcano.

The mesozooplankton community was analyzed over a 6-year period (2013-2018) during the post-eruptive stage of the submarine volcano Tagoro, located south of the island of El Hierro (Canary Archipelago, Spain). Nine cruises from March 2013 to March 2018 were carried out in two different seasons, spring (March-April) and autumn (October). A high-resolution study was carried out across the main cones of Tagoro volcano, as well as a large number of reference stations surrounding El Hierro (unaffected by the volcano). The zooplankton community at the reference stations showed a high similarity with more than 85% of the variation in abundance and composition attributable to seasonal differences. Moreover, our data showed an increase in zooplankton abundance in waters affected by the volcano with a higher presence of non-calanoid copepods and a decline in the diversity of the copepod community, indicating that volcanic inputs have a significant effect on these organisms. Fourteen different zooplankton groups were found but copepods were dominant (79%) with 59 genera and 170 species identified. Despite the high species number, less than 30 presented a larger abundance than 1%. Oncaea and Clausocalanus were the most abundant genera followed by Oithona and Paracalanus (60%). Nine species dominated (>2%): O. media, O. plumifera, and O. setigera among the non-calanoids and M. clausi, P. nanus, P. parvus, C. furcatus, C. arcuicornis, and N. minor among the calanoids. After the initial low abundance of the copepods as a consequence of the eruption, an increase was observed in the last years of the study, where besides the small Paracalanus and Clausocalanus, the Cyclopoids seem to have a good adaptive strategy to the new water conditions. The increase in zooplankton abundance and the decline in the copepod diversity in the area affected by the volcano indicate that important changes in the composition of the zooplankton community have occurred. The effect

INTRODUCTION
Zooplankton plays a key role in the functioning of marine food webs that control the biogeochemical ocean cycles (Longhurst, 1985) where usually the copepods are dominant (Longhurst, 1995;Vinogradov, 1997). They are the major prey of pelagic fauna and a key component of the biological pump, transporting organic matter to the deep ocean throughout their diel or seasonal migration (Longhurst, 1995;Hernández-León and Ikeda, 2005). Zooplankton population structure can be affected by the productivity of the marine area, as well as the hydrographic conditions where they live (Beaugrand and Ibanez, 2004). Thus, analyzing the zooplankton community and its species over time will provide valuable information, essential for diversity studies and the efficient management of the marine pelagic ecosystems, in particular the most vulnerable. In addition, identifying species morphologically is a challenging task. Consequently, there are fewer studies in this area of expertise because such studies require considerable time, effort, and years of specialization.
Marine areas subjected to volcanic activity are very vulnerable and unfortunately poorly analyzed in comparison to other ecosystems. They can be used as a template for studies of global climate change to predict and understand global scenarios of future oceans (Dahms et al., 2018). The Canary archipelago is located in the north-eastern part of the Subtropical North Atlantic Ocean as part of an intraplate hotspot chain offshore NW Africa, with decreasing ages of seamounts and volcanic islands from east to west (Carracedo et al., 1998;Klügel et al., 2020). El Hierro, situated on ∼155 Ma old ocean crust, is the youngest island along the chain and the island with the most recent submarine volcanic activity in the last 500 years of geological history (Tagoro submarine volcano -2011; Klügel et al., 2020). The marine ecosystem of El Hierro before 2011 was not linked to a recent hydrothermal process. El Hierro is also the westernmost island of this archipelago (Figure 1), open to the oligotrophic waters of the North Atlantic Subtropical Gyre where nutrient concentrations are usually undetectable in the epipelagic water column (Santana-Casiano et al., 2013). Due to the offshore position from the NW African upwelling, El Hierro borders some of the clearest waters surrounding the Canary Islands (Barton et al., 1998), and the epipelagic zone is characterized by low mesozooplankton biomass and well oxygenated waters (Arranz, 2007). Despite the low productivity of El Hierro waters, fisheries resources are abundant, particularly in the Mar de las Calmas marine protected area, where only artisanal fishing is allowed (Tuya et al., 2006).
Submarine volcanoes constitute a significant source of gases, solutes, and heat into the ocean, discharging hot and acidic fluids, reducing oxygen concentrations, and exerting significant controls on the composition of seawater and the marine biota (Sedwick and Stuben, 1996). In October 2011, after more than 15000 earthquakes, a volcanic tremor was registered in the south of the island of El Hierro, indicating the beginning of a submarine eruption. The submarine volcano Tagoro was located 1.8 km south of El Hierro, and the main cone grew from an initial depth of 355 m to a maximum height of 88 m below the sea surface (Fraile-Nuez et al., 2012Sotomayor-García et al., 2019). Tagoro emitted molten material during the first 6 months, affecting a wide area around the island, with different peaks of intensity. During the eruptive stage (October 2011-March 2012), dead fish were observed floating on the surface waters (Fraile-Nuez et al., 2012). The emission of magma, gases, and volcanic particles drastically changed the physical-chemical properties of the water column, producing temperature and salinity anomalies of +3 • C and −0.3 respectively, with a maximum temperature anomaly of +18.8 • C (Fraile-Nuez et al., 2012). Furthermore, the intense bubbling and degassing caused a marked increase of CO 2 and a decrease of pH in more than 2.9 units (Santana-Casiano et al., 2016). Large amounts of metals and inorganic nutrients were emitted (Fraile-Nuez et al., 2012;González-Vega et al., 2020), as well as other compounds (Santana-Casiano et al., 2013), contributing to the decrease of dissolved oxygen, occasionally reaching anoxic levels. These conditions led to a drastic reduction of the marine biota both in the benthic and pelagic ecosystem (Santana-Casiano et al., 2013).
Following the eruptive stage (March 2012), the Tagoro submarine volcano entered in a hydrothermal moderate phase involving the release of gases, heat, metals, and nutrients from multiple vents dispersed around the main cone (Santana-Casiano et al., 2016;Fraile-Nuez et al., 2018;Sotomayor-García et al., 2019). Currently, the Tagoro volcano exhibits new benthic habitats thriving around the main and secondary craters, colonized by small hydrozoan colonies with a high diversity of annelids, arthropods, cnidarians, and mollusks as the first colonizers (Sotomayor-García et al., 2019). Although in other submarine volcanic areas the hydrothermal vent emissions have been observed to have a considerable effect on the diversity, abundance, biology, and ecology of the planktonic communities (Tarasov, 2006), no clear impact has yet been proved on the local phytoplankton community at Tagoro volcano (Gómez-Letona et al., 2018). Monitoring of the new submarine volcano could provide a unique opportunity to study the resilience of the pelagic ecosystem in El Hierro after an acute volcanic disturbance as well as to depict the effect of the subsequent diffuse volcanic emissions on the zooplankton community. (C) High-resolution bathymetry of the volcanic edifice and high-resolution transect ("subnet;" stations 50-59 and 61) carried out across the main cones. All maps were generated using Matlab 8.5 R2020a.
Preliminary analysis of main zooplankton groups showed similar composition in those stations unaffected by the volcano and around El Hierro island (de Vera et al., 2017;Fernandez de Puelles et al., 2017). During the eruptive stage and close to the volcano, a clear shallowing of the deep acoustic scattering layer was observed with a disruption of the diel vertical migration of mesopelagic organisms, which may affect the distribution of the pelagic communities (Ariza et al., 2014). Besides these studies, no previous reports exist on the structure and composition of the zooplankton communities in the area or about their changes in response to volcanic activity.
Due to the difficulty of sampling the area after the drastic volcanic eruption, which completely covered the surface layers of the water column with ash, a zooplankton monitoring effort started in March 2013 and continued to March 2018. During these 6 years, the zooplankton community was analyzed in the epipelagic waters around El Hierro (Figure 1) and more intensively in a high-resolution transect across the main cones of the Tagoro submarine volcano ( Figure 1B). Because of the lack of zooplankton community data in the subtropical area, our main goals were (i) to improve our knowledge regarding biomass, abundance, and species composition with special emphasis on its dominant group, the copepods to completely describe the seasonal variation of the pelagic zooplankton inhabiting this area and (ii) to test the hypothesis that post-eruptive volcanic emissions affect the zooplankton community.

MATERIALS AND METHODS
The zooplankton community was sampled between 2013 and 2018 through nine oceanographic cruises, almost twice per year during two different seasons: (i) Spring (March-April), with five oceanographic surveys (2013, 2014, 2016, 2017, and 2018) and (ii) Autumn (October), with four oceanographic surveys (2013, 2015, 2016, and 2017). During 2013-2014, the sampling area was not only focused around the submarine volcano but also around the whole island of El Hierro ( Figure 1A). After 2015, the sampling area was reduced to provide a higher resolution around the submarine volcano (1 km radius) with the implementation of a transect across the main volcanic cones (hereafter, subnet) ( Figure 1B) consisting of stations located a few meters apart, a feat achievable using the dynamic positioning system of the vessels. In addition, several unaffected stations were sampled as reference stations during each cruise outside the influence of the volcano.

Sample Collection
Vertical profiles of temperature, conductivity, and pressure from the surface to 200 m depth were collected using a Seabird 911 plus CTD equipped with dual temperature and conductivity sensors, with accuracies of 0.001 • C and 0.0003 S/m respectively, continuously recording data with a sampling interval of 24 Hz. The CTD also carried sensors for fluorescence (Wet Labs ECO-Fl) and dissolved oxygen (Seabird 43). All CTD sensors were calibrated at the Seabird laboratory before and after the cruises. Discrete water samples of inorganic nutrients were obtained using a Rosette equipped with 24, 12-liter Niskin bottles. The dynamic positioning (DP) system of the research vessels allowed precise positioning of the stations separated by only a few meters of distance. Moreover, for continuously sampling the affected area close to the seabed over the main craters, several towyo transects were carried out following the methodology of our research team in the area (Santana-Casiano et al., 2016) which consist in continuously lowering and raising the rosette between 1 and 40 m above the seabed and with the ship moving at 0.2-0.4 kn in DP, obtaining a sawtooth-shaped dataset with a high spatial resolution close to the source of the emissions (González-Vega et al., 2020).

Zooplankton Sampling and Statistical Treatment
For the zooplankton analysis, a double WP2 net (56 cm frame) equipped with a 200 µm mesh was used to carry out a total of 160 tows (160 × 2 samples) in the study area. Vertical hauls were done from the surface to 200 m depth, or 5 m above the seafloor in those stations where the total depth of the water column was shallower than 200 m. The volume of water sampled was estimated as that of a cylinder where the height corresponds to the maximum depth of the corresponding haul and the base equals the area of the net opening. One of the nets was used for taxonomic studies and preserved with formaldehyde (4%) and the data given in individuals per cubic meter (ind.m −3 ). The other net was used for zooplankton biomass estimation using traditional methods (UNESCO, 1997), drying the sample to 60 • C for 24-48 h and the data given in mg of dry weight m −3 . During the March 2016 cruise, biomass samples were not collected.
In the laboratory, a Folsom plankton splitter was used to analyze the sample for taxonomic studies and at least two aliquots were counted, representing the total organism abundance. On each sample at least 500 ind.m −3 were counted. All zooplankton groups were identified and standardized to ind.m −3 . Copepods adults were identified, to species level and the juveniles to genus, following the most relevant taxonomic bibliography (Boltovskoy, 1981;Shmeleva, 2006, 2010;Razouls et al., 2020). Due to the complexity of the small size of Oncaea media, juveniles and adults were together named as O. media group. The Shannon-Wiener index was calculated to represent the copepod diversity (Shannon and Wiener, 1963).
Cluster and non-metric multidimensional scaling (NMDS) analyses were used to examine patterns in the community structure. The Bray-Curtis similarity index was applied coupled with group-average linkage. The same methodology was applied to the copepod species composition data to define copepod species assemblages involved in the study and the two sampled seasons. The similarity percentage (SIMPER) routine was then applied to identify the copepod species with higher contributions to the significant groups among pairs of samples selected by the factor "seasons" (spring and autumn) using Bray Curtis similarity resemblance after square root transformed abundance data. Significant differences in community structure between the two groups were tested by similarity analysis of one Anova way (ANOSIM). When the zooplankton data were related to the environment, they were firstly normalized. All procedures were performed using the Primer-7 software package for the above analyses from the Plymouth Laboratory (Clarke et al., 2014).

Other Statistical Treatment
Although the emissions during the post-eruptive stage of Tagoro volcano have been observed to reach an approximate diameter of 0.5 km around the main cone (Santana-Casiano et al., 2013;González-Vega et al., 2020), not all samples within this area necessarily show volcanic influence due to the intermittent behavior of the emissions as well as the effects of local currents (Fraile-Nuez et al., 2018). We used two physical-chemical parameters as tracers to determine the presence or absence of volcanic influence: temperature and concentration of silicate (González-Vega et al., 2020). The stations which presented anomalies in one or both of these parameters as compared to a reference profile were considered as affected. The reference profile was established for each cruise as the average of several stations outside of the affected area. Therefore, previous to the statistical treatments, all stations were classified as affected or not affected by the volcanic emissions.
Principal Component Analysis (PCA) was conducted to explore correlation patterns of environmental variables considering those stations affected and unaffected by volcanic influence; using mean values of temperature 0-200 m (T), salinity 0-200 m (S), dissolved oxygen 0-200 m (DO) and fluorescence 0-100 m (F) on each sampling station conducted during early spring vs. autumn season. Redundancy Analyses (RDA; Borcard et al., 2011) were also used to test the potential effects of the environment (T, S, DO, and F) on the dominant copepod species (>1.5% of the total abundance). The potential variance within sampling stations (consistently monitored during each cruise) was controlled including this co-variable as a condition factor. The significant effect of each environmental variable was assessed using the permutation procedure implemented in the ANOVA function. The goodness of RDA fitted was ensured after testing the linear dependencies among explicative variables using variance inflation factors (VIF) obtaining values > 3 (Borcard et al., 2011). RDA was completed using the RDA function of the R software package vegan (Oksanen et al., 2014).
In addition, generalized linear mixed models (GLMMs. fitted using R lme4 library; Bates et al., 2015) were used to test for potential differences in diversity (H'; Shannon index) and species abundance among the different cruises over the 6 year study period. In this sense, response variables H' (boxcoxtransformed) and abundance (square root transformation) were individually tested related to the numerical co-variable year and the categorical co-variables season (spring vs. autumn) and volcanic-effect (affected vs. unaffected). Considering the potential variability within sampling stations, the GLMMs incorporated the Stations as a random factor. After fitting the models, the distributions of the residuals were checked for normality.
The effect of volcanic emissions on the structure of volcanic emissions was tested using the adonis2 permutational analysis of variance using distance matrices (Oksanen et al., 2014) on distance matrices constructed using the Bray-Curtis index from abundance tables normalized using the decostand function with permutations constrained to samples of the same cruise in the seasonal analysis.

Environmental Scenery
Vertical profiles of temperature, salinity, fluorescence, and dissolved oxygen in the upper 200 m depth were averaged for each Frontiers in Marine Science | www.frontiersin.org Data from tow-yo transects and in situ samples were used to observe differences between affected and not affected stations in the water column. This sampling strategy allowed us to focus on the close surroundings of the hydrothermal vents, where all physical-chemical parameters showed significant anomalies with respect to the reference stations calculated as the average of all unaffected stations, separated by season (spring and autumn; Figure 2). The temperature and salinity anomalies (Figures 2A,B) presented a different pattern in spring (mostly located below 130 m depth) than in autumn (60-120 m depth, not seen in salinity). The concentration of silicate was studied using a similar approach (Figure 2C). Anomalies of silicates were found above 140 m depth, with a peak around 127 m regardless of seasonality. These results confirm the adequacy of these two parameters (temperature and silicates) as tracers for the classification of stations affected/unaffected by volcanic emissions. Moreover, fluorescence and dissolved oxygen were evaluated in the water column with the same criteria (Figures 2D,E). Surprisingly, fluorescence showed significatly different anomalies in spring from 80 to 170 m depth and shallower and less intense during autumn from 40 to 120 m depth. Similarly, dissolved oxygen showed an increase in its concentration just during spring from 80 to 170 m depth.
Principal Component Analysis (PCA) explained 81% of the total variance in the first two axes, indicating clear differences between the two seasons (Figure 3). The first axis was mostly influenced by temperature and fluorescence, whereas the second axis was mostly influenced by salinity and dissolved oxygen. Samples clustered primarily by season with small differences between years, e.g., the warmest autumn in 2015 and the coldest spring with high fluorescence in 2014. This indicates that seasonal changes in environmental variables exert a stronger influence on the composition of the zooplankton community than volcanic emissions when considered jointly.

Abundance and Biomass Zooplankton Interannual Changes During the Study
For the abundance and taxonomy, 102 zooplankton samples during early spring and 58 in autumn were analyzed, corresponding to more than 225,000 individuals counted and identified. Almost 80% of these organisms were copepods followed by appendicularians (9%), chaetognaths (3%) and ostracods (2% ; Table 1). Another nine holoplanktonic groups were recorded at low abundances below 1% as well as four groups of meroplankton larvae with decapod larvae (<1%) being the most abundant. Zooplankton abundance was higher in spring (755 ± 365 ind.m −3 ) than in autumn (506 ± 299 ind.m −3 ). The overall contribution of copepods was similar in both seasons and always exceeding 75%. Other groups were usually more abundant in spring than in autumn. The seasonal pattern of zooplankton biomass remained stable throughout the years, displaying very low values never higher than 8 mg.m −3 . Values were higher in spring than in autumn, with no clear pattern regarding the effect of the volcano ( Figure 4A).
On the other hand, the abundance of zooplankton throughout the period was variable, ranging from the lowest value of 436 ind.m −3 in March 2013 to the highest value of 1541 ind.m −3 in March 2018 ( Table 2). During 2013, similar zooplankton abundance was observed in all stations, both in spring and autumn. In March 2014, higher zooplankton abundance was observed in the whole area, coinciding with the high fluorescence and cool temperature of that spring. In October 2015, the zooplankton abundance was low and closer to the initial years of the study. Slightly higher values were observed in spring 2016 than in autumn 2016. Nevertheless, the last two surveys (October 2017 and March 2018) displayed a high increase in the zooplankton abundance, particularly in the area affected by the volcano. This increase was mainly due to small copepods, particularly the non-calanoids (Figures 4B-D) accompanied by modest increases in the abundances of appendicularians and chaetognaths.
After averaging the zooplankton abundance of each year, the hierarchical clustering analysis revealed two main assemblages, showing high similarity among the first 4 years (80%) and

General Composition and Diversity of the Copepod Community Around El Hierro
Of the 59 different copepod genera found throughout zooplankton samples, four genera (Oncaea, Clausocalanus, Oithona, and Paracalanus) contributed 60% of the total copepod population, forming the main bulk of the copepod population. Calanus, Corycaeus, Temora, Mecynocera, Calocalanus, and Acartia followed them being the other genera lower than 1%.
Calanoid copepods showed higher abundances around 65% in the initial years of the study as compared to usually comprising a larger share of the community in spring (65-70%) but their relative contribution declined over the 6 years. Conversely, noncalanoid copepods, particularly in autumn, increased from about 30% to around 50% in the last years of the study. Similarly, the contribution of Oncaea, Oithona, and Corycaeus increased over the years (Figure 6). Other non-calanoids rarely exceeded 2% of the total community. Clausocalanus and Paracalanus were the dominant calanoid copepods showing great variability in their abundances throughout the years. Calanoids had their highest abundance in spring 2014, showing a decline in later years. Nevertheless, other copepods such as Acartia and Temora after 2016 had a higher contribution in the whole area than previous years (Figure 6).
Of the 174 species of copepods identified, only 53 had a contribution > 2%. O. media was the most abundant (15%), followed by juveniles of Oithona and Clausocalanus. The copepods most frequently found in the samples were Mecynocera clausi, Nannocalanus minor, Clausocalanus furcatus, Lucicutia flavicornis, and Mesocalanus tenuicornis which were always present although they never dominated. Interannual variability was observed in the dominant taxa throughout the study (see Table 2). The dominant group of Oncaea media was usually constant with a drastic increase in later years (almost twice of the previous years). Oithona copepodites were more abundant during autumn and while O. setigera widely increased during the study, other congeneric species decreased such as O. tenuis and O. plumifera. Copepodites of Corycaeus also were more abundant in the later years (with more than 12 species identified), as well as C. arcuicornis and C. furcatus. Among the Calanids, N. minor was dominant with a high contribution during the spring of 2014 when the fluorescence was maximum, declining later. Differences were also observed among other congeneric species, P. parvus and P. denudatus increased regularly, in opposition to P. nanus.
Closer observations were seen with A. danae which declined along the study, meanwhile, A. negligens increased in the last period as well as T. stylifera. On the contrary, others such as L. flavicornis, frequent in all samples, had similar abundance throughout the study.
The dominant copepods as well as their contribution at the different seasons during the six years can be observed in Table 3. The majority of the calanoids predominated in spring as well as many Clausocalanus and Paracalanus. However, C. furcatus and P. denudatus predominated in autumn as well as most of the Oithona (such as O. plumifera, O. setigera, and O. tenuis) and the Corycaeus group. Great differences were not observed for other copepods such as M. clausi or M. tenuicornis which predominated in both periods.
The fourth main assemblages of copepods were observed by hierarchical cluster analysis, Oncaea media group and Oithona (similarity of 90%; group a) with the highest abundance dominating in both periods, due particularly to their high increase at the end of the study. Following them, the other 15 species were also very abundant (group b; similarity of 78%) represented by P. parvus, P. denudatus, M. clausi,  were predominant in autumn (Figure 7). All of them formed the main taxa of the copepod community during the 6year study. The copepods abundance experimented a significant increase from the beginning of the study (GLMM p < 0.001). The diversity   Average abundance (Av. Abund.) and relative contribution (Contrib.%). Juveniles (j).
seems to decline, particularly after 2016 (Figure 8), however this trend is not significant. The seasonality had a clear effect on the diversity and abundance but with an inverse pattern. Meanwhile, the diversity was higher in autumn and lower during spring (p < 0.05) the abundance showed higher values during spring (p < 0.001; Figure 8). Although no significant differences were observed on abundance between affected and unaffected stations, the diversity was significatively lower in the stations affected by the Tagoro volcano (p = 0.067; Figure 8).

Environmental Influence on Copepod Assemblages
In our study, the RDA analysis accounted for 38% of the total variance when we considered all canonical axes. Of this 38%, 75.88% was explained by the first two axes (p < 0.001). Temperature, salinity, fluorescence, and dissolved oxygen were found to be significant explanatory factors (p < 0.001). When the abundance of the different copepod species were related to the environmental variables by RDA analysis, the bulk of the most frequent copepods were located mainly along the temperature axis (Figure 9). And the distribution of the species along the two axes may indicate their preference for several environmental conditions. O. media and A. danae, on the right side with a cold temperature and high salinity are typical of earlyspring. On the other side, F. gracilis, Lubbockia spp., P. denudatus, C. furcatus, Macrosetella spp. and O. setigera thrived in warmer waters with lower salinity, conditions more usually observed in autumn.
In addition, PCA analysis of the most frequent copepod genera (> 10%) considering all stations in spring and autumn, indicated that often the copepods affected by the volcanic emissions diverged from those found in the surrounding waters (Figure 10). All genera were considered as variables for the analysis regardless of whether they were dominant or not, except for those which only appeared in a small number of stations (<20%), whose presence was considered accidental. This analysis shows a clear differentiation between the two groups of stations (affected and unaffected) concerning their abundance, clearer in spring. The affected stations seem to present a stronger presence of some of the dominant genera such as Oncaea, Temora, Acartia, and Corycaeus, but also some less abundant genera such as Centropages. Indeed, permutational multivariate analysis of variance based on Bray-Curtis distance matrixes of main copepod assemblages revealed that not only did they diverge between spring and autumn seasons (p < 0.01) but also that significant differences were found between affected and unaffected waters (p < 0.01), although these differences were only significant in spring (p < 0.01) not in autumn (p > 0.05).

DISCUSSION
The zooplankton community was analyzed over 6 years, during the post-eruptive process of the submarine volcano Tagoro, 17 months after the sizeable eruption of the volcano in a marine area affected by its influence. At that time, magmatic emissions had ceased but the area was still affected by a moderate hydrothermal activity characterized by diffusive flows of CO 2 , reduced compounds, and inorganic nutrients including nitrate, phosphate, silicate, and iron (Fraile-Nuez et al., 2012;Santana-Casiano et al., 2013;González-Vega et al., 2020). Our study was focused on the copepod community as the main zooplankton group and the areas affected and unaffected by the volcano were compared. The waters surrounding El Hierro and particularly those affected by the volcano were sampled in two different seasons, spring and autumn, analyzing the main variations of the abundance, biomass, and structure of the zooplankton community. Our dataset represents the first and most extensive study of the mesozooplankton community (>200 µm) carried out for this area of the ocean after the eruption of the Tagoro submarine volcano.

Seasonal and Interannual Zooplankton Abundance
During the first years of the study (2013-2014), low zooplankton abundance was found in the surrounding area of the volcano Tagoro, as was expected after the initial devastation produced by the intense submarine eruption (Fraile-Nuez et al., 2012;Santana-Casiano et al., 2013Ariza et al., 2014). After that initial phase, a significant increase in the abundance of the mesozooplankton was observed, particularly after 2017, suggesting that probably some environmental conditions FIGURE 7 | Cluster dendrogram of the dominant copepod taxa (>1%) over the 6 year study taken in the study area of El Hierro considering their abundance as ind m −3 after square root transformation and resemblance based on Bray Curtis similarity. could affect the abundance and structure of the zooplankton community. Previous to the eruption, low productivity was observed around El Hierro Island (Arístegui et al., 1997;Hernández-León and Ikeda, 2005) as well as during the first years of the post-eruptive phase (Gómez-Letona et al., 2018). However, continuous injection of inorganic nutrients during the post-eruptive phase from Tagoro submarine volcano, particularly silicates, has been reported in the area (González-Vega et al., 2020), yielding nutrient fluxes similar to those measured in the eastern Canary Islands more influenced and closer to the NW-African coastal upwelling. Due to the physical disturbances caused by the heat released from the hydrothermal emissions, the enhanced nutrient concentrations could be recirculated together with small planktonic organisms in the surrounding waters of the volcano (Speer and Rona, 1989).
Following our cluster analysis, the zooplankton community in the initial years of the study was negatively affected with very low abundances. But the drastic increase in abundance observed after the eruptive stage indicates a recovery of the zooplankton community. The declining trend of the diversity with regards to changes in the structure of the zooplankton community and in particular for the abundance of several copepods, would show us the capacity of adaptation and resilience of this community to a different environment (Mazzocchi et al., 2012).
Overall, the highest abundance and biomass detected in spring as compared to autumn, is in agreement with a higher productivity usually found in early spring in canarian waters (Fernández de Puelles, 1986;Hernández-León, 1988;Arístegui et al., 2001;Hernández-León et al., 2004), favored by the island mass effect usually found in the leeward of the Canary Islands (Barton et al., 1998;Hernández-León et al., 2004). During the autumn, a marked thermocline characterizes the water column in this area and the high stratification limits the input of nutrients to the epipelagic waters (Arístegui et al., 1997). Conversely, in late winter, higher productivity is usually observed in the Canary Islands as a result of deeper water mixing. Accordingly, the largest zooplankton abundance recorded in March 2018 could be a result of good conditions that favored the phytoplankton growth in the euphotic zone and the subsequent increase.
In the present study, this increment in zooplankton abundance was not followed by a proportional increase in zooplankton biomass. The biomass always showed low values (<8 mg m −3 ), much lower than those reported by other studies, years before the eruption (6-14 mg m-3; Arranz, 2007) although the values were close to those reported by Ariza et al. (2014), just after the eruption. Nevertheless, it is very difficult to compare these data when different meshes were used, collecting different biomass . On the other hand, the relatively constant interannual pattern of the zooplankton biomass with a higher abundance at the end of the study would indicate the predominance of smaller size organisms.
As in other oligotrophic areas of the Canary Islands, the copepods were the dominant group of the mesozooplankton community (Fernández de Puelles, 1977, 1986Hernández-León, 1988). In the present study, the bulk of the mesozooplankton community was composed mainly of copepods which contributed to more than 70% of the total abundance. They were followed by appendicularians and chaetognaths commonly found in neighboring areas (Hernández et al., 2014;de Vera et al., 2017). In the last cruises, there was an increase in the number of appendicularians in agreement with the large aggregations reported in other marine volcanic areas of the Pacific (Lindsay et al., 2015;Pietro et al., 2019). This increase in the number of appendicularians could be related to their high ability to capture suspended picoplancton as good filter feeders since they are suited to exploit bacteria as an energy source as well as phytoplankton (Alldredge, 1981;Gorsky et al., 1999).

Diversity of Copepods and Volcanic Influence
Overall, the copepod taxa identified during the 6 years studied were close to the historical data described in the epipelagic water of the nearby island of Tenerife and El Hierro (Corral, 1970;Fernández de Puelles, 1977, 1986, where the dominant copepod genera were also Oncaea, Oithona, Clausocalanus, and Paracalanus. All the dominant species that are quantitatively important in El Hierro in more or less constant numbers are considered also cosmopolitan and ubiquitous in waters of the large subtropical open ocean (Fernández de Puelles et al., 2019). Cyclopoids as well as Clausocalanus and Paracalanus were usually characteristic of oligotrophic areas (Bucklin et al., 2010). The small calanoids (Clausocalanus, Paracalanus, Calocalanus, and Acartia) as well as the non-calanoids, Oncaea, Oithona, and Corycaeus, dominating the epipelagic waters, are generally observed at other locations around the Canary Islands (Corral, 1970;Fernández de Puelles, 1986). Although high diversity is another relevant feature of these oligotrophic areas (Corral, 1970;Fernández de Puelles, 1977), the decline observed in diversity throughout the period could be attributed to the particular conditions of the submarine volcanic areas during the posteruptive stage (Tarasov, 2006;Chan et al., 2016). A higher proportion of the cyclopoids Oncaea and Oithona found in the last years of the study could be another relevant characteristic of the Tagoro post-eruptive phase. Accordingly, hydrothermal vent fluids in the proximity of the volcanic area in islands of the Pacific have a considerable effect on the abundance, diversity, biology, and ecology of the planktonic communities both in the phytoplankton and zooplankton communities (Kosihina and Malakhov, 1991). In addition, during the post-eruptive phases, an increase of the copepod abundance with a decline of their diversity was also observed (Tarasov, 2006;Mantha et al., 2013). In those areas, the abundance of heterotrophic microplankton and zooplankton of shallow vents depended on the distribution of the organic matter and in particular on the bacterioplankton communities. In waters of the Tagoro volcano, these analyzed communities suffered drastic changes during the eruption, but when this stage ceased, the evidence of the induced changes from the volcano was no longer observed (Ferrera et al., 2015).
Copepods have been able to colonize all extreme aquatic habitats, including the deep ocean (Fernández de Puelles et al., 2019) as well as the active hydrothermal vents (Dahms and Hwang, 2013). Compared to coastal waters, copepods at venting and volcanic areas have to cope with hydrodynamic stress more than in deeper water (Nielsen et al., 2017), tending to produce a higher abundance and a lower diversity (Tarasov, 2006) as could be the case in our study. As previously mentioned, Oncaea, as well as Oithona, were in a higher density at the end of the study, peaking in spring and autumn respectively, as well as other small calanoids such as P. parvus. C. furcatus, C. arcuicornis or even A. negligens. It is noteworthy that at the beginning of the study the calanoids predominated while at the end the noncalanoids increased in their abundance. Maybe the increase of inorganic nutrients in the volcanic area favored the environment for a higher presence of several copepods, more adapted to the new water conditions. This food source appears to benefit this population of omnivorous and detritivores zooplankton feeders, such as small Clausocalanus, Paracalanus, and Acartia among the calanoids copepods as well as the cyclopoids, Oncaea and Oithona (Benedetti et al., 2016). These last copepods have been also related to the search for an additional food sources in the form of vent bacteria and microzooplankton (Burd and Thomson, 1994). Thus, we assume that these non-calanoids copepods (e.g., Oncaea, Oithona and Corycaeus) are very active ambush feeders able to thrive on the enhanced concentrations of small microzooplankton feeding on bacteria near volcanoes and the hydrothermal vents.
During the study of the post-eruptive stage of the Tagoro volcano, we did not observe any particular change in the overall environmental conditions of the epipelagic waters around El Hierro. However, the continuous injection of inorganic nutrients, particularly silicates, in the area affected by the volcano (González-Vega et al., 2020), and the differences observed between the copepod community inhabiting the volcanoaffected areas and the unaffected waters observed in this study, suggest that the Tagoro submarine volcano is indeed changing its marine environment. This fertilization process is further supported by the fact that changes occurred in the zooplankton community and its population structure, in particular detectable in spring, when the water column is well mixed down to the depths where volcanic emissions take place. This situation, however is not observed in autumn where practically the entire photic layer is isolated from the volcanic inputs by a strong thermocline. The sharper decline in diversity of the copepod community in the waters affected by the volcano compared to the surrounding waters clearly shows that the volcanic emissions might alter the composition of the zooplankton community by modulating the productivity in the area. Moreover, changes observed in the contribution of the dominant species (%) were detectable between affected and non-affected waters, suggesting that although volcanic emissions could help the increase of the zooplankton abundance, several species were favored to live under those particular environmental conditions. Further research and longer monitoring efforts in this natural laboratory located on the westernmost side of the Canary Islands will help gain new insights into the mechanisms that control the abundance and diversity of the zooplanktonic community.

CONCLUSION
The birth of a new submarine volcano Tagoro in October 2011 together with the singular location of the island of El Hierro, immersed in the most western oligotrophic oceanic waters of the Subtropical North Atlantic Ocean far away from direct anthropogenic influences, allowed us to analyze the variation of the zooplankton community for a period of 6 years and assess its response to diffuse volcanic emissions.
During the whole sampling period, from 2013 to 2018, the zooplankton community at the reference stations around the Tagoro submarine volcano (unaffected stations), showed a high similarity with more than 85% of the variation in abundance and composition attributable to seasonal differences. Main copepod assemblages were identified in two contrasting periods of the year, spring and autumn, as good indicators of the zooplankton community in an oligotrophic area of the subtropical North Atlantic Ocean.
Conversely, the zooplankton community around Tagoro submarine volcano displayed, between stations affected and nonaffected by the hydrothermal emissions, important changes in its abundance and composition mainly driven by the dominant group, the copepods. Furthermore, a low zooplankton density was seen during the 1st years of the study, followed by an apparent recovery of the zooplankton community mainly indicated by an increase in the copepod abundance, when the cyclopoids predominated. In addition, a decline in the diversity was observed particularly in the area affected by the volcanic emissions, showing that specific species were favored under their influence.
The opportunity to monitor the zooplankton community during the post-eruptive phase over 6 years (2013)(2014)(2015)(2016)(2017)(2018) represented a unique natural experiment, allowing us to determine for the first time the effects of volcanic activity in the subtropical North Atlantic Ocean. This study also sets a baseline that will be instrumental for the efficient management of El Hierro National Park for future research (UNESCO, 2000). The faunistic richness and singularity of the pelagic ecosystem of the island of El Hierro and the necessity of knowledge for its protection stress the need for further studies of this unique shallow submarine volcanic ecosystem.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.