Harmful algae diversity from a coastal upwelling system detected by high-throughput sequencing

Introduction In recent years, environmental DNA (eDNA) amplicon sequencing has been used to unveil plankton diversity in the field. Nevertheless, molecular methods, such as this, are rarely used in harmful algal bloom (HAB) monitoring programs, which mainly rely on morphological identification by conventional light microscopy. Methods The present study focused on a shallow marine environment (Ría de Vigo, Northwest Spain), where sediment and plankton samples were collected from 2016 to 2018. Results The application of eDNA amplicon sequencing allowed us to simultaneously detect 25 potential harmful species (mainly diatoms and dinoflagellates) included in the IOC-UNESCO Taxonomic Reference List of Harmful Microalgae. Among these, causative agents of amnesic shellfish poisoning (Pseudo-nitzschia spp.) paralytic shellfish poisoning (Gymnodinium catenatum and Alexandrium minutum), azaspiracid producers (Azadinium poporum) and ichthyotoxic haptophytes (Chrysochromulina leadbeateri), were identified. Some toxic microalgae were better represented in sediment (e.g., Pseudo-nitzschia pungens, Gymnodinium catenatum) or planktonic fractions (e.g., Pseudo-nitzschia, Gymnodinium smaydae), confirming the importance of including both sediment and plankton fractions in eDNA monitoring studies. Despite the limitations of sequencing short amplicons, it was possible to discern in this study six Pseudo-nitzschia species and associate each of them with each seasonal peak produced in summer periods. Furthermore, several species previously unreported in Ría de Vigo (Pseudo-nitzschia turgidula, Chrysochromulina leadbeateri, Azadinium poporum) could be detected. Discusion These results point out the application of eDNA amplicon sequencing to expand our knowledge about harmful species in HAB monitoring programs and early warning systems for low abundant and rare taxa.


Introduction
Photosynthetic organisms (e.g., microalgae and cyanobacteria) that inhabit benthic and planktonic marine ecosystems are essential constituents of food webs. Under certain environmental conditions, due to natural and/or anthropogenic-induced changes, their uncontrolled growth turns in some cases into harmful algal blooms (HABs), which can be deleterious for the aquatic ecosystems and cause severe economic impacts in the aquaculture industry (Berdalet et al., 2015;Grattan et al., 2016;Gobler, 2020;Wells et al., 2020;Hallegraeff et al., 2021). In recent decades, several factors, such as the expansion of human activities in coastal ecosystems, the implementation of monitoring systems, and the consequences of environmental changes (eutrophication, marine transport, global warming), could explain the apparent worldwide increase in the number, frequency, and intensity of HABs (Hallegraeff, 2004;Anderson et al., 2012;Grattan et al., 2016;Hallegraeff et al., 2021).
Almost 300 phytoplankton species have been listed as harmful (Lundholm et al., 2009, onward), and approximately 80 are known to produce toxins (Hallegraeff, 2004;Berdalet et al., 2015). Often, HABs are associated with toxin-producing species that can lead to the accumulation of these compounds in aquatic organisms, typically filter-feeding molluscs. Nonetheless, marine biotoxins can affect aquatic organisms from invertebrates to fish, mammals and seabirds, and their accumulation in the food web represents a potential risk for public health through the consumption of contaminated shellfish and fishery products (Glibert et al., 2014;Bresnan et al., 2021). Otherwise, HABs can lead to oxygen depletion in the water column during the decline of their populations, or cause physical damage (mechanical injury/irritation) to fish gills (Anderson et al., 2012).
Episodes of HABs are among the most critical environmental factors affecting the Galician shellfish industry due to the closure of harvesting installations (Rodrıǵuez et al., 2010), risking the loss of 118 M€ (production value in 2016; Avdelas et al., 2021). In Galicia, the development and effects of toxic algal blooms remained unnoticed until 1976, when a paralytic shellfish poisoning (PSP) episode developed in the Rıás Baixas, hindering mussel aquaculture. Approximately 170 people were affected in Spain and several European countries due to the consumption of contaminated mussels (Trainer et al., 2010). In 1981, over 5000 consumers in Spain were affected by diarrhetic shellfish poisoning (DSP) associated with blooms of Dinophysis sp. (Reguera et al., 2008;Trainer et al., 2010). Monitoring for the presence of toxins in shellfish by bioassay was officially adopted, and the scarce phytoplankton surveys that were initiated in the 1950s increased in frequency and focused on potentially toxic species. The relevance of the aquaculture sector pushed efforts toward the development of an intensive monitoring system for biotoxins, phytoplankton and oceanographic conditions. Moreover, the European Reference Laboratory for Marine Biotoxins (EURLMB) has also been based in the Vigo region since 1993.
Most national HAB monitoring systems carry out regular samplings at fixed stations and cell counts of potentially toxic genera/species under an optical microscope using standard procedures (Anderson, 2009). However, this technique is timeconsuming and requires highly trained staff with taxonomical skills to analyze and report the results (Karlson et al., 2010). In recent years, the advent of new technologies, such as autonomous imaging flow cytometry (IFC) (Dashkova et al., 2017), robotic platforms (drones, autonomous surface vehicles) (Wu et al., 2019), satellites (Palenzuela et al., 2020;Detoni et al., 2023) and molecular methods (e.g., qPCR, in situ hybridization) (Ebenezer et al., 2012;Zamor et al., 2012;Medlin and Orozco, 2017), have provided multidisciplinary approaches to collect information about harmful phytoplankton species and population dynamics, which are valuable for monitoring systems.
Although molecular technologies are widely and routinely applied in research, their transfer to laboratories concerned with monitoring requires a change in paradigm and the choice of suitable methods (i.e., allowing multiple detection of target organisms, providing quantitative results useful for early warning and risk assessment, etc.; FAO et al., 2023). Among the several techniques currently available, highthroughput sequencing constitutes an interesting and valuable approach. DNA metabarcoding, based on massive sequencing of eDNA, gives higher taxonomic coverage, allowing the identification of thousands of species in a single sample (Deiner et al., 2017). This methodology has been successfully used for unveiling unknown plankton diversity (Stoeck et al., 2010;de Vargas et al., 2015), detecting novel genera/species (Dzhembekova et al., 2017;Sze et al., 2018), and determining the composition and dynamics of cooccurring species in microalgal blooms (Buchan et al., 2014;Needham and Fuhrman, 2016;Zhou et al., 2018;Hou et al., 2020).
In this study, we report the results of metabarcoding targeting the 18S V9 ribosomal RNA (rRNA) gene region, with a special focus on harmful algae. Data collected in the Rıá de Vigo (Northwest Spain), in a production area of bivalve mollusks during a three-year period were used to detect toxic HAB species in surface water and surface sediment samples. This area is highly influenced by the wind-driven upwelling and downwelling events, inducing seasonal changes in the diversity of organisms (Rosoń et al., 2002;Arıśtegui et al., 2009). The results allowed to gather new information in the area about dominant and rare phytoplankton including toxic bloom forming taxa. This approach, if conducted regularly within the monitoring system, would allow to track seasonal changes in the planktonic communities, providing an early warning of potentially harmful species overlooked by conventional microscopical techniques.
Meira, a bivalve production area in the Rıá de Vigo, Northwest Spain (42°17´6.72´´N/8°43´18.80´´W) ( Figure 1A). Samples were collected every three months from July 2016 to June 2018. Surface sediment was collected from three different areas of the intertidal zone (distance of 5 m among sampling points) and then mixed. Seawater was collected at approximately 10 meters from the seashore, and successively filtered through several pore sizes to obtain different planktonic fractions according to Sieburth et al. (1978). For this purpose, 75 m 3 of seawater was filtered through a 200-µm plankton net to obtain the mesoplankton fraction. Subsequently, a volume of 40 L was filtered again through a 65µm pore size polycarbonate filter to obtain the microplankton fraction. Finally, 2 L of seawater were filtered using a 0.22-µm pore size polycarbonate filter to obtain the nanoplanktonpicoplankton fraction ( Figure 1B). All samples were kept at -20°C until DNA extraction. A total of 28 samples of sediment and water fractions were processed.

Metabarcoding assay and data analysis
The metabarcoding assay of eukaryotic communities was performed as detailed in Rıós-Castro et al. (2021). This methodology included DNA isolation, amplicon sequencing of the V9 region of the 18S rRNA gene, sequence analysis using bioinformatic tools (Microbial Genomics Module, CLC) to obtain the eukaryotic OTU (Operational Taxonomic Unit)  (A) Location of the sampling area in Rıá of Vigo (Meira, NW Spain). (B) Summary of the general methodology applied in this study. Superficial sediment was collected from the intertidal and seawater was collected at approximately 10 meters from the seashore, and it was successively filtered by different pore sizes to obtain the mesoplankton, microplankton and nanoplankton fraction. The DNA metabarcoding procedure was performed according to Rıós-Castro et al., 2021 to obtain the eukaryote OTU table. This study was focused in the microalgae community, particularly in the study of harmful microalgae.
Dinophyceae, Dictyochophyceae, Haptophyceae, Pelagophyceae and Raphidophyceae were used for further analysis since harmful algal species in the IOC-UNESCO list belong to these taxonomic groups.
To decrease the number of nonrepresentative taxa, OTUs with a combined abundance in all samples ≤ 5 were excluded from the analysis. The total number of reads and OTUs from each class was represented as well as the number of reads obtained for each class in each sediment and plankton sample. The community composition of microalgae was represented as relative abundances according to the following criteria. Overall, all genera representing >1% of the total reads were represented. When the sum of the relative abundance of all OTUs within a genus was between 0.5% and 1%, these OTUs were classified at the class level. Finally, genera with abundances lower than 0.5% were grouped as "other" because of their low representativeness in the microalgae composition. The top 20 most abundant OTUs from the classes Bacillariophyceae and Dinophyceae were also represented in sediment and plankton compartments. Microalgal OTUs considered by the IOC-UNESCO as harmful algae were selected for further analyses. To avoid the inclusion of OTUs with ambiguous identity, representative sequences of all OTUs associated with the toxic HAB species were also verified with BLASTn, and only OTUs with 100% similarity with the associated single species were considered.

Environmental parameters
Abiotic parameters are represented in Supplementary Table 1. Seawater temperature followed the same patterns as air temperature, except for July 2016, when the air temperature was 6°C higher than the water temperature. The low water temperature recorded in July 2016 coincided with the highest upwelling index recorded in comparison with other summer periods. The highest precipitation values were recorded in February 2018 and June 2018. The lowest fluorescence values were observed in winter, especially in February 2018, which coincided with the lowest dissolved organic carbon (DOC) value. The lowest dissolved organic nitrogen (DON) was obtained in summer, and maximum values occurred in winter, mainly in February 2018. The highest level of DOC was observed in September 2017.

Sequencing results and microalgae communities
The eukaryotic sequencing results were detailed in Rıós-Castro et al. (2021). After the eukaryotic OTU table subsampling procedure was performed in the present study, 6,819 eukaryotic OTUs were obtained. After filtering the main microalgal classes that may include harmful species, 444 microalgal OTUs were retained, representing 18.73% of the total eukaryotes in the sediment and 26.28% of the total eukaryotes in the plankton (Supplementary Figure 1). In terms of OTU number, the most diverse microalgae were Dinophyceae (259 OTUs) and Bacillariophyceae (146 OTUs), followed by minor contributions from Haptophyceae (24 OTUs), Dictyochophyceae (7 OTUs), Raphidophyceae (4 OTUs) and Pelagophyceae (4 OTUs) (Figure 2A). Dinophyceae and Bacillariophyceae represented 6% (66,348 reads) and 5.6% (61,654 reads) of the total eukaryotic reads ( Figure 2B) after removing those OTUs with a total abundance lower than 5 reads.
In sediment, the class Bacillariophyceae was the most represented, especially in the samples from January 2017 (11,656 total eukaryotic reads) and June 2017 (10,881 total eukaryotic reads) ( Figure 3A). Benthic diatoms of the genera Navicula (27% of total microalgal reads) and Amphora (12% of total microalgal reads) dominated the microalgae composition in sediment, with some OTUs of Navicula phyllepta, Navicula arenaria and Amphora laevis ranked in the top 20 most abundant OTUs of their class (Supplementary Figure 2). Other minority reads from Bacillariophyceae belonged to the genera Craticula, Berkeleya, Cymbella, Sellaphora and Pseudo-nitzschia ( Figure 3B). The genera Amphidinium (3.14%), and Proterythropsis (2.22%), both from Dinophyceae, were also detected in sediment ( Figure 3B). In general, seasonal variability in the microalgal community was not clearly observed among the different sediment samplings ( Figure 3B). However, in the October 2016 sampling, OTUs from the class Dinophyceae represented almost 37% of the total microalgae, with the genus Amphidinium reaching 16.37% of total reads ( Figure 3B).
Conversely, Dinophyceae was the most representative class in the planktonic fractions, with 24,137 total eukaryotic reads in September 2017 ( Figure 3C). The most representative genera were Oxyrrhis, Amphidinium, Alexandrium and Parvodinium, which accounted for 29%, 5.96%, 2.69% and 1.5% of the total microalgal reads, respectively. Other classes like the Haptophyceae and Dictyochophyceae, with Chrysochromulina sp. and Pteridomonas sp., were also well represented in the plankton ( Figure 3D). The microalgal community composition was highly variable, and the massive proliferation of certain organisms was occasionally detected, such as Oxyrrhis marina, which reached 84% of total microalgal reads, among the top 20 OTUs with the highest abundance in the Dinophyceae ( Figure 3D; Supplementary Figure 2).

Harmful algal communities
A total of 26 species listed as harmful algae according to the IOC-UNESCO list were detected in the marine environment (Figure 4). Within the Bacillariophyceae, the main genus detected was Pseudo-nitzschia, whereas Gymnodinium and Azadinium were the dominant genera in the class Dinophyceae. Minor harmful taxa such as Amphidinium carterae, Dinophysis sp., Alexandrium minutum and Karenia mikimotoi were also detected. The descriptions of harmful microalgae detected within the classes Bacillariophyceae and Dinophyceae are detailed below.

Bacillariophyceae
The most abundant potentially harmful algae detected was the genus Pseudo-nitzschia (Bacillariophyceae). It appeared in all samplings both in sediment and planktonic fractions. Despite the highest contribution of Bacillariophyceae in sediment, Pseudonitzschia was especially dominant in water samples. The highest number of reads obtained for this genus corresponded with summer months and coincided with the peaks in temperature and upwelling Relative abundance of each microalgae class (A) and genera (B) found in sediment. Relative abundance of each microalgae class (C) and genera (D) found in plankton. Relative abundance was calculated respect the total of microalgae selected in this study. B, Bacillariophyceae; D, Dinophyceae.
index. Moreover, their total number of reads decreased from summer to winter samplings during the studied period. Therefore, a seasonal pattern of Pseudo-nitzschia sp. could be detected ( Figure 5A). Several species were identified and confirmed with 100% similarity in BLASTn: P. pungens, P. fraudulenta, P. delicatissima, P. australis, P. heimii and Pseudonitzschia sp.* (100% BLASTn similarity with P. turgidula, P. cuspidata and P. lineola). All Pseudo-nitzschia spp. were associated with planktonic samples, with the exception of P. pungens, which was mostly recorded in sediment samples ( Figure 5B). Furthermore, each summer peak in Pseudo-nitzschia could be associated with a dominant species among other secondary ones ( Figure 5B). For example, P. australis and P. delicatissima dominated in the peak of July 2016, while P. pungens was dominant in summer 2017, and Pseudo-nitzschia sp.* (P. turgidula, P. cuspidata or P. lineola) in summer 2018.

Dinophyceae
Potentially harmful microalgae belonging to the genera Amphidinium, Gymnodinium and Azadinium were detected in both sediment and planktonic samples ( Figure 6). Amphidinium sp. was present in sediment samples without a clear seasonal abundance pattern associated with environmental factors (temperature and upwelling index). However, a massive proliferation of this genus occurred in October 2016 in comparison with other samplings ( Figure 6A). The species Amphidinium cupulatisquama, A. massartii, A. steinii, A. carterae and A. operculatum were detected, with A. steinii being the main species responsible for the abundance peak in October 2016 in sediment. Although almost all species of Amphidinium prevailed in sediment, A. carterae was better represented in plankton in June 2018. Furthermore, A. operculatum was exclusively detected in plankton in September 2017. Other OTUs similar to Total reads of potential harmful microalgae species found according the Taxonomic reference lis of harmful microalgae (IOC-UESCO).
Amphidinium sp. but not confirmed at 100% similarity after BLASTn analysis were also detected (Supplementary Table 2).
Gymnodinium sp. was detected in the sediment and water column in all sampling periods without a clear seasonal abundance pattern ( Figure 6C). The highest number of reads was obtained in June 2018, mainly associated with the species Gymnodinium smaydae in the water column, followed by the abundance of G. catenatum in sediment in July 2016 ( Figure 6D). Although five Gymnodinium species were detected, G. catenatum and G. smaydae prevailed, and they were associated with sediment and plankton fractions, respectively. Gymnodinium dorsalisulcum was mainly represented in planktonic samples from September 2017 and February 2018, and minority species such as G. impudicum and G. aureolum were also found in planktonic samples.
Azadinium sp. was mainly recorded in sediment, with an abundance peak associated with A. poporum. Azadinium trinitatum was also detected, but with low representativeness and exclusively in sediment ( Figures 6E, F).

Discussion
The impact of HABs in coastal ecosystems includes marine organisms, human activities (aquaculture, fisheries, and tourism) and public health issues (Grattan et al., 2016;Gallardo-Rodrıǵuez et al., 2019). HAB records have globally increased in recent decades, mostly associated with environmental changes and the expansion of human activities and parallel pressure in these areas (Davidson et al., 2014;Gobler, 2020;Hallegraeff et al., 2021). Thus, several strategies for their mitigation, prevention and control in coastal waters have been assessed thus far (e.g., Anderson, 2009;Anderson et al., 2012).
Despite this, the advent of new technologies in recent decades for the molecular identification of taxa (or genes associated with toxin production) have rarely been incorporated into routine monitoring programs (e.g., McGirr et al., 2022). Among the several techniques proposed for qualitative or semiquantitative detection of toxin-producing species, eDNA sequencing appears to be useful for the description of marine diversity, including harmful microalgae (Hii et al., 2021;Jacobs-Palmer et al., 2021) of interest in surveillance and monitoring programs.
To explore the possibilities of this approach, the Galician rias are particularly appropriate due to the intense aquaculture and fishery activities, exposed to the recurrent development of HABs and highly influenced by upwelling and downwelling events (Álvarez-Salgado et al., 2008;Dıáz et al., 2014;Broulloń et al., 2020). Marine biotoxins produced during these events cause periodical harvesting closures during the year. Thus, improving our knowledge of oceanographic conditions and toxic phytoplankton species can provide critical information for monitoring systems, stakeholders, and management purposes.
In the present study, the application of eDNA amplicon sequencing allowed us to simultaneously detect 25 potential harmful species from several algal groups included in the IOC-UNESCO Taxonomic Reference List of Harmful Microalgae (Lundholm et al., 2009, onward) in both sediment and water samples and in different seasons over a two-year period.
The most abundant harmful species (in terms of a higher number of reads) belonged to the genus Pseudo-nitzschia (Bacillariophyceae). Diatoms from this genus are well documented in the Galician rias (Fraga et al., 1998), including species associated with the production of domoic acid (DA), the causative agent of amnesic shellfish poisoning (ASP); (Trainer et al., 2010)). The accumulation of this toxin in seafood causes severe economic impacts in the aquaculture sector and human health (Reguera et al., 2008). Despite in Galician rıás P. australis is known as the main causative agent of DA production, the results from DNA metabarcoding allowed us to identify multiple Pseudonitzschia species related to DA production in other areas (Dıáz et al., 2014;Bates et al., 2018;Broulloń et al., 2020). The species P. australis, P. delicatissima, P. fraudulenta and P. pungens have been detected before in the Galician rias (Fraga et al., 1998;Trainer et al., 2010;Zapata et al., 2011;Dıáz et al., 2014;Palenzuela et al., 2020).
We also detected a sequence that could be associated with P. turgidula, P. lineola or P. cuspidata, as we obtained the similarity percentage, coverage, and e-values with all these species. The first two species have not been detected previously in Spain (Bates et al., 2018); however, P. cuspidata was previously detected on the Atlantic coasts of Spain and Portugal (Fraga et al., 1998;Churro et al., 2009). Although maritime transport could influence the geographical expansion of these species, additional studies should be performed to confirm the presence of new species in the Rıá de Vigo. Total reads of harmful Dinophyceae in total sediment and plankton and in each sediment and plankton sampling in: (A, B) Amphidinium sp., (C, D) Gymnodinium sp. and (E, F) Azadinium sp. S, sediment; P, plankton. -Castro et al. 10.3389/fmars.2023.1200135 Frontiers in Marine Science frontiersin.org

Ríos
The analysis of both sediment and planktonic fractions allowed us to reveal differences in species distribution within the same genera between environmental compartments. For example, P. pungens was more associated with sediment samples, whereas the other Pseudo-nitzschia species recorded were mainly present in water samples. This could be explained by the vertical flux and flocculation or presence of resistant forms in sediment, as previously described in Pseudo-nitzschia (Dortch et al., 1997;Parsons and Dortch, 2002;Trainer et al., 2010).
One of the main benefits of HTS-based methods is the possibility of easily monitoring changes in a community across different seasons over a period (Gilbert et al., 2012;Hernańdez-Ruiz et al., 2018;Rıós-Castro et al., 2021). Rıá de Vigo is highly influenced by the wind-driven upwelling and downwelling events, inducing seasonal changes in the diversity of organisms (Rosoń et al., 2002;Arıśtegui et al., 2009). In upwelling events, nutrient fertilization associated with the northerly winds and the increment of solar irradiation results in a high primary production, which is jeopardized due to the massive proliferation of toxin producing HABs such as, Dinophysis acuminata, Pseudonitzschia spp. and Alexandrium minutum (Figueiras and Rıós, 1993;Broulloń et al., 2020). In contrast, Gymnodinium catenatum blooms occurs in downwelling events (Broulloń et al., 2020). In this study, a seasonal pattern highly related with upwelling events was only clearly observed in Pseudo-nitzschia spp. This genus was always detected in sediment and plankton, with abundance peaks associated with summer periods from 2016 to 2018. These results agree with previous associations of its blooms with upwelling periods in late springsummer (Figueiras and Rıós, 1993;Trainer et al., 2010;Trainer et al., 2012;Broulloń et al., 2020). Furthermore, each peak in abundance of Pseudo-nitzschia spp. detected in summer months could be attributed to different species: P. delicatissima in summer 2016, P. pungens in summer 2017, and P. turgidula/P. cuspidata/P. lineola in summer 2018. In the Rıá de Vigo, two Pseudo-nitzschia bloom events, from June-July 2016 and from August 2018-October 2018, were included in the Harmful Algae Events Database (IOC-UNESCO, HAEDAT, http:// haedat.iode.org), corresponding with our sampling periods. However, they could not be identified at the species level using traditional light microscopy. It is known that Pseudo-nitzschia species can be distinguished by genetic means, but they are often morphologically identical (cryptic species), or only a few morphological differences can be found after a detailed examination (pseudocryptic species) (Trainer et al., 2012). The specific results obtained with the application of molecular methods in this study indicated that P. delicatissima and P. australis could be the main species involved in the bloom recorded in summer 2016, whereas three possible species (P. turgidula/P. cuspidata/P. lineola) could be involved in the peak in summer 2018. Although we detected an abundance peak of P. pungens in summer 2017, this possible HAB event was not recorded in the HAEDAT database. As we mentioned above, this genus was mainly detected in sediments, and given that phytoplankton monitoring is only performed in the water column, it is likely that this episode was not detected in water samples during that period.
The eDNA metabarcoding allowed the detection of harmful dinoflagellates well known to occur in the Rıá de Vigo, such as G. catenatum and A. minutum, the main toxin producers associated with paralytic shellfish poisoning (PSP) in the Northwest Iberian Peninsula. In these cases, information on their temporal trends was obtained during the study. Gymnodinium catenatum blooms usually occurred in winter periods (Fraga et al., 1993;Figueroa et al., 2008), and they were present in HAEDAT records in February-March 2017 and September 2017. However, it was not possible to associate those records with our data. One interesting finding was the high abundance of G. catenatum obtained in sediment samples in comparison with planktonic fractions, which could suggest the presence of resting cysts (Figueroa et al., 2008;Liu et al., 2020). Alexandrium minutum was detected but in general with very low abundance, despite a proliferation that caused an exceptional red tide recorded in June-July 2018 in the Rıá de Vigo (Ben-Gigirey et al., 2020;Nogueira et al., 2022). In addition, this study provided the first records of azaspiracid (AZA) producers, such as the Amphidomataceae genus Azadinium (A. poporum), in the Northwest Iberian Peninsula. These toxins have been detected in minor amounts by INTECMAR weekly monitoring analyses and should be surveyed in the near future. These emergent toxins constitute a potential risk given the harvesting closures they cause in shellfish aquaculture areas in other European countries such as Ireland. Furthermore, these small dinoflagellates cannot be identified by optical microscopy to the species level. In fact, in 2018, the Marine Institute (Galway, Ireland) validated and accredited for the first time the method for A. spinosum detection by qPCR in a routine monitoring program (Wietkamp et al., 2020).
In the context of the present study, where some potentially harmful taxa were recorded for the first time, it worths to be mentioned the absence of some benthic dinoflagellates from warmer waters, such as the genera Gambierdiscus and Ostreopsis, the latter in apparent expansion in neighboring areas of the Iberian Peninsula (Drouet et al., 2021), as well as the low number of reads (21 reads) of Dinophysis acuminata, the main producer of diarrhetic shellfish toxins (DSTs) in this area (Velasco-Senovilla et al., 2023).
Finally, one of the haptophytes detected in this study with a representative number of reads, Chrysocromulina leadbeateri, deserves also particular mention as it has been responsible for large fish kills in Norway due to its ichthyotoxicity (John et al., 2022). Even if such fish mortalities have seldom been reported in the area, the presence of potentially ichthyotoxic haptophytes and raphidophytes, such as those reported in this work, would represent another target organism readily detected by HTS analyses.

Conclusion
In conclusion, the use of metabarcoding methodology based on the sequencing of the V9 region of the 18S rRNA gene has provided valuable and new information on the composition of harmful microalgae which can directly affect marine fauna and human health. Despite the recent interest in including sediment in metabarcoding studies (Forster et al., 2016;Holman et al., 2019;Rıós-Castro et al., 2021), few studies targeted their respective communities in shallow-water habitats. In fact, this approach revealed the presence of distinct harmful microalgae representative in sediment (e.g., P. pungens, G. catenatum) and/or planktonic fractions (e.g., P. delicatissima, G. smaydae), highlighting the importance of including sediment in eDNA monitoring studies. Furthermore, despite the limited taxonomical resolution inherent to short amplicons, it was feasible to characterize species of particular genera, such as Pseudo-nitzschia, and even species previously unreported in the area (P. turgidula). Moreover, each summer abundance peak could be associated with particular Pseudo-nitzschia species. This suggests the importance of including eDNA metabarcoding as an integrative and systematic monitoring tool for harmful algae and other organisms of interest. This molecular approach can complement morphology-based analysis, which usually lacks resolution at species level due to subtle morphological differences and/or the size of some harmful organisms, previously overlooked by conventional methods which rely on time-consuming analysis by highly trained staff with taxonomic expertise.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material. The dataset generated for this study can be found in the Sequence Read Archive (SRA) http://www.ncbi.nlm.nih.gov/sra. Under the accession number PRJNA761019.

Author contributions
RR-C: Methodology, formal analysis, investigation, data curation writing-original draft, and visualization. BN: Conceptualization, resources, methodology, writing-review & editing, supervision, and funding acquisition. JH-U: Writing-review & editing, supervision, FR: Writing-review & editing, supervision, AF: Conceptualization, resources, methodology, formal analysis, writing-review & editing, supervision, and funding acquisition. All authors contributed to the article and approved the submitted version.