The Artificial Substrate Units sample (hereafter “ASUs sample”) was made up of three nylon pan scorers tied together with a cable tie and attached to a stainless-steel rod that was anchored to the substrate. ASUs samples were deployed by scuba diving at a depth ranging from 14 to 25 m in three locations across the Central Mediterranean Sea: Rovinj (Croatia, North Adriatic Sea), Palinuro (Italy, South Tyrrhenian Sea) and Livorno (Italy, North Tyrrhenian Sea). Within each location, three sites were randomly selected at some kilometers of distance from each other: Rovinj: Piccola Figarola (FIG), Bagnole (BAG), San Giovanni (SAN); Palinuro: Punta Spartivento (SPA), Cala del Ribatto (RIB), Costa Azzurra (COS) and Livorno: Calafuria (CAL), Sonnino (SON), Isola (ISO) ( Figure 1 , Table S1 ). In each site, three ASUs samples were randomly deployed a few meters apart from each other, in summer 2018 and recovered in summer 2019 (Year 1), while new ASUs samples were deployed in summer 2019 and recovered in summer 2020 (Year 2). Each ASUs sample was recovered by placing it inside a plastic bag for preventing the loss of organisms, paying attention not to touch the collected sample to avoid contamination. Due to exceptional rough maritime conditions, two ASUs samples were lost in winter 2019 (one in Livorno, SON site, and one in Palinuro, SPA site) and two in winter 2020 (both in Palinuro, RIB site). Thus, a total of 50 ASUs samples were processed across the two years (25 ASUs samples per each year).

In laboratory, ASUs samples were placed in 1 L plastic jars with molecular-grade absolute ethanol, and the water content of the plastic bag was added after sieving at 40 µm sterile stainless-steel sieve. The jars were then stored at 4°C until processing. At each step, all the consumables were bleached to prevent sample cross-contamination.

The morphological identification of specimens was performed on one ASUs sample (replicate) randomly chosen per each site, for a total of nine samples per each time. The COI metabarcoding was only performed on all 2019 ASUs samples. See Figure S1 for a general outline about sample processing.

In laboratory, the three pan scourers that made up each ASUs sample were separated. Each pan scourer was shaken vigorously in a bucket containing 5 L deionized water to remove loose material and organisms. Then, each pan scourer was unraveled to pick out the remaining organisms that remained stuck in the mesh, under the stereoscope.

The bucket content was then sieved through embedded 500 µm and 40 µm sterile stainless-steel sieves. The macrofaunal fraction (≥ 500 µm) was stored apart in molecular-grade absolute ethanol in a 250 mL sterile plastic jar to perform morphological identification of the organisms. The remaining sample fraction (40-500 µm) was stored in another sterile 250 mL plastic jar filled with 95% absolute ethanol. To avoid any possible external contamination the bench and all laboratory material used were bleached between one ASUs sample and the other.

The organisms were initially separated under stereoscope into four major taxonomic groups: Annelida, Arthropoda, Echinodermata, Mollusca. All the organisms that did not fit into these groups were placed in the “Others” group. For each taxonomic group, the organisms were identified using stereoscope and light microscope to the lowest taxonomic level possible with the help of the dichotomous keys and other relevant scientific literature, and then counted. The organisms belonging to the group “Others” were only counted, but not identified (see Table S2 for the complete list of taxa found in ASUs samples over the two sampling years). After identification, all the macrofaunal organisms from each ASUs sample were merged with their respective 40-500 µm fraction for further genetic analysis (hereafter “DNA bulk sample”).

To obtain a well-homogenized product, each DNA bulk sample was transferred into an ice blender and blended at full speed for 15 seconds. The homogenate was poured into the 40 μm mesh net and washed with molecular-grade absolute ethanol several times. The resulting wet homogenate was gently squished to let the aqueous phase go through the net to dry the fraction as much as possible. This dried homogenate was then placed in a sterile 15 ml falcon tube with molecular-grade absolute ethanol, shaken vigorously and stored at -20°C.

The DNA extraction on about 1 ml of each sample was carried out using the Macherey-Nagel NucleoSpin Soil kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) (Mahmoudi et al., 2011) following the protocol provided by the manufacturer. Each DNA extract was quantified using a Qubit 2.0 fluorometer (Thermo Fisher Scientific Inc., Waltham, MA, USA) using the Invitrogen Qubit dsDNA BR Assay Kit. Samples were amplified using the newly designed degenerate primer set targeting a 313-bp fragment of cytochrome c oxidase subunit 1 (COI) (forward IIICRrev: GGNTGAACNGTNTAYCCNCC; reverse HBR2d: TAWACTTCDGGRTGNCCRAARAAYCA) that has been developed to amplify a wide range of phyla across eukaryotes and algae. Forward and reverse primers were indexed with 11–13 nucleotide long tags differentiated by at least three different nucleotides (for details see: Corse et al., 2017), and a 0 to 3 nt heterogeneity spacer (none/N/NN/NNN) was added to mitigate the issues caused by low sequence diversity in Illumina amplicon sequencing (Fadrosh et al., 2014). A unique combination of tags was used for each sample. A high sensitivity Taq polymerase was employed (Qiagen Multiplex PCR Kit, Qiagen, Hilden, Germany), and the PCR reaction was performed in 15 ml of final volume (7.5 μl of Multiplex Mastermix, 2.5 μl molecular grade water, 1.5 μl of each  10 mM primer, 2 μl of DNA) using a 2720 thermal cycler by Applied Biosystems (Thermo Fisher Scientific Inc., Waltham, MA, USA). Extraction controls consisting of negative (aerosol and  nuclease-free water samples) and positive (known DNA pools of  marine and terrestrial multi-species mock communities, data available under request to Dubut V.) controls were included in all PCR reactions. The methodology is described in Thomasdotter et al. (in press). The used amplification protocol was as it follows: (i) 15 min 95°C, (ii) 5 cycles (30 s 95°C, 40 s 45°C, 1 min 72°C), (iii) 30 cycles (30 s 94°C, 40 s 48°C, 1 min 72°C), (iv) 10 min 72°C. Three independent PCR technical replicates were made to perform post-sequencing quality checks. Every single PCR plate content was pooled in a single 1.5 ml sterile Eppendorf tube, according to relative sample quantification using a Qubit 2.0 fluorometer. Bulk pools were purified using Macherey-Nagel NucleoSpin Gel and PCR Clean-up kit (Thermo Fisher Scientific Inc., Waltham, MA, USA), according to manufacturer’s protocol. Single indexing PCRs using a couple of individual P5, and P7 Illumina adapters were performed on each purified technical replicate pool. The total indexing PCR volume was 60 μl (19.8 µl H2O, 30 µl Mastermix, 2.1 µl of each Illumina index, 6 μl cleaned PCR pool) and was performed using Thermo Scientific Phusion High-Fidelity PCR Master Mix (Thermo Fisher Scientific Inc., Waltham, MA, USA). The indexing PCR protocol consisted in (i) 30 s 98°C, (ii) 12 cycles (40 s 98°C, 45 s 55°C, 1 min 72°C), (iii) 10 min 72°C. An electrophoresis run using 1.5% agarose gel was then performed for quality checking and correct indexed DNA bands were cut out with a scalpel and purified using Macherey-Nagel NucleoSpin^® Gel and PCR Cleanup kit, following the protocol provided by the manufacturer. The three technical replicates were then quantified (Qubit 2.0 fluorometer), properly diluted and pooled at final 20 µM into a sterile 1.5 ml Eppendorf tube. A final quality check was then performed using a Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA). The libraries were then sent out to a sequencing facility (Cirad, https://www.cirad.fr/en, France) and sequenced via Illumina MiSeq benchtop system using v2 chemistry (250 bp paired-end sequencing).

Raw FastQ files from Illumina sequencing were supplied to us by the sequencing facility. Using FastQC software (Andrews, 2010), the quality of raw sequencing files was checked. OBITools package (Boyer et al., 2016) was used for trimming and joining forward and reverse reads, for removing unaligned sequences and for samples demultiplexing. The resulting 6 demultiplexed fastq files were then merged, and the amplicon length of 313 bp ±10% was selected. Chimeric sequences were sorted out using VSEARCH (Rognes et al., 2016), and non-chimeric ones were clustered (97%) into Molecular Operational taxonomic Units (MOTUs) using SWARM (Mahé et al., 2015) with d = 1 following indication by Brandt et al. (2021) that observed comparable OTU numbers for COI clustering at d = 1–13. Moreover, clustering at d-values higher than 1 can led to the loss of true species diversity, particularly for taxa known to be poorly resolved, as in the case of those living in marine natural habitats (Brandt et al., 2021). Singletons were removed through Metabarpark scripts (https://github.com/metabarpark/R_scripts_metabarpark), and the taxonomic assignment involved RDP Classifier (Lan et al., 2012) with default parameters. The COI barcoding reference database employed was COInr-Med+, which includes marine-only Mediterranean taxa, basing on OBIS (www.obis.org) and including new barcodes from central-western Mediterranean Sea (Mugnai et al., 2023). Possible contamination of samples was checked removing MOTUs present in negative and positive controls employed. All the above commands were run individually per each of the three technical replicates, and a final MOTUs table including mean values of the three runs was obtained. The final MOTUs table was normalized to the minimum library size using the online platform MicrobiomeAnalyst (https://www.microbiomeanalyst.ca, MDP protocol), with default normalization parameters. From this normalized MOTUs table, taxonomic lineages were assigned and 0.8 bootstrap threshold (Wang et al., 2007) was chosen for the taxonomic assignments of MOTUs per each taxonomic level (Claesson et al., 2009) ( Table S3 ).

Descriptive faunal analyses of the main taxonomic groups were done for morphology and metabarcoding datasets. Alpha diversity as taxa/MOTUs richness (S) and faunal abundance (N) (in terms of number of individuals and number of recovered reads per MOTUs, respectively) were calculated for each sample. Permutational multivariate analysis of variance (PERMANOVA; Anderson, 2017) based on Euclidean distance matrices of untransformed data, was used to evaluate (i) spatial-temporal variability of richness and abundance of benthic communities based on morphology (two factors: Location and Year crossed), and (ii) variability of richness and abundance of benthic communities at different spatial scales based on metabarcoding data (two factors: Location and Site nested in Location). To visualize similarity in community structure among samples for each method and sampling design, non-metric multidimensional scaling (nMDS) was performed. Then, differences in the community structure for each of the two-factor models were analyzed by PERMANOVA based on Bray-Curtis similarity of fourth root transformed data.

The similarity percentages analysis (SIMPER; Clarke, 1993) was performed on morphologic data to identify taxa that most contributed to dissimilarities in community structure detected trough morphology among Locations and Years.

To evaluate the effect of taxonomic sufficiency, the following analyses were performed on 2019 datasets. First, the data for the same 9 replicates used for the morphologic analysis were extracted from the whole metabarcoding dataset. Then, a preliminary analysis to evaluate whether the “Others” group could influence the results was done. For each method two abundance matrices with and without the group “Others” were built. The corresponding similarity matrices obtained by Bray-Curtis index, were compared by RELATE analysis, which compares pairs of similarity matrices through the Spearman rank correlation coefficient (Rho) (Somerfield and Clarke, 1995). Owing to the very high correlation between matrices with and without the group “Others” for both methods, new matrices were obtained from the initial species-per-sample matrix by aggregating raw data at increasingly higher taxonomic levels (genus, family, order, class, and phylum) without the “Others” group. Specifically, for the metabarcoding dataset, we performed the aggregation matrices starting from the species matrix obtained keeping the MOTUs by filtering for bootstrap values ≥ 0.80. For each taxonomic level for both methods, similarity matrices were obtained by Bray-Curtis similarity after the fourth root transformation. The similarity in community structure among Locations for each matrix was displayed by nMDS and formal significance tests for differences in community structure among Locations were done by one-way PERMANOVA followed by pairwise comparison when Location factor resulted significant. The resulting similarity matrices were analyzed by RELATE routine. Spearman rank correlation (Rho) was calculated between the corresponding elements of each pair of similarity matrices for species abundances and similarity matrices derived from abundances of higher taxonomic levels, separately for the morphological and metabarcoding datasets, and between morphological and metabarcoding similarity matrices at the corresponding taxonomic levels. The inter-matrix correlations matrix was used as input similarities in the creation of the ‘second stage’ nMDS ordination plot (Somerfield and Clarke, 1995) to visualize the degree of similarity among resemblance matrices obtained using various levels of taxonomic resolution for a given dataset, and between matrices obtained by the two methods at the same taxonomic resolution.

All the PERMANOVA analyses were carried out using unrestricted permutation of raw data and 9999 permutations. When significant main factors or interactions were detected, the specific procedure provided within PERMANOVA was used for a posteriori pairwise comparisons. If less than 900 unique values in the permutation distribution were available, asymptotical Monte Carlo p-values (pMC) were used instead of permutational p-values. All statistical analyses were carried out using the software PRIMER (version 7.0.21) (Clarke and Gorley, 2015, www.primer-e.com), with significant criterion for all tests set at α = 0.05.

Table comparing the species found by the morphological and metabarcoding methods in 9 samples from 2019 was created. To understand why some species were found only with one method and not with the other, for the species detected only by morphology, research of their COI sequences in NCBI (https://www.ncbi.nlm.nih.gov/nuccore) and BOLD (https://boldsystems.org/index.php) genetic databases was performed.

A total of 1,775 and 4,037 specimens of macrofaunal invertebrates were retrieved in 2019 and 2020 ASUs samples respectively, belonging to 211 taxa from the four phyla: Annelida, Arthropoda, Echinodermata, Mollusca. In general, the abundance was higher in 2020 than in 2019 ( Figure 2A ). In 2019 Mollusca was the most abundant phylum in Rovinj, while Annelida was the most abundant in Livorno and Palinuro. In 2020 mollusks were the most abundant in all three locations, followed by annelids in Livorno and Palinuro, and arthropods in Rovinj. Echinodermata had low abundance overall with higher values in 2020. Taxa richness was also higher in 2020 than in 2019 ( Figures 2B , 3 ). Annelida was the richest phylum overall, followed by Mollusca. Echinodermata had the lowest taxa richness across all locations and years, except for Livorno in 2019 where Arthropoda was the least rich group. Taxa richness was significantly different between Locations, Years and for their interaction ( Table 1 ). The results of the post hoc analysis carried out on the interaction revealed that Taxa richness in 2019 was significatively higher in Rovinj compared to Livorno and Palinuro (ROV > LIV; ROV > PAL) ( Table 1 ). Moreover, in 2020 higher taxa richness were detected in Palinuro and Rovinj in respect to 2019 ( Table 1 ). Structure of the faunal assemblages was significatively different between Locations, Years, and for their interaction. Post-hoc tests evidenced significative differences in structure of the faunal assemblages between the locations Palinuro and Rovinj in 2019, as well as between the two Years in Palinuro and Rovinj ( Table 2 ). The nMDS plot confirmed differences in the structure of the faunal assemblages between Years and between Locations within Years ( Figure 4 ). In 2019, SIMPER analysis ( Table S3 ) evidenced the highest average dissimilarity in community structure between Palinuro and Rovinj (68.63%). This difference was mostly due to the higher abundance of the polychaete annelids Filograna sp. and Vermiliopsis infundibulum (Philippi, 1844) in Palinuro, and higher abundance of the bivalve Hiatella arctica (Linnaeus, 1767) in Rovinj, but also due to the exclusive presence of several species of polychaetes, crustaceans and molluscs in Rovinj. Taxa that mostly contributed to the dissimilarity in community structure between Livorno and Palinuro in 2019 (63.79%) were polychaete annelids Filograna sp., V. infundibulum and Spirobranchus cf. polytrema (Philippi, 1844) that were more abundant or exclusive (S. cf. polytrema) in Palinuro, and echinoderm Ophiotrix fragilis (Abildgaard in O.F. Müller, 1789) and the bivalve H. arctica that were more abundant in Livorno. Dissimilarity between Livorno and Rovinj in 2019 (59.95%) was mostly due to the exclusive presence of decapod crustaceans Galathea sp. and Pilumnus spinifer H. Milne Edwards, 1834 and bivalves Striarca lactea (Linnaeus, 1758) and Hiatella striata (Fleuriau de Bellevue, 1802) in Rovinj, as well as to the higher abundance of amphipod Jassa sp. also in Rovinj. Several other polychaete species present exclusively in Rovinj contributed to this difference.

In 2020, the average dissimilarities in community structures between the three locations were lower than in 2019 ( Table S4 ). The highest average dissimilarity occurred between Palinuro and Rovinj (60.41%). Bivalves H. arctica and Limaria tuberculata (Olivi, 1792) were more abundant in Rovinj, while Arcuatula senhousia (W. H. Benson, 1842) was more abundant in Palinuro. Moreover, bivalve Nuculidae ind. and crustaceans Tanidacea ind. were exclusively present in Rovinj. Average dissimilarity in community structures between Livorno and Palinuro in 2020 (59.58%) was mostly due to the higher abundance of bivalve H. arctica and exclusive presence of bivalve Mimachlamys varia (Linnaeus, 1758) in Livorno and higher abundance of gastropod Jujubinus sp. and exclusive presence of annelid polychaete Thelepus cincinnatus (Fabricius, 1780) in Palinuro. The lowest dissimilarity in community structure in 2020 was between Livorno and Rovinj (56.10%) and was mostly due to the higher abundance or exclusive presence of several mollusks and crustaceans in Rovinj. Crustaceans Tanaidacea ind. and Isopoda ind., gastropode Jujubinus sp. and bivalve H. arctica were more abundant in Rovinj, while bivalves Nuculidae ind. were exclusively present in Rovinj. Comparing community structure between 2019 and 2020 at each location, the highest difference resulted in Palinuro (72.86%) ( Table S4 ) because of the presence of the polychaete Filograna sp. only in 2019, and the presence of crustacean Jassa sp., the gastropod Jujubinus sp. and the polychaetes Hesiondae ind. only in 2020. Moreover, gastropod Bittium reticulatum (da Costa, 1778) was more abundant in 2020. Numerous other species appeared exclusively in 2020, contributing to the difference in community structure between the two years in Palinuro Location. In Livorno, the average dissimilarity in community structure between 2019 and 2020 resulted 66.54% and was mostly due to the higher abundance or exclusive presence of several species of polychaetes, mollusks and crustaceans in 2020. The lowest average dissimilarity in community structures between 2019 and 2020 was noticed in Rovinj (62.18%) and was particularly due to the higher abundance or exclusive presence of several species of mollusks and crustaceans in 2020.

A total of 13,258 MOTUs were obtained after contamination checking and bootstrap threshold filtering procedures, dropping to 596 MOTUs after data normalization ( Table S3 ). To better describe the communities of the ASUs samples, MOTUs were then grouped into the same five categories considered for the morphological approach (Annelida, Polychaeta, Mollusca, Echinodermata, and Others; see Table S5 for reads abundance per each taxonomic level). Annelida resulted to be the most abundant phylum in all Sites and Locations, followed by Arthropoda ( Figure 5A , Table S5A ). The exception was for RIB site (Palinuro), where Arthropoda resulted the most abundant phylum, followed by Annelida. Mollusca and “Others” had similar abundances, and their differences in abundance depended on sites. The “Others” group was mainly composed of Cnidaria (60-80%) and Porifera (18-35%), followed by Chordata (1-8%) and Bryozoa (1-3%). A very low presence of Myzozoa, Nematoda, Nemertea and Platyhelminthes were detected as well, albeit not accounting for more than 1% in total. After removing the “Others” group, Arthropoda resulted the richest phylum in all sites and locations, with higher values in Livorno, followed by Annelida which showed the similar pattern ( Figure 5B ). Mollusca and Echinodermata resulted evenly distributed across all Sites, with the former accounting for 3-4 times more species than the latter.

Although Livorno accounted for 20% MOTUs more than Palinuro and Rovinj (which were comparable in terms of MOTUs richness), no significant differences in MOTUs richness were revealed between Locations and Sites ( Table 3 , Figure 6 ).

The nMDS plot showed a relatively clear separation of the communities’ structure among Locations, but also a high variability within each Location ( Figure 7 ). In fact, PERMANOVA evidenced statistically significant differences in the structure of benthic communities based on metabarcoding analyses across Locations and Sites nested in Locations. Post-hoc tests indicated significant differences between Livorno and Rovinj, and between Livorno and Palinuro ( Table 4 ).

Before analyzing the correlation between morphological and metabarcoding methods in detecting a spatial variability of the communities’ structure, a preliminary analysis was carried out. A RELATE analysis considering the 4 target phyla (Annelida, Arthropoda, Echinodermata, Mollusca) with and without the “Others” group was done both for morphological and metabarcoding matrices. The correlation between “with Others” and “without Others” matrices for the morphological method resulted highly significant with Rho value of 0.894. Similarly, the correlation between “with Others” and “without Others” matrices obtained by metabarcoding approach resulted significant as well, with Rho value of 0.763. Owing to these results, comparisons among matrices obtained at different taxonomic resolution within each method (morphological and metabarcoding) and between the two approaches at each corresponding taxonomic resolution were done without the group “Others”.

The MDS plots obtained from fourth root transformed abundances data of the different taxonomic levels for morphological data appeared similar at level of species, genus, and family ( Figure 8 ). This was confirmed by the high correlation coefficient between species and genus and genus and family ( Table 5 ). The MDS plots obtained by abundance data aggregated to higher taxonomic levels (order, class, and phylum) appeared to be less similar compared to those at lower taxonomic levels and the correlation coefficients between species and higher taxonomic levels resulted lower ( Table 5 , Figure 8 ). PERMANOVA results showed significant differences in faunal structure among Locations for all taxonomic levels (species level: Pseudo-F = 2.83, P(MC) = 0.0163; genus level: Pseudo-F = 2.542, P(MC) = 0.023; family level: Pseudo-F = 2.721, P(MC) = 0.025; order level: Pseudo-F = 3.245, P(MC) = 0.023; class level: Pseudo-F = 4.113, P(MC) = 0.022; phylum level: Pseudo-F = 17.068, P(MC) = 0.0002) and pairwise comparisons were statistically significant between Palinuro and Rovinj for each taxonomic level, while at order, class and phylum levels differences also occurred between Livorno and Rovinj ( Table S6 ). MDS ordinations of matrices derived from metabarcoding data, after fourth root transformed abundances of various taxonomic levels are reported in Figure 8 . They appeared very similar at the levels of species, genus, and family, indicating that at these levels the overall patterns of spatial variation of community structure were retained, as confirmed by high correlation coefficients ( Table 5 ). High correlation values (> 0.80) occurred only between species and genus, species and family and genus and family. The MDS plots obtained by abundance aggregated to order, class, and phylum levels, appeared to be less similar compared to those at lower taxonomic level and the correlation coefficients between species and higher taxonomic levels resulted <0.80 ( Table 5 , Figure 8 ). PERMANOVA revealed significant differences in faunal structure for the Location factor at species (Pseudo-F = 2.051, P(MC) = 0.049), genus (Pseudo-F = 2.162, P(MC) = 0.042), and phylum (Pseudo-F = 2.812, P(MC) = 0.049) levels. However, pairwise comparisons were not statistically different ( Table S6 ).

When comparing the MDS plots obtained from the morphological and metabarcoding methods at corresponding taxonomic level, the spatial pattern of the community structures resulted different ( Figure 8 ). This was confirmed by the very low correlation coefficients ( Table 5 ), below the 0.80 threshold suggested by Somerfield and Clarke, 1995.

Patterns among matrices were visualized through 2nd stage MDS plot ( Figure 9 ), where each symbol represents a similarity matrix, and the relative distance between symbols reflects the degree to which the similarity matrices are correlated. The plot showed a clear separation of the two methods on the first dimension and a separation of groups for taxonomic levels on the second dimension. Regarding morphological approach, matrices were arranged from species to phylum in a progressive disposition with species, genus, and family closer to each other compared to phylum, class, and order levels. On the other hand, the order, class, and phylum matrices from metabarcoding data resulted arranged around the cluster species, genus, and family.

The above observed patterns may be explained by comparing the species found by morphological and metabarcoding methods in the four chosen phyla (Annelida, Arthropoda, Echinodermata, Mollusca) from 9 samples ( Table 6 ). Seventy-three species were identified using morphology and 81 using metabarcoding, while only 21 species were detected by both methods ( Table 6 ). Fifty-two species were detected only through morphology but not with metabarcoding. Among them, Annelida represented 62%, followed by Mollusca (28%), Arthropoda (6%) and Echinodermata (4%). Of these species, 26 (50%) do not have COI sequences available neither in NCBI nor in BOLD databases. On the other hand, 21 species have COI sequences present in both online databases, and some of them also with a high coverage in terms of number of sequenced individuals (e.g., arthropods Jassa marmorata Holmes, 1905 and Pilumnus spinifer H. Milne Edwards, 1834; polychaetes Thelepus cincinnatus (Fabricius, 1780) and Spirobranchus triqueter (Linnaeus, 1758) which is present in the databases with its older synonym Pomatoceros triqueter (Linnaeus, 1758); ophiuroid Ophiothrix fragilis (Abildgaard in O.F. Müller, 1789); bivalve Arcuatula senhousia. Five of the species detected only by morphological analyses have COI sequences available in only one of the two databases. 60 species were detected only through metabarcoding but not with morphological analyses. Of these, 24 species (40%) are annelids, particularly representatives of the family Syllidae (13 species), 17 species (28%) are arthropods, 12 species (20%) are mollusks, and 7 species (11.7%) are echinoderms. The most unbalance between the species exclusively identified with a single approach can be attributed to Annelida phylum in metabarcoding.