A total of 18 samples were collected in 2024 during routine river biomonitoring in the Ohře River basin (Czech Republic) conducted by the state enterprise Povodí Ohře. Sampling followed a standardized biofilm collection protocol CEN/TR 17245 (2018). Negative field controls were collected regularly throughout the sampling campaign and processed alongside environmental samples. Samples were fixed in 96% ethanol immediately after collection and stored at −20 °C until DNA extraction.

Physicochemical data were obtained from the routine monitoring program of Povodí Ohře following their standard operating procedures. Annual mean values of water chemistry parameters for the studied river profiles were provided by the authority. Detailed information on sampling dates, site coordinates, and physicochemical parameters is given in Suppl. material 1: table S1.

Total genomic DNA was extracted using a modified CTAB/isopropanol precipitation protocol (Urbánková and Štěpka 2025). DNA extraction was performed in two replicates per sample, with negative extraction controls included in each batch. Artificial internal spike-in DNA was added prior to extraction but was not recovered at sufficient levels for further use. Sequences are provided in the Suppl. material 1: table S3.

Two rbcL markers differing in length (263 bp and 331 bp) were targeted. Here, we use the term “marker” to refer to an amplified region defined by a specific primer pair within the coding region of the same gene. This usage is common in the metabarcoding literature, including studies of diatoms (Pérez-Burillo et al. 2022) and freshwater macroinvertebrates (Vamos et al. 2017). PCR amplification was performed in duplicate for each extraction.

The 263-bp marker was amplified using three forward primers (Diat_rbcL_708F1: AGGTGAAGTAAAAGGTTCWTACTTAAA, Diat_rbcL_708F2: AGGTGAAGTTAAAGGTTCWTAYTTAAA, Diat_rbcL_708F3: AGGTGAAACTAAAGGTTCWTACTTAAA) and two reverse primers (R3_1: CCTTCTAATTTACCWACWACTG, R3_2: CCTTCTAATTTACCWACAACAG) (Vasselon et al. 2017a). The PCR reaction mix contained 0.3 μL of each forward primer and 0.5 μL of each reverse primer (10 μM stocks).

The 331-bp marker was amplified using the forward primer rbcL-646F (ATGCGTTGGAGAGARCGTTTC) and reverse primer rbcL-998R (GATCACCTTCTAATTTACCWACAACTG) (Kelly et al. 2020), with 1 μL of each primer (10 μM) per reaction.

All primers used for the first PCR included TruSeq adapters. The reaction mixture contained 12.5 μL of 2× Combi PPP Master Mix (Top-Bio, Czech Republic), 1 μL of template DNA, primers as specified above, and nuclease-free water to a final volume of 25 μL.

The thermal cycling profile was identical for both markers and consisted of an initial denaturation at 94 °C for 2 minutes, followed by 35 cycles of denaturation at 94 °C for 20 seconds, annealing at 55 °C for 45 seconds, and elongation at 72 °C for 1 minute, with a final extension at 72 °C for 7 minutes.

PCR replicates were pooled and purified using magnetic beads (Mag-Bind TotalPure NGS, Omega Biotek, USA) with bead-to-sample ratios of 1.4 and 1.2 for the 263-bp and 331-bp markers, respectively. Library preparation and paired-end sequencing (2 × 300 bp) were performed by Eurofins Genomics (Ebersberg, Germany) using an Illumina MiSeq platform. Raw sequence data were deposited in the European Nucleotide Archive under project number PRJEB106074.

Raw paired-end FASTQ files were processed in R using DADA2 (Callahan et al. 2016). The 263-bp and 331-bp rbcL amplicons were analysed separately. Primers were removed with cutadapt (Martin 2011); reads lacking both primers were discarded.

Reads were split by sequencing run (n = 2 per amplicon) and processed independently to account for run-specific sequencing error profiles. After quality filtering (maxEE = c(2,2), truncQ = 2, maxN = 0), reads were truncated (263-bp: 260/240 bp; 331-bp: 280/250 bp), with truncation lengths selected based on visual inspection of per-base quality profiles in FastQC (Andrews 2010). ASVs were inferred with DADA2 (including error model learning with default parameters, dereplication and denoising), with error models estimated separately for each sequencing run. Paired reads were merged (maxMismatch = 1), sequence tables were combined across runs, and chimeras were removed using removeBimeraDenovo (consensus).

ASVs with < 2 reads were removed, and only sequences within expected lengths (263–266 bp; 331–334 bp) were retained. Reads were partitioned into off-target green plastids, internal standards (see Suppl. material 1: table S3 for sequences), and target red-lineage ASVs. To minimise stochastic noise from rare sequences, only red-lineage ASVs with a relative abundance ≥ 0.01% in at least one sample were retained for downstream analyses.

Off-target red plastid lineages were evaluated by BLAST against NCBI (accessed April 2026). The top 10 hits were inspected; if all were unclassified, the first hit with a confident assignment at class level was used. Off-target red-lineage ASVs were retained.

Consistency of technical replicates was evaluated after chimera removal. The relationship between replicate detection and ASV abundance was assessed using a binomial detection model based on total read counts. To account for sampling stochasticity, this relationship was additionally validated by repeated rarefaction (Suppl. material 1: fig. S1).

The effect of rarefaction to the minimum depth (89,023 reads per sample) was evaluated using Bray–Curtis dissimilarity and NMDS with Procrustes analysis (Suppl. material 1: fig. S2) and showed no meaningful effect; rarefaction was therefore not applied. All datasets were normalised to relative abundances.

For direct comparison, 331-bp ASVs were truncated to the overlapping 263-bp region (3' end) and collapsed to core sequences. Core sequences were classified as shared or amplicon-specific; prevalence and maximum relative abundance were computed and visualised (Suppl. material 1: fig. S3).

Correlation of relative abundances of shared core sequences across samples was examined using Spearman’s rank correlation.

Binding sites of the 263-bp forward primers were evaluated on 331-bp ASVs to assess potential primer bias. Mismatch number and position were computed for all primer variants (IUPAC), retaining the best match. Mismatch counts were compared between shared and amplicon-specific ASVs (Wilcoxon test), and positions were visualised (Suppl. material 1: fig. S4).

Reference sequences were obtained from Diat.Barcode v15.2 (Rimet et al. 2016). A 331-bp reference was prepared from aligned rbcL sequences using curated species names provided in Diat.Barcode. The full reference alignment was truncated to generate 331-bp and 263-bp reference datasets using AliView (Larsson 2014).

Taxonomy was assigned using the DADA2 naive Bayesian classifier (minBoot = 85). The 263-bp dataset and corresponding cores derived from 331-bp ASVs were classified jointly to ensure consistent taxonomic annotations. Taxa assigned only to genus level were retained as Genus sp.

The full reference dataset was screened for mismatches in both forward and reverse primer binding sites, with particular focus on marker-specific species. Special attention was given to mismatches occurring in the last three nucleotides of the primer binding region, which are critical for PCR annealing.

ASV tables were aggregated to species level. Abundances were corrected using CFv2 values from the trait table of the Diat.Barcode database to account for species-specific differences in rbcL copy number (Vasselon et al. 2018). Species composition was compared using Bray–Curtis (abundance) and Jaccard (presence–absence) distances. Correlation of community composition between amplicons was assessed using a Mantel test. Amplicon effects were tested by PERMANOVA with a paired design (stratified by sample), all implemented in the vegan R package (Oksanen et al. 2026).

IPS was calculated using IPSS and IPSVV values extracted from the Diat.barcode database and converted to the 20-point scale. Ecological status classes were assigned following AFNOR (2007): High (17 ≤ IPS ≤ 20), Good (13 ≤ IPS < 17), Moderate (9 ≤ IPS < 13), Poor (5 ≤ IPS < 9), Bad (1 ≤ IPS < 5). IPS values were compared between amplicons using a paired Wilcoxon test and Pearson correlation.

To disentangle the sources of IPS differences between markers, IPS was recalculated under four conditions: (i) full comparison between 263-bp and 331-bp datasets; (ii) comparison after harmonising sequence length (263 vs. 263core) to account for annotation effects; (iii) comparison after removal of marker-specific taxa to assess the contribution of taxon turnover; and (iv) comparison restricted to taxa shared within each sample to isolate the effect of relative abundance differences. Differences between conditions were summarised and compared across samples.

The R scripts and input files used for sequence processing and downstream analyses are archived on Zenodo: https://doi.org/10.5281/zenodo.20166294.

In total, 18 epilithic biofilm samples from streams within the Ohře River basin were analysed. The sites spanned a pronounced environmental gradient, particularly in nutrient concentrations (N–NO₃^–: 0.52–5.20 mg l^–1; P–PO₄: 0.02–0.18 mg l^–1), conductivity (11.8–105.7 mS m^–1), and suspended solids (2.5–53.8 mg l^–1), whereas pH varied within a relatively narrow range (7.17–8.20).

Sequencing yielded a median of 206,000 and 196,000 reads per sample for the 263-bp and 331-bp amplicons, respectively. After quality filtering and chimera removal, a median of 171,000 and 142,000 reads were retained, corresponding to median retention rates of 84% and 74%. Read retention across processing steps is summarised in Suppl. material 1: table S2.

Off-target amplification was negligible for both markers. Off-target green lineages accounted for 0.1% (119 ASVs) and 0.31% (13 ASVs) of raw reads for the 331-bp and 263-bp datasets, respectively. Off-target red lineages represented 0.14% (13 ASVs) and 0.04% (9 ASVs) of reads in the final datasets.

Technical replicates showed high consistency, with concordance increasing with ASV abundance and approaching saturation for abundant ASVs (Suppl. material 1: fig. S1). The 263-bp marker showed slightly higher consistency for low-abundance ASVs.

Rarefaction had a negligible effect on community composition, as indicated by very low Bray–Curtis dissimilarity and nearly identical ordination structure (Procrustes correlation = 0.9999; Suppl. material 1: fig. S2). The full dataset was therefore retained for downstream analyses.

A total of 603 and 595 red-lineage ASVs (including off-target sequences) were detected in the 331-bp and 263-bp datasets, respectively. Of these, 414 (~70%) were shared based on identical 263-bp core sequences.

Shared ASVs accounted for most reads (98.77% and 99.35% for the 331-bp and 263-bp datasets, respectively).

Relative abundances of shared ASVs were strongly correlated across samples (Spearman ρ = 0.65–0.87, median = 0.77; Fig. 1), and overall community composition was highly consistent between datasets (Mantel test, r = 0.90, p = 0.001).

Dataset-specific ASVs were generally rare and occurred in few samples (Suppl. material 1: fig. S3). Only a single 331-bp-specific ASV reached high abundance (~10%) in one sample.

Analysis of primer binding sites revealed no significant differences in mismatch number or positional distribution between shared and dataset-specific 331-bp ASVs (Wilcoxon test; Suppl. material 1: fig. S4), indicating no systematic primer bias explaining ASV specificity.

A total of 4,787 sequences from the Diat.Barcode v15.2 database passed filtering criteria and were used as the reference dataset for the 331-bp marker. The 263-bp reference was derived by trimming these sequences to the overlapping region.

Taxonomic assignment (minBoot = 85) yielded 535 (84.5%) and 497 (83.7%) annotated ASVs for the 331-bp and 263-bp datasets, respectively (including genus-level assignments treated as Genus sp.). Exact matches to the reference database were obtained for 199 ASVs (31.4%) in the 331-bp dataset and 195 ASVs (32.8%) in the 263-bp dataset.

The number of detected taxa was similar between datasets (185 for 331-bp and 181 for 263-bp), with 159 taxa shared. The remaining 26 and 22 taxa were initially dataset specific.

After accounting for differences in sequence length (i.e. comparison with truncated 331-bp sequences), the number of genuinely marker-specific taxa decreased to 21 (331-bp) and 17 (263-bp). These taxa were consistently of low abundance and low prevalence.

Six species detected only by the 263-bp marker showed mismatches in one or both primer binding regions of the 331-bp marker (e.g. Caloneis amphisbaena, Epithemia sorex, Hantzschia amphioxys). Similar primer mismatches were also observed in taxa such as Cocconeis spp., which were recovered by both markers at comparable relative abundances.

Species composition differed significantly between datasets (PERMANOVA, p = 0.001), although marker identity explained only a small proportion of variation (R² = 0.035). Paired samples were consistently more similar to each other than to other samples (p < 0.001; Suppl. material 1: fig. S5).

IPS values were highly consistent between markers (r = 0.94), although the 263-bp marker tended to yield higher IPS(20) scores (median difference = 0.890 IPS(20) units). This shift resulted in a change of ecological status class in three samples (Fig. 2).

Sequential decomposition of IPS differences showed that most of the difference was already captured when only relative abundance differences among shared taxa were considered (median effect = 0.568 IPS(20)). Adding sample-specific and marker-specific taxa caused only minor additional changes (median shifts = -0.010 and -0.002 IPS(20), respectively). In contrast, annotation bias related to amplicon length produced a larger additional shift (median shift = 0.334 IPS(20)), indicating that annotation-related effects represented the main secondary driver of IPS differences after abundance variation. Overall, these results confirm that marker-specific taxa contributed only minimally to IPS variation (Fig. 3).

This study compared two partially overlapping regions of the rbcL gene (263-bp and 331-bp) as candidates for routine DNA-based biomonitoring under the European Water Framework Directive. Although the dataset was relatively small, it covered a meaningful gradient of key physicochemical drivers of diatom communities. Both markers have previously been validated across broader gradients (Pérez-Burillo et al. 2020; Kahlert et al. 2021; Kelly et al. 2024), suggesting that their performance is not restricted to the conditions represented here. The dataset is therefore suitable for evaluating methodological consistency between markers.

Across analytical levels—from ASVs to taxonomic assignment and index calculation—both markers recovered a highly consistent dominant community signal. This dominant fraction, which primarily determines biomonitoring outcomes, was stable between datasets, resulting in similar species composition patterns and consistent IPS-based ecological classifications.

Agreement between datasets at the ASV level was strongly dependent on relative abundance. Nearly all highly abundant ASVs were shared, whereas discrepancies were largely confined to low-abundance ASVs occurring in few samples. This pattern is consistent with previous studies showing that rare ASVs disproportionately contribute to differences between metabarcoding datasets (Pérez-Burillo et al. 2021).

Rare ASVs, and the taxa they represent, are inherently more sensitive to methodological variation and stochastic detection. Three main sources contribute to this variability: primer bias, PCR stochasticity, and sampling effects (Leray and Knowlton 2017). While primer bias can in principle affect both abundant and rare taxa, PCR and sampling stochasticity disproportionately influence the detection of low-abundance taxa.

Primer bias is a well-documented source of variation between markers. Mismatches in primer binding sites can reduce amplification efficiency and may lead to underrepresentation or failure to detect affected taxa (Elbrecht and Leese 2015; Piñol et al. 2015). Analysis of primer binding sites revealed mismatches at the 3' end of the 331-bp primers for several taxa in the Diat.Barcode reference database, some of which were also present in our samples.

In our dataset, primer bias did not result in the loss of abundant taxa but instead manifested as reduced PCR efficiency, leading to minor shifts in relative abundance. For example, taxa of the genus Cocconeis were recovered by both markers at comparable abundances despite the presence of mismatches.

In contrast, primer mismatches are likely more critical for rare taxa and can contribute to their apparent absence. However, for low-abundance taxa, the effects of primer bias cannot be clearly separated from stochastic PCR and sampling processes. In our study, this concerned species such as Hantzschia amphioxys, Caloneis amphisbaena, and Epithemia sorex, which can become locally abundant under specific conditions and are relevant for IPS calculation (i.e. have defined IPSS and IPSVV values). Although their influence on IPS was negligible in this dataset, their potential underrepresentation should be considered, particularly under environmental conditions that favour their proliferation.

The multiplex primer design used for the shorter amplicon reduces primer bias and may enhance detection of low-abundance ASVs, thereby increasing sensitivity to the rare biosphere (Vasselon et al. 2017a). Nevertheless, the low prevalence of primer mismatches in our data indicates that most observed discrepancies are better explained by stochastic processes than by systematic primer bias.

Differences in species composition and IPS values can also arise during taxonomic assignment. These are primarily attributable to differences in marker length, limitations of reference databases, and the inherent behaviour of the Naive Bayesian classifier (Pérez-Burillo et al. 2022).

Together, these factors lead to differences in taxonomic resolution (both loss and gain) and uncertainty introduced by classification thresholds (Pérez-Burillo et al. 2022). Species identification is therefore strongly influenced by the coverage and quality of reference databases.

Reference coverage has improved substantially in recent years, with dominant taxa now being well represented, as reflected in robust IPS calculations. In particular, taxonomic work using integrative approaches (e.g. Kollár et al. 2026; Pinseel et al. 2017) can further improve both the content and structure of reference diatom databases. However, these approaches are not yet a gold standard in diatom taxonomy, and most taxonomic studies still rely primarily on morphological identification (Kollár et al. 2025).

Uncertainty associated with taxonomic assignment can be reduced by calculating ecological indices directly at the ASV level. This approach better exploits the full potential of NGS data by capturing both species-level patterns and intraspecific variation (Tapolczai et al. 2021; Kelly et al. 2024). However, fully decoupling from species is not desirable, as species-level traits underpin ecological interpretation and regulatory frameworks.

Despite differences at lower levels of organisation, IPS-based ecological assessment remained highly consistent between markers, reflecting the dominance of high-abundance taxa in index calculation.

Decomposition analysis provides a mechanistic explanation for this pattern: IPS differences were primarily driven by quantitative variation in shared taxa, with annotation bias related to marker length as a secondary contribution, while marker-specific and sample-specific taxa had negligible effects.

An important factor influencing species abundances is the correction factor (CFv2), which accounts for interspecific differences in rbcL copy number per cell and converts read counts to species abundances. The correction factor applied in our study was originally established for the 263-bp marker (Vasselon et al. 2018). However, differences in primer design (single pair vs. multiplex) as well as amplicon length may warrant marker-specific calibration. Deriving a marker-specific correction factor for the 331-bp region could therefore improve agreement between markers and with morphology-based assessments.

Overall, our results demonstrate that differences in marker length, which determine taxonomic resolution, do not necessarily translate into differences in ecological interpretation. Instead, the robustness of IPS depends primarily on the dominance structure of communities and on how strongly individual taxa contribute to index calculation (Tapolczai et al. 2019).

From an applied perspective, the longer 331-bp amplicon appears more suitable for routine biomonitoring. Its higher taxonomic resolution improves discrimination among abundant taxa, which are the main drivers of index values. This advantage can be further enhanced by appropriate calibration to better translate read counts into species abundance estimates.

For biomonitoring purposes, rare taxa are generally not critical for index calculation (Vasselon et al. 2017a). However, for studies focusing on the rare biosphere, it may be beneficial to adapt primer design (e.g. using multiplex strategies like those applied for the 263-bp marker) to improve detection of low-abundance taxa, while still retaining the advantages of the longer marker.

Taken together, our results indicate that the 331-bp marker provides a robust and reliable basis for routine biomonitoring under the Water Framework Directive, with performance comparable to the shorter marker and sufficient consistency for ecological status assessment. While further optimisation, such as marker-specific calibration, may improve quantitative accuracy, our findings suggest that the 331-bp marker can be applied in practice without major methodological adjustments.