We built up a species-specific 18S rRNA GCN database for marine protists that contained species specific GCN per cell values, along with the estimated cellular biovolume and carbon content (C-content), which were acquired from several studies published in the literature and are listed in Suppl. material 1. The absolute GCN per cell in the original research articles were achieved through different methodological approaches such as qPCR (Zhu et al. 2005), real-time PCR (Godhe et al. 2008), single-celled qPCR (Gong et al. 2013), bioinformatics pipeline (Gong, Marchetti 2019) and single cell digital PCR (Yarimizu et al. 2021). We assumed that the different 18S rRNA GCNs listed in the GCN database (Suppl. material 1) are the absolute 18S rDNA gene copy numbers per cell and that each GCN is read once during one metabarcoding sequencing event, and therefore that ‘1 read count’ is equivalent to ‘1 GCN’.

From the GCN database, we calculated the mean, median, and standard deviations GCN per cell for each of the four established groups: Bacillariophyta, Dinoflagellata, Ciliophora and flagellated cells. Median GCN cell^-1 values were used due to data limitation and presence of outliers. These values per group were used as key components for the correction factors and development of the mathematical equations.

The equation (i) estimates the corrected metabarcoding cell relative abundances (MTB_CFcell) for each defined plankton group (g) within a single sample. Metabarcoding number of reads (expressed in GCN) are divided by the corresponding plankton group CF (GCN cell^-1), and the following division by the total sum of the four groups estimates ultimately the proportional contribution to the community (%).

The equation (ii) estimates the cellular GCN:C-content ratios (CC ratio) for each plankton group (values reported in Table 1). Median cellular C-content values (pg C cell^-1) were calculated from the microscopy (MCP) dataset for each plankton group, in order to be consistent with cellular C-content of the species found in the area (BCZ). Finally, the equation (iii) estimates the corrected C-biomass relative abundances (MTB-CFbio).

g = (Bacillariophyta, Dinoflagellata, Ciliophora, Flagellated cells)

$MTB\_CFcel{l}_g[\%]=\frac{\frac{MT{B}_g\left[GCN\right]}{C{F}_g[GCN.cel{l}^{-1}]}}{{\sum }_g\left(\frac{MT{B}_g\left[GCN\right]}{C{F}_g[GCN.cel{l}^{-1}]}\right)}$ (i)

$CCrati{o}_g[GCN.pg{C}^{-1}]=\frac{C{F}_g[GCN.cel{l}^{-1}]}{Cconten{t}_g[pgC.cel{l}^{-1}]}$ (ii)

$MTB\_CFbi{o}_g[\%]=\frac{\frac{MT{B}_g\left[GCN\right]}{CCrati{o}_g[GCN.pg{C}^{-1}]}}{{\sum }_g\left(\frac{MT{B}_g\left[GCN\right]}{CCrati{o}_g[GCN.pg{C}^{-1}]}\right)}$ (iii)

The relative abundances of the uncorrected metabarcoding results and the corrected ones from equation (i) and (iii), MTB_CFcell and MTB_CFbio respectively, were directly compared to MCP. This comparison was performed for the fifteen environmental samples.

Seawater samples were collected at 3 m depth using 4 L Niskin bottles connected to a CTD sensor (Sea-bird SBE25). The monitoring in the Belgian Coastal Zone (Lapeyra Martin et al. 2022) took place from March 2018 to June 2019 (see Suppl. material 4: Table S2) aboard the RV Simon Stevin (Vlaams Instituut voor de Zee) at the St. 330 (51°26.05'N, 02°48.50'E). Both in 2018 and 2019, throughout the spring-summer months, one extra monthly cruise was undertaken and samples were collected to closely follow the evolution of the phytoplanktonic blooms that occur during this time period in Belgian waters (Gypens et al. 2007).

The DNA samples for the study of the protistan community were collected vacuum filtering 500–800 mL of water (from Niskin) through 0.22 µm polycarbonate filters (47 mm) and storing the samples immediately at -20 °C. Total DNA was extracted from filters using NucleoSpin Soil extraction Kit (Macherey-Nagel, Düren, Germany) following manufacturer’s protocol. For a maximum efficiency of the extraction from the filters, a sample lysis step was added using 10 mL cryotubes (using a high velocity bead beater for 10 min). Up to three filters were pooled and used for DNA extraction when no sufficient biomass was found on a single filter. Standard polymerase chain reactions (PCR) were performed to amplify the universal eukaryote SSU 18S rRNA gene. Primers TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′), TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′) were used to target the V4 region of the 18S rRNA gene (Stoeck et al. 2010). PCR reactions performed had a total volume of 25 μL, containing 2.5 μL of microbial DNA (5 ng µL^-1), 5 μL of both amplicon forward and reverse primers (1 µM) and 12.5 μL of high-fidelity polymerase HotStart ReadyMix (Kapa Biosystems). Plates were sealed and the following PCR-program was run in a thermal cycler: initial denaturation at 95 °C for 3 min, followed by 25 cycles of 95 °C for 30 s, annealing at 55 °C for 30 s; extension at 72 °C for 30 s final extension at 72 °C for 5 min. All PCR products (480 bp, ~383 bp + 97 bases of primers) were verified on a 1.5% agarose gel. The following library preparation of 18S ribosomal RNA gene amplicons was performed: PCR clean-up 1, index PCR, PCR clean-up 2, library quantification, normalization and pooling following the 16S Metagenomic Sequencing Library Preparation guide (Illumina 2013). Library denaturing and sample loading to the Illumina MiSeq system was performed to perform a 300 bp paired-end sequencing using V2 chemistry.

The reads were denoised and merged with DADA2, v1.16. (Callahan et al. 2016) using default cut-off parameters and annotation reads were subsequently classified with assignTaxonomy, the DADA2 implementation of the naive Bayesian classifier method described in Wang et al. (2007), against the Protist Ribosomal Reference Database PR2 v4.12.0 (Guillou et al. 2013). As this study is focusing on the protists, all reads assigned to Metazoans were excluded from the processed ASV-tables. Preparation of the ASV-tables was done in R v4.0 (R Core Team 2018) and assigned ASV taxonomy was used as proxy to classify the number of reads into described groups. The number of final sequencing reads ranged between 23,299 and 52,256 per sample. Subsampling to equal read number was done with the function rarefy from the ’vegan’ package v2.0-10 (Oksanen et al. 2013). The metabarcoding dataset used for this study is available in Suppl. material 2 and can be found in online repositories: at DDBJ/EMBL/GenBank under the accession KFLC00000000. The raw data corresponding to the raw fastq files can be found on the online repository ZENODO, with the accession number: 10.5281/zenodo.6827112.

Seawater samples from St. 330 (100 mL) were fixed with Lugol’s solution (1% final concentration), stored in the dark until microscope analysis using Utermöhl-type sedimentation (Edler and Elbrächter 2010). Cell counting and taxonomic identification down to genus/species level were performed under inverted light microscopy (LM) for all samples. If cells could not be identified down to the species or genus level using LM, higher taxonomic levels were used.

The cellular C-biomass was calculated based on biometric factors determined for each species and cell density. Mean species-specific biometric values were obtained from published data (Olenina et al. 2006; Nohe et al. 2018) and converted to biovolumes. Calculation of biomass was performed using carbon to volume relationships for diatoms, dinoflagellates and other protist plankton found in Menden-Deuer and Lessard (2000). In the case of ciliates, cell sizes were converted into cell volumes using geometric formulae from Peuto-Moreau (1991). Biovolumes were converted using the conversion factor 0.14 pg C µm³ according to Putt and Stoecker (1989). For Phaeocystis globosa, a characteristic bloom forming species of the BCZ, three different life stages were identified and distinguished: flagellate, free-living and colonial. Biovolumes and cellular carbon content were estimated following species-specific carbon content data (Rousseau et al. 1990).

For each sample, proportional cell and C-biomass abundances from the microscopy dataset were calculated (percentage of the total cells and C-biomass per defined plankton group) and directly compared to the corrected and non-corrected metabarcoding data (i.e., relative abundance of each taxonomic group, estimated as the percentage of sequencing reads assigned to a taxonomic group compared to the total reads per sample). The microscopy results dataset used for this study is available in Suppl. material 3.

All statistical analyses, data processing and plotting were performed using R version 4.0.2 (R Core Team 2018) and ggplot2 (Wickham, 2011) was used for visualization. Kruskal-Wallis and Wilcoxon tests were used to check differences in 18S rRNA GCN and C-biomass among different taxonomic groups or methods. The significance level was set at p < 0.05. Log-log plot was used to reflect the relationship over many orders of magnitude between GCN and biomass of marine protists. To check for differences in the community composition produced by the CF and assess its effects, we applied and compared generalized linear modeling based on beta distribution using the beta regression package betareg v3.1-4 (Cribari-Neto and Zeileis 2010). Two beta regression analyses were performed to assess the efficiency of the CF to improve cell relative abundances. The first one included the proportions of microscopy cell abundances and metabarcoding reads. The second model, instead of metabarcoding reads percentages, included the corrected results; after the application of the correction factor. As our dataset also assumes the “zero” values (referring to 0% of relative abundance), thus a useful transformation in practice was used for the models: (y · (n − 1) + 0.5)/n where n is the sample size (Smithson and Verkuilen 2006). The percentage of relative abundances of the eukaryotic groups were the response variable (y) and were tested against the explanatory variables (x), methods used (i.e.: MCP, MTB, MTB-CFcells, MTB-CFbio), and interaction of the methodology with each of the examined taxonomic groups.

The 18S rRNA GCN database comprised a total of 65 species’ absolute GCNs per cell, biovolumes and estimated biomasses (see Suppl. material 1, 3). Species were classified in four taxonomic groups: Bacillariophyta (Stramenopiles), Dinoflagellata (Alveolata), Ciliophora (Alveolata) and flagellated cells. GCN data collected and sorted by taxonomic group are visualized in Fig. 1 and values reported in Table 1. Ciliophora showed the highest GCN cell^-1 followed by Dinoflagellata and Bacillariophyta, and flagellated cells being the group with the lowest GCN per cell on average by several orders of magnitude. Indeed, GCN per taxum showed significant differences among all taxonomic groups (Kruskal-Wallis, p < 0.01).

Figure 1.

(A) 18S rRNA gene copy number per cell belonging to the major single celled eukaryotic plankton groups: Bacillariophyta (Stramenopiles), Dinoflagellata (Alveolata), Ciliophora (Alveolata) and flagellated cells. See Supplementary Table 1 for the median, mean and standard deviation values. (B) Log-log correlation between the number of 18S rRNA copies and cellular C-content (pgC cell^-1). The linear regressions were significant (p < 0.05) for all groups. Data were obtained through different methodological approaches: qPCR (Zhu et al. 2005), real-time PCR (Godhe et al. 2008), single-celled qPCR (Gong et al. 2013), bioinformatics pipeline (Gong, Marchetti 2019) and single cell digital PCR (Yarimizu et al. 2021).

Regarding the cellular carbon content, a significant positive correlation (R² = 0.79; p < 0.001) was found in the GCN database dataset between the gene copy numbers and the estimated C-biomass of the single cell eukaryotes. When the same analysis was performed taking into consideration the taxonomic groups established (Fig. 1B), all groups showed a significant positive relationship with Ciliophora being the least correlated one (R² = 0.47, p = 0.014) (See Suppl. material 1 for references).

The comparison of the data produced by the three datasets used in this study, sorted by taxonomic groups, is shown in Suppl. material 4: Table S1. The comparison of the number of species per taxum included in metabarcoding (n = 1884) outnumbered microscopy (n = 165) and GCN database contained the smallest number of species (n = 65). Suppl. material 4: Fig. S1 displays a comparison of cellular C-content estimates between microscopy and GCN database. No significant differences were found for Bacillariophyta cellular C-content between both datasets. In fact, Bacillariophyta median cellular C-biomass was found to be almost equal in microscopy dataset and GCN database (37.6 and 37.9 pg C.cell-1 respectively). Regarding flagellated cells, even though the median and mean C-contents per cell were higher in GCN database, there was no significant difference found (Wilcoxon, p > 0.05). Ciliophora showed the highest difference of all groups both in median and mean biomass values. Dinoflagellata C-biomass per cell showed significant differences among microscopy and GCN database (Wilcoxon, p < 0.001). In Table 1, cellular 18S rRNA GCN values are shown (displayed in Fig. 1A). The difference between the median and mean in GCN values was noted (emphasized specially in Bacillariophyta), where means were dominated by the outliers rather than the typical values. The median GCN values per taxum were used in this study as key components for the application of CF.

The cell relative abundances for the described plankton groups belonging to the 15 samples set from the BCZ are displayed in Fig. 2. Relative abundances refer to the evenness of distribution among defined groups in the community. Whereas the relative abundances in microscopy (cells) throughout all samples were dominated by the group flagellated cells, in metabarcoding proportions (reads) Dinoflagellata was the most abundant group (avg. = 51.4 ± 19.1%).

Figure 2.

Comparison between relative abundances (%) of the cell counts from microscopy (MCP; blue), reads from DNA metabarcoding (MTB; yellow) and corrected read abundance percentages (MTB_CFcell; green) per sample. Right panels: box plots summarizing the comparison percentages of cell counts (MCP), metabarcoding reads abundance (MTB) and the corrected (MTB_CFcell) out of the 15 marine water samples per plankton group: (A) Bacillariophyta (Stramenopiles), (C) Dinoflagellata, (B) Ciliophora and (D) flagellated cells. The vertical line inside the box plots represents the median, the top and the bottom hinges correspond to the interquartile range, and the whiskers show the minimum and maximum non-outlier values. Note the scale difference between plankton groups.

Bacillariophyta relative abundances from microscopy ranged between 1.4 – 61.9%, with an average (avg.) of 16.3 ± 14.7% throughout the total set (Fig. 2A). Corresponding metabarcoding proportions presented similar values with a mean 14.8 ± 7.8%, but with a maximum of 30.15% (Sample 15). Sample 3 showed the highest difference among microscopy and metabarcoding relative abundances (produced by blooming diatom Minutocellus scriptus, see Suppl. material 2). It was noted that the relative abundances between metabarcoding and microscopy results did not follow the same trend between samples. Dinoflagellata microscopy values were generally low ranging 0.08 – 2.8% (Fig. 2B), metabarcoding relative abundances were significantly higher varying between 77.2% (Sample 3) and 19.0% (Sample 9); reaching values higher than 50% in 9 out of 15 samples. The average of all samples for Dinoflagellata was the highest among all groups and attained up to 51.4%. The case of the other Alveolates group studied, Ciliophora, was similar (Fig. 2C) with averaged microscopy values less than 1% and microscopy varying between 28.4% (Sample 12) and 0.53% (Sample 5). The group of flagellated cells was formed by organisms from different taxonomic groups generally < 20 µm. Despite the limitations for the adequate characterization of small organisms by the microscopy methodology, flagellated cells percentages (Fig. 2D) reached the highest values among the four groups ranging from 35.6% (Sample 1) to 98.4% (Sample 9) (avg. = 82.7 ± 15.4%). Microscopy values were higher than metabarcoding for all samples, metabarcoding values ranging between 5.6 – 47.9% (avg. = 23.8 ± 12.9%).

When the correction factor was implemented, the estimation of relative abundances decreased considerably (avg. = 2.7 ± 2.6%) for Bacillariophyta (Fig. 2A) in all samples. The application of the CF for Dinoflagellata, which strongly dominated the community structure in the metabarcoding dataset (Fig. 2B), was able to emulate more accurately the microscopy proportions decreasing the mean to ~ 0.4% (range = 0.04 – 1.2%). The relative abundance of metabarcoding reads attributed to Ciliophora was considerably higher than percentages obtained through MCP. Corrected values’ average was found to be 0.03 ± 0.02%, getting much closer to the proportions found in MCP. When the CF was applied, the values for flagellated cells strongly approached the microscopy cell proportion values (avg. = 96.8 ± 2.8%).

Distributions of the two beta regression analyses, performed to assess the efficiency of the CF to improve cell relative abundances, were significant (p < 0.01) and detailed results are shown in Suppl. material 4: Table S3. According to the beta regression analysis containing the corrected data, GCN-CFcell resulted the best fit out of the two models with a higher R² and Phi precision parameter (R² = 0.84 to 0.96 and Phi (φ) = 14.96 to 33.23; displayed in Suppl. material 4: Fig. S2). The models significantly explained the variance in a different way for each phylum; this could be observed as well when examining the relative abundance estimates for each group in Fig. 2.

Equation (ii) was applied to calculate the cellular GCN:C-content ratios for the species present in the microscopy dataset of the BCZ and are reported in Table 1. The ratios varied from 0.94 for flagellated cells to 27.17 for Dinoflagellata. Bacillariophyta and Ciliophora displayed similar values both lower than 5 GCN:pg C.

The comparison of the biomass estimates (in relative abundances) from microscopy together with metabarcoding values and the average effect of the correction approach (MTB-CFbio) are displayed in Fig. 3. Sample per sample comparison for the entire sample set is shown in Suppl. material 4: Fig. S3. Whereas the cell relative abundances were dominated by the group flagellated cells (Fig. 3D), the C-biomass in microscopy was generally dominated by diatoms (Bacillariophyta, avg. = 46.8 ± 24%). They were followed by flagellated cells (ranged from 3 – 80.9%; avg. = 33 ± 25.2%) led by the bloom forming species Phaeocystis globosa, outweighing diatoms occasionally in specific samples (Samples 1,10,11,12). Dinoflagellata biomass outweighed the other groups in 3 out of 15 samples, Noctiluca scintillans being the main constituent of these samples (i.e.: Samples 2, 5; max. = 71.9%) with an average biomass of 19.7 ± 22.8%. Ciliate biomass estimates throughout all the samples did not attain more than 2% and was on average 0.42 ± 0.68%.

Figure 3.

Boxplots summarizing the comparison of the 15 marine water samples biomass proportions from inverted microscopy (MCP), with DNA metabarcoding number of reads (MTB) and the corrected C-estimation values from metabarcoding reads (MTB_CFbio) using median cellular C-content values per taxum from microscopy dataset. (A) Bacillariophyta, (B) Dinoflagellata, (C) Ciliophora, and (D) flagellated cells. The vertical line inside the box plots represents the median, the top and the bottom hinges correspond to the interquartile range, and the whiskers show the minimum and maximum non-outlier values.

Equation (iii) was applied using CC ratios to estimate C-biomass proportions from metabarcoding reads (MTB_CFbio), which included median cellular C-content per taxum values from the microscopy dataset. Results showed that Bacillariophyta relative biomass corrections were underestimated from 42.3 ± 23.9% to 11.1 ± 9.2% (median) when compared to microscopy (Fig. 3A). Dinoflagellata estimated abundances were enhanced (Fig. 3B) from 10.4 ± 22.8% in metabarcoding to 5.6 ± 6.0% in MTB-CFbio, however, it was found that samples characterized by monospecific blooms (i.e: Samples 2, 13 and 15) were not adequately corrected (Suppl. material 4: Fig. S3). MTB-CFbio estimations of Ciliophora were similar to metabarcoding biomass proportions (MTB = 10.0 ± 8.9%; microscopy = 8.4 ± 7.6%). Both metabarcoding and MTB-CFbio strongly overestimated Ciliophora biomass compared to microscopy (Fig. 3C). Similarly, flagellated cells overall biomass was strongly overestimated by the application of the correction factor (Fig. 3D), from 29.9 ± 25.2% to 72.3 ± 12.6%.

The two beta regression analyses were performed to assess the effect of the correction factor to improve the C-biomass relative abundances. The first one compared microscopy and MTB. The second model compared MTB-CFbio using microscopy taxa specific median C-biomass dataset values (Suppl. material 3). Distributions of the two models were significant (p < 0.01) and detailed results are shown in Suppl. material 4: Table S4. When comparing the results of the beta regression analyses performed for microscopy proportions (displayed in Suppl. material 4: Fig. S3), and MTB-CFbio, MTB-CFbio resulted in the best fit of the beta regression model with the highest R² and Phi precision parameter (R² = 0.70 and Phi (φ) = 8.8), against the first model (R² = 059 and Phi (φ) = 7.8). The models significantly explained the variance in a different way for each phylum; this could be observed as well when examining the relative abundance estimates for each group in Suppl. material 4: Fig. S4.

Cell number

Semi quantification by metabarcoding across groups has proven to be problematic in eukaryotes larger than 2 μm (Santoferrara et al. 2019) and good correlations between the proportion of reads and absolute species abundances in published literature are rare (Aylagas et al. 2018). To date, unfortunately, there is still little consensus on the combination of the methods for both identifying and measuring abundance of microbial populations.

Our results regarding the contrast of the two methods — microscopy (cells) and metabarcoding (sequencing reads) — over the set of fifteen environmental samples from the Belgian Coastal Zone (Fig. 2) displayed different types of mismatches in the four groups studied: Bacillariophyta proportions between microscopy and metabarcoding were not significantly different (Fig. 2A). Conversely, flagellated cells were underrepresented in metabarcoding against microscopical observations (Fig. 2D). In the case of Alveolates (Ciliophora and Dinoflagellata), a strong overestimation of the relative abundances was observed in our environmental samples when sequence-based methodologies were used (Fig. 2B,C), but this phenomenon is well observed in the field (Elferink et al. 2017, 2020; Bruhn et al. 2021; Santi et al. 2021).

When the CF was applied, cell relative abundances were strongly improved in both alveolate groups and flagellated cells (Fig. 2). In fact, the hypothesis that one cell with many gene copies gives many reads during sequencing is addressed by applying the correction factor. However, for Bacillariophyta the application CF worsened the estimations, underestimating them. The factors explaining both microscopy and metabarcoding similarities and the CF effect worsening in Bacillariophyta may lie beyond the differences in GCNs (Fig. 1). In fact, Bacillariophyta biovolume range is the largest one of the four taxonomic groups studied (from 10¹ to 10⁹ μm³) (Vasselon et al. 2018), 100 times greater than for dinoflagellates (Harrison et al. 2015). Given the fact that the 18S rRNA GCN in diatoms has been significantly correlated to biovolume (Godhe et al. 2008), this huge natural variability makes the development and application of a CF challenging for Bacillariophyta. Thus, it might be necessary to apply correction factors at lower taxonomic levels in groups such as diatoms with an important morphologic diversity (i.e.: Vasselon et al. 2018).

This variability was not so obvious in the other groups studied. In accordance with the GCN database results (Fig. 1), alveolates dominated the metabarcoding dataset due to higher rRNA copy numbers than other groups. This has been already shown in previous studies (Massana et al. 2015; Käse et al. 2020; Lapeyra Martin et al. 2022). Flagellated cells underrepresentation in metabarcoding dataset is best explained by the low GCN per cell (Table 1, Fig. 1). Dinoflagellata discrepancies among morphological and molecular relative abundances might be also affected by the lack of resolution of microscopical approaches to identify < 20 μm small dinoflagellates (Käse et al. 2020). Especially for ciliates, the observed discrepancy between microscopy and metabarcoding might be produced by the preservation method used.

C-Biomass

There is a realization of the fact that number of reads is not generally well suited to determine absolute abundances (Medinger et al. 2010; Mäki et al. 2017). However, a meta-analysis of 22 metabarcoding studies targeting a wide variety of organisms found a weak quantitative relationship between carbon biomass and number of sequence reads produced (Lamb et al. 2019). The GCN database results displayed in Fig. 1B showed a strong positive correlation between GCN and cellular C-content for the four defined groups. This is consistent with Godhe et al. (2008) and Zhu et al. (2005) studies, that showed a significant correlation of cell length and biovolume in some marine protist taxa (Dinoflagellata and Bacillariophyta).

As regards the entire marine plankton community, Santi et al. (2021) performed a comparison between the output of microscopy analysis (in relative biomass abundances) and DNA metabarcoding, by using the same grouping of organisms as the one used in this study. Their results revealed that alveolates (Dinoflagellata and Ciliophora) displayed differences in proportions between the two methods tested whereas Bacillariophyta results did not vary significantly. The group flagellated cells showed even higher inconsistency between the two methodologies. In contrast, our results in Fig. 3 (and Suppl. material 4: Fig. S3) showed significant differences in all studied groups but flagellated cells.

The application of the CF to estimate biomass using GCN:C-content ratios, which included GCN and median C-content per taxum from the microscopy dataset, only improved the biomass estimations in two out of the four groups studied: Dinoflagellata and Ciliophora. This might be happening due to dinoflagellates that showed the highest GCN:C-content ratio and therefore the application of the correction factor resulted in the biggest impact correcting the relative abundance estimates.

The main weaknesses behind the use of GCN:C-biomass ratios for estimating biomass relative abundances from metabarcoding reads presented in this study are certainly attributed to the limitation of the number of species representing each of the plankton groups described. Particularly in the case of Ciliophora. Since there was a single ciliate species identified on our microscopy dataset, we were aware of the misrepresentation of the C-content values for the taxum Ciliophora. However, we decided to include this group in our biomass estimation analysis (application of equations ii and iii) in order to be able to perform the generalized linear modeling based on beta distribution for the entire community and assess the general effect of the correction factor.

In addition, there are more sources of bias related to the method that should be considered. For example, the estimations of C-content proportions in microscopy might be partially caused by errors in cell size measurements (and use of mean species-specific values), which are cubed to give volumes, and the effect of fixatives, which may cause shrinking or swelling of cells (Godhe et al. 2008). Additionally, there exists a difference between small cells of diatoms and dinoflagellates, which have a ∼60% higher carbon density (carbon per unit cell volume) than large cells (Harrison et al. 2015). These reasons add up inaccuracies that go far beyond the 18S rRNA GCN bias in metabarcoding, and constrain the use of the correction factor to estimate relative biomass abundances.

DNA metabarcoding still has practical limitations (inherent to the techniques) that our GCN correction factor did not address. Some technical biases are more important than others, and amplification by PCR together with primer biases are one of the largest sources (Nichols et al. 2018). Van der Loos and Nijland (2021) suggested that PCR-free methods such as direct sequencing of rRNA genes without amplification is the obvious solution. Nonetheless, PCR dependence on metabarcoding is not likely to disappear shortly.

Two major limitations were encountered as regards the GCN database and the application of the correction factors to estimate cell and C-biomass relative abundances from metabarcoding reads. First, due to the fact that few marine protists’ 18S rRNA GCNs have been assessed (Yarimizu et al. 2021) (see Suppl. material 1 and supplementary figure W3 in de Vargas et al. (2015)), we found ourselves with the necessity to use broad taxonomic phyla or plankton groups (i.e. Bacillariophyta and flagellated cells). Along with this first limitation, it must be highlighted that high variances can be found among organisms of same phyla. Nevertheless, in order not to let the outliers control the mean values, medians were used among the available values. Intra-genomic polymorphisms were used as well, which is well documented in protists (Gong et al. 2013; Pawlowski et al. 2007). Apart from its possible negative effects for species identification, it might play a role biasing the proposed CF and reducing its effectiveness.

Secondly, the GCN values used for the application of equations (i), (ii) and (iii) (see Table 1) were strictly associated with the specimens found in our GCN database, which might not necessarily coincide with a high taxonomic level resolution with the species found in our sampling area (BCZ). Accordingly, the representativeness of the defined plankton groups based on a limited number of species for each of the plankton group is constraining, likewise the replicability of the approach. Does this mean that many more gene copy number measurements of planktonic eukaryotes will enable the creation of an ideal correction factor? Fortunately, there is a considerable improvement in the measurement of 18S rRNA GCN thanks to the emergence of bioinformatic approaches that use whole genome next generation sequencing such as used by Gong and Marchetti (2019) or Sharma et al. (2021). Indeed, this tool presents promising advances — fast and cost-effective — to overcome the organismal limitations of the 18S rRNA GCN dataset. On this basis, the presented approach may become more and more robust with the unceasing increment of marine protists’ GCN assessment. The GCN per cell information provided in Suppl. material 1 proceed from different methodological approaches, and confidence values were rarely reported in these studies. We believe that for the future GCN quantification studies on protists it should be considered adding exactitude parameters since some values might be more precise than others. Apart from this, it would also be interesting to consider comparing the different technical approaches and/or standardizing robust methods for future GCN quantification studies on protists. Nevertheless, other biological biases must also be pointed out due to the possible effect on the efficiency of the CF. For instance, GCN per cell will change according to the growth phase and physiological status of the cell. Indeed, Gonzalez-de-Salceda and Garcia-Pichel (2021) found that the number of 18S RNA genes per cell follows an allometric power law of cell volume with an exponent of 2/3. In addition, information regarding ribosome number variation would help to anticipate potential biases that can occur in environmental 18S rRNA metabarcoding dataset from marine environments.

Traditionally, trend analyses and dynamics of autotrophic plankton biomass have been often based on chlorophyll-a (Chl-a) pigment concentration by fluorimetry (Strickland, Parsons, 2002). To date, combination methods such as cell counting or quantitative PCR remain the only means to estimate absolute abundances (Weber and Pawlowski 2013; Canesi and Rynearson 2016; Vasselon et al. 2018; Santoferrara et al. 2020). The GCN-CF does not calculate the absolute cell numbers or C-mass per water volume, but provides relative abundances within the community. The standardization of a protocol that combines corrected metabarcoding relative abundances and Chl-a values might contribute to move forward towards a qualitative and quantitative monitoring of marine protistan plankton, using metabarcoding as core methodology. This way, the approach presented in this work enables us to bridge the gap between a gold standard technique in plankton research (metabarcoding) and traditional phytoplankton biomass assessment methods (Chl-a).

At the present time, community relative abundances are useful and reliable in the context of ecological interpretations (Piwosz et al. 2020). Even so, we believe that acknowledging the 18S rRNA GCN differences among taxa (the highest taxonomic resolution, the better) should become necessary in plankton metabarcoding over the next few years, since the application of a robust gene copy number CF will result in more accurate representations of the eukaryotic community structure. Even if clade specific values for corrections might be used as well, we must be aware that this presents its own problems due to high variances between GNC among groups (i.e.: ciliates).