A total of 561,941 raw sequences were obtained from MiSeq paired-end sequencing, and all the raw sequences were generated as contig sequences ( Supplementary Table 1 ). The sequences filtered and trimmed via Mothur were 179,827 unique sequences, and 465,118 reads were obtained. Among the sequences, 164,789 unique sequences and a total of 433,983 reads of the 18S rDNA V8–V9 region were obtained. After removal of chimeras and ‘pre.cluster’, 35,897 unique sequences and a total of 332,895 reads were passed through filtering. After singleton sequences were removed, 10,224 unique sequences and a total of 307,222 reads were obtained ( Supplementary Table 1 ).

Eight supergroups of Eukaryota were identified from metabarcoding sequences using BLAST search with PR² database ( Supplementary Table 2 and Supplementary Figure 2 ). The total read number identified by the BLAST search in the 19 sediment samples was 234,173 reads, with the highest number of 139,129 reads for Alveolata in the TSAR ( Supplementary Table 2 ). Total reads of the supergroups and divisions derived from the 19 sediment samples were varied among 5,124 reads at URA18 to 17,703 reads at URA02 ( Supplementary Table 3 ). Supergroup and division showing the most abundant reads was TSAR and Alveolata, respectively, in URA01–URA14 and URA19 samples (ranged from 33.79% at URA19 to 79.83% at URA09, Supplementary Table 4 and Supplementary Figure 2 ). In contrast, Obazoa and Opisthokonta were the most dominant supergroup and division in URA15–URA18 (ranged from 33.86% at URA15 to 42.59% at URA17), respectively.

Ten HAB species belonging to the genera Alexandrium, Azadinium, Chattonella, Fibrocapsa, Heterocapsa, and Heterosigma were identified in this study ( Table 1 ). In the case of the genus Skeletonema, the total numbers of unique sequences and reads of the genus were also shown in Table 1 because not only Skeletonema costatum but also many species of the genus are considered to cause nori bleaching in Japan. The HAB species Aureococcus anophagefferens, Azadinium dexteroporum and A. spinosum, Dictyocha fibula, Karenia mikimotoi, Pseudo-nitzschia delicatissima, P. multiseries, P. pseudodelicatissima, and P. pungens were also detected. However, the total read numbers of those HAB species were less than fifty which was difficult to examine for vertical distribution in the sediment samples, so those species were excluded from the subsequent analysis.

In the molecular phylogenetic tree of the genus Alexandrium, unique sequences of the genus Alexandrium obtained this study were separated into five major clades ( Supplementary Figure 3 ). Four clades contained sequences of A. affine, A. leei, A. pacificum (group IV) and A. tamiyavanichii, whereas the remaining clade contained two species (A. hiranoi and A. pseudogonyaulax) in the molecular phylogenetic tree. A. pacificum (Group IV) and A. tamiyavanichii have been reported as paralytic shellfish toxin (PST) producers (A. pacificum: Brown et al., 2021; Li et al., 2011; Wang et al., 2005, A. tamiyavanichii: Hashimoto et al., 2002; Lim et al., 2006; Sagara et al., 2010, Table 1 ). A. affine and A. leei have been reported to be associated with fish mortality (A. affine: Nguyen-Ngoc, 2004, A. leei: Shikata et al., 2020; Tang et al., 2007, Table 1 ). A. hiranoi and A. pseudogonyaulax produce goniodomin A, which is an antifungal polyether macrolide that inhibits actin reorganization and angiogenesis (Murakami et al., 1988, 1998; Triki et al., 2016, Table 1 ).

Regarding the toxic species of the genus Azadinium, three unique sequences belonging to the clade A. poporum/A. dalianense were obtained ( Supplementary Figure 4 ). Since A. dalianense has not been found and A. poporum has recently been found in Japanese coastal waters via many water samples collected from various locations in Japan (Takahashi et al., 2021b), these unique sequences may be those of A. ‘poporum’, known as an azaspiracid producer (Krock et al., 2015; Luo et al., 2017, 2018; Tillmann et al., 2011, 2018, Table 1 ). The azaspiracid is known to cause symptoms similar to diarrheal shellfish poisoning, including nausea, vomiting, abdominal pain, and severe diarrhea symptoms in humans who ingest mussels or scallops contaminated with azaspiracids (Twiner et al., 2008).

The unique sequences belonging to the Chattonella marina complex (C. marina var. antiqua, C. marina var. marina, C. marina var. ovata and C. minima) were detected ( Supplementary Figure 5 ). Several strains of the C. marina complex produce reactive oxygen species (ROS), which may affect the gills of fish during red tide outbreaks and cause fish mortality (Shikata et al., 2021, Table 1 ). Regarding phylogeny, the C. marina complex detected in the sediment sample may also have the potential for ROS production because sequences of the productive ROS strains and those of the C. marina complex detected in this study appeared in the same clade on the tree ( Supplementary Figure 5 ). The unique sequences assigned to Fibrocapsa japonica were detected ( Supplementary Figure 5 ). Fish mortality caused by F. japonica has been reported (Iwasaki, 1971; Landsberg, 2002, Table 1 ), and some studies have reported that this species produces ROS and hemolysins, which may cause cell clogging in the gills (de Boer et al., 2009, 2012). The unique sequence of Heterosigma ‘akashiwo’, which has been reported as a HAB that causes fish mortality, was detected ( Table 1 and Supplementary Figure 5 ). The genus Heterosigma contains a new species, H. minor, described in 2016 (Engesmo et al., 2016). These two species of Heterosigma are indistinguishable in the 18S rDNA V8–V9 region analyzed in this study. However, H. minor has only been described from a strain isolated from Virginia, USA (strain ARC HA0504-1), and no Japanese strain has been isolated to date (Engesmo et al., 2016; Imai et al., 2021; Mehdizadeh Allaf, 2023). Therefore, we determined the unique sequences in the phylogenetic tree as H. ‘akashiwo’ ( Supplementary Figure 5 ).

Three unique sequences were found in the genus Heterocapsa, but the position of those sequences in the phylogenetic tree of the genus Heterocapsa was difficult to identify at the species level, since those sequences belonged to one clade along with several other species of the genus Heterocapsa ( Supplementary Figure 6 ). This is because the 18S rDNA V8–V9 region, which is the target region for metabarcoding in this study, cannot identify each species of the genus Heterocapsa. However, considering that there are reports of bivalve mortality caused by H. circularisquama in Uranouchi Inlet (Horiguchi, 1995; Imai et al., 2006; Matsuyama, 2012; Sakamoto et al., 2021; Shiraishi et al., 2008), these three sequences were considered as H. circularisquama and named as H. ‘circularisquama’ ( Table 1 ).

The blooms of Skeletonema have been reported to be responsible for the color bleaching of nori (Pyropia spp.) in Japan (Imai et al., 2006, Table 1 ) and gill lesions and mortality of Atlantic salmon Salmo salar in British Columbia (Kent et al., 1995), unique sequences assigned to this genus were treated as HAB species in this study ( Supplementary Figure 7 ).

The HAB monitoring data provided by the Kochi Prefectural Fisheries Experiment Station for fishery damage in Uranouchi Inlet between 1984 and 2017 showed that the largest fishery damage on record occurred in 2001 in this inlet, with the damage amounting to 60 million yen ( Table 2 ).

Water quality survey results from 1980 to 2020 showed that DIN-N, DON-N, DOP-P, T-N and T-P stayed high in 2005, while inorganic nitrogen (NH₄-N, NO₂-N, NO₃-N and DIN-N) and inorganic phosphorus (PO₄-P) were intermittently high from 1980 to 2020 ( Table 3 ).

The vertical distribution of HAB species analyzed via a heatmap revealed that the eleven HAB species could be divided into three groups ( Figure 1 ). The first group was composed of six species found in samples from almost all the sediment layers: A. hiranoi/pseudogonyaulax, C. marina complex, F. japonica, H. ‘circularisquama’, H. ‘akashiwo’ and Skeletonema spp.

The second group of A. affine, A. pacificum (Group IV), and A. tamiyavanichii was not detected or was detected at low abundance in the upper layers of the core sample ( Figure 1 , URA01–08), whereas these species were detected in the deeper layers (under URA09 and URA10, whose years of occurrence were estimated to be 1980s).

The third group was composed of two HAB species, A. leei and A. ‘poporum’, which were not detected in the deep layers, unlike the URA10 and URA11 samples, respectively, but were detected in the upper layers ( Figure 1 , URA01–URA11).

Among the detected HAB species, Alexandrium hiranoi is known to form resting cysts in the winter and can remain dormant until the environmental conditions become suitable for germination (Kita et al., 1985). A. leei (Anderson et al., 2012; Shikata et al., 2020), A. pseudogonyaulax (Anderson et al., 2012; Lassus et al., 2016; Montresor, 1995; Montresor and Marino, 1996), A. affine, A. pacificum (Group IV) and A. tamiyavanichii are also known to form cysts (Anderson et al., 2012; Nguyen-Ngoc, 2004). The following other HAB species are also known to form resting cysts in the sediment: A. poporum (Gu et al., 2013a; Luo et al., 2018), C. marina complex (Imai and Yamaguchi, 2012; Ishikawa et al., 2022; Katano et al., 2014), F. japonica (Cucchiari et al., 2010; Yoshimatsu, 1987), H. circularisquama (Shiraishi et al., 2008), H. akashiwo (Mehdizadeh Allaf, 2023; Mehdizadeh Allaf and Trick, 2023) and Skeletonema spp (Ellegaard and Ribeiro, 2018; McQuoid and Hobson, 1996; Stenow et al., 2020). Therefore, those DNAs detected in the sediment samples have the possibility of originating from cysts or from DNA adsorbed on humic acids (Lewis et al., 2023; Pedersen et al., 2015; Saeki et al., 2011) in the sediment samples.

From the vertical distribution of eighteen HAB species analyzed by heatmap, the first six HAB species group (A. hiranoi/pseudogonyaulax, C. marina complex, F. japonica, H. ‘circularisquama’, H. ‘akashiwo’ and Skeletonema spp.) which appeared in almost every sediment sample suggest that these six cyst-forming species occurred continuously throughout the ages, and DNA derived from the cysts or cells were deposited in the bottom sediment resulting in their detection in all sediment samples. Since 1984, when records began to be kept by the Kochi Prefectural Fisheries Experimental Station, four genera, Chattonella, Fibrocapsa, Heterocapsa, and Heterosigma, have repeatedly formed red tides and caused economic damage to aquaculture (~60 million JPY, Table 2 ). Although there are no records prior to 1984, these records after 1984 seem to correspond to the existence of those four genera in the sediments during that period ( Figure 1 , URA01–URA09).

According to Kochi Prefectural Fisheries Experimental Station, four Alexandrium species (A. affine, A. pacificum (Group IV), A. tamiyavanichii and A. leei) detected by metabarcoding have not been reported by direct cell counting, but three Raphidophyta (genera Chattonella, Fibrocapsa, and Heterosigma) have been reported ( Table 2 ). The reason for this difference is the number of copies of 18S rDNA per cell in Alexandrium spp. and Raphidophyta. It is known that the copy number of 18S rDNA per cell in the genus Alexandrium is higher than that in Raphidophyta (Yarimizu et al., 2021). Therefore, it is likely that metabarcoding would have detected Raphidophyta but Alexandrium was not. This is a weakness of the occurrence analysis by metabarcoding and should be analyzed by quantitative metabarcoding in the future (Ushio et al., 2018; Sato et al., 2021; Tsuji et al., 2022).

The reason why A. ‘poporum’ was detected by metabarcoding but its occurrence has not been recorded may be because the size of this species is small and difficult to identify under microscopic observation (Takahashi et al., 2021b; Tillmann et al., 2011). Considering these issues, analysis of HAB species in bottom sediments by metabarcoding is useful for clarifying the occurrence history of HAB species.

In relation with the second HAB species group (A. affine, A. pacificum (Group IV), and A. tamiyavanichii) in the heatmap, Alexandrium spp. blooms occur on a large scale in Osaka Bay and Hiroshima Bay in the Seto Inland Sea, Japan, when nutrient concentrations such as DIN and PO₄ are less than 12.8 μM and 0.4 μM, respectively (Itakura et al., 2002; Yamamoto, 2019). Considering these issues, it has been reported that conditions suitable for Alexandrium spp. proliferation may involve oligotrophic or mesotrophic waters rather than eutrophic waters (Itakura et al., 2002; Natsuike et al., 2018a, 2018b; Yamamoto, 2019). Therefore, the bloom formation of the Alexandrium species in Uranouchi Inlet may have been suppressed by eutrophication caused by aquaculture after 1977–1988 ( Table 3 ), as discussed later.

There are two possibilities for why the third HAB species group (A. leei and A. ‘poporum’) in the heatmap were not detected in deeper samples than the URA11 and URA12 samples, whose years estimated by radiometric dating (Takahashi et al., 2021a) were 1977 and 1954. First, these two species did not occur in Uranouchi Inlet prior to 1977 or 1954 and appeared after those years; second, the cysts that formed prior to 1977 or 1954 might have decomposed since the cysts of these two HAB species are not highly durable. Regarding the first possibility, yellowtail aquaculture started in Uranouchi Inlet in 1959, and the aquaculture industry may have caused environmental changes such as eutrophication, as discussed later in the enclosed inlet, which may have enhanced bloom formation in these two species. The cyst durability of these two species should be investigated to clarify which possibility is conceivable.

Possible causes of the change in HAB species composition after 1977–1988 (URA09–11) include the following two possibilities: first, eutrophication of Uranouchi Inlet due to the start of aquaculture, and second, climate change, such as global warming. Regarding the first possibility, the total N (T-N) and total P (T-P) concentrations in Uranouchi Inlet in 1985 were 0.260 mg/L and 0.070 mg/L, respectively ( Table 3 ). These values suggest that the seawaters in Uranouchi Inlet were eutrophicated or polluted at that time, based on the criterion of eutrophication (T-N: 0.220– 0.650 mg/L, T-P: 0.03–0.09 mg/L, Tavakoly Sany et al., 2014) and pollution in coastal waters (‘low-level’ T-N pollution: 0.252–0.308 mg/L; and ‘very high level’ T-P pollution: > 0.031 mg/L, Smith, 2003). It is possible that the eutrophication in Uranouchi Inlet may have occurred as a result of the proliferation of aquaculture in Uranouchi Inlet from the 1960s to the 1970s, which resulted in changes in HAB species composition from 1977–1988. However, information on nutrient concentrations in the years prior to the 1980s is not available. This hypothesis is supported by the finding that Skeletonema spp., a known indicator of eutrophication in the Seto Inland Sea of Japan (Itakura and Yamaguchi, 2007; Nishikawa et al., 2010; Yamada et al., 1980a, 1980b, 2011), have been detected in greater numbers in Uranouchi Inlet since 1977–1988 ( Figure 1 ). On the other hand, there is a gap of more than a decade between the start of aquaculture in Uranouchi Inlet (1959) and the transition timing of the HAB species groups occurred (1977– 1988). This gap may be due to eutrophication caused by nutrients increasing beyond the environmental capacity while the aquaculture was gradually expanded since the start of aquaculture was not sudden on a large scale.

Second, the sea surface temperature (SST) around Japan has increased due to global warming, and the annual average SST in the northwestern Pacific around Japan has shown an increasing trend (1.24 °C/100 years, Japan Meteorological Agency, 2023). The annual average SST in Uranouchi Inlet increased by 0.19 °C/decade from the 1970s to the 2010s (Tanaka et al., 2012). Such SST trends in Uranouchi Inlet and nearby sea areas may increase blooms of HAB species such as A. leei, which have been reported to occur in tropical and subtropical areas, such as Malaysia (Usup et al., 2002), Singapore (Tang et al., 2007), Thailand (Kodama et al., 1982), Vietnam (Nguyen-Ngoc, 2004), and China (Gu et al., 2013b). This hypothesis was supported by the report that A. leei recently formed blooms and caused fish mortality in Nomi Bay in 2017 (Shikata et al., 2020), which is close to Uranouchi Inlet. Although it is difficult to directly link such changes in the bloom formation of HAB species with global warming, it may be important to consider global warming as at least one of the factors that caused the changes in the occurrence trends of HAB species in Uranouchi Inlet observed in this study.

In this study, we revealed the presence of eleven HAB species, with notable shifts from approximately 1977–1988 in Uranouchi Inlet, Kochi, Japan. This shift corresponded to two hypotheses: the development of aquaculture and the resulting eutrophication, or sea surface temperature rising due to global warming. Moreover, metabarcoding using vertical core sediment samples provides footprints of how HAB species composition has changed and maybe be affected by anthropogenic environmental changes.