18S rRNA (eukaryotes)

The 18S small-subunit rRNA gene is the eukaryotic counterpart of 16S: one conserved gene, present in every eukaryote, with variable regions (V4, V9) that carry enough signal to separate protists, microalgae, fungi and small metazoans. In MetaFlux it is the marker of choice for surveying eukaryotic diversity in soil, water, gut and plant material.

Read proportions are not cell proportions, and 18S is worse than 16S here

Every rRNA-gene survey is distorted by gene copy number, but the scale of the problem differs by orders of magnitude between the two markers. Bacterial genomes carry roughly 1–20 copies of the 16S gene, about a tenfold spread. Eukaryotic rDNA is arranged in tandem arrays whose copy number ranges from a handful to tens of thousands per genome, with ciliates and some dinoflagellates at the extreme end.

A ciliate carrying thousands of copies and a small flagellate carrying a few can be equally abundant as cells and differ by three orders of magnitude in reads. Treat an 18S table as a comparison of the same taxon between samples, not as a census of relative cell abundance within one sample. Copy-number correction for protists is not generally available, so compare composition between samples processed in the same way.

From raw FASTQ to ASVs the path is the same as for every other amplicon marker — primer removal, truncation, DADA2 denoising, merging, chimera removal. Two things are specific to 18S, and both are handled by the marker pack workflow/markers/18S.yaml:

• the length probe and the taxonomy come from different databases;
• the taxonomy reference ships two files that disagree on how many ranks they carry.

Neither needs manual intervention beyond setting the marker type and the taxon filter.

Turning it on

Set type: 18S and supply the 18S primers. A minimal amplicon block:

amplicon:
  type: 18S
  primers:
    fwd: /path/to/fwd_primer.fasta
    rev: /path/to/rev_primer.fasta
    orientation: fixed

  expected_length: auto
  probe_length_stat:
    18S: p95

  extraction:
    enabled: false          # 18S has no region extractor

  taxonomy:
    method: sintax          # rdp | sintax — both are available for 18S
    filter:
      enabled: true
      keep:    [d__Eukaryota]
      discard: []

Setting	Why it is here
type: 18S	Selects the marker pack, and with it both reference databases, the rank models, and the "no extractor" fact.
expected_length: auto	Lets the in-silico PCR probe size the amplicon window from the primers actually supplied. Fine for every V4 pair on a long enough read budget; not for V9 — see the primer table below. auto resolves to the statistic named in probe_length_stat (p95 by default), not to the typical amplicon length. For TAReuk V4 on SILVA-Euk that is ~494 bp against a 381 bp median, so the merge-overlap constraint asks for truncLen_R1 + truncLen_R2 ≥ 506 bp. A 2×300 run clears that comfortably; a 2×250 run has only ~462 bp left after primer trimming and stops at pick_trunclen. On a tight read budget, either set probe_length_stat.18S: median (or q3), or give expected_length a number of its own.
probe_length_stat.18S	Which statistic of the probed length distribution becomes the resolved expected_length — the single amplicon size used for the truncLen/overlap constraint. p95 is the shipped default: 18S V4 has a ragged upper tail, and p95 keeps it permissive without chasing the outliers that p99 would include. The ASV length-filter window is a separate calculation — in auto mode it uses the probe's q1 and p95 when a probe has run, and falls back to expected_length ±15% when none has, see below. The key is required for 18S whether or not expected_length is auto — if it is missing the workflow stops at start-up with a message naming it.
extraction.enabled: false	There is no 18S region extractor. Left true it is forced off with a warning, so this is a tidiness setting rather than a functional one.
taxonomy.filter.keep	The keep/discard lists have no built-in default. A 16S list here would discard every eukaryotic ASV, so the marker's own tokens have to be set explicitly.

Both classifiers work for 18S: rdp (DADA2 assignTaxonomy against the PR2 DADA2 file) and sintax (VSEARCH against the PR2 UTAX file, with its domain moved to the SINTAX letter d: so that VSEARCH classifies the domain first; see below). See Keeping and discarding taxa for how the two paths are made to behave identically.

Two databases, two jobs

18S is the one marker where the database that measures amplicon length is not the database that assigns names.

Database	Role in an 18S run	Why
SILVA-Eukaryotic 18S v132	In-silico PCR probe substrate only — never used for taxonomy	It recovers 96–97% of its references for every V4 primer pair tested, against 66–75% for PR2. A probe that finds nearly every reference sizes the length window from the real amplicon distribution rather than from whichever subset happened to match.
PR2 v5.1.1	Taxonomy reference (both rdp and sintax)	The better-curated eukaryote/protist reference. It is not eukaryote-only — see below.

PR2 also carries prokaryotic and organellar references, which is useful

Counted directly in the file the workflow downloads (240,201 references):

Domain field	References
Eukaryota	223,210
Bacteria	8,028
Eukaryota:plas (plastid)	6,765
Eukaryota:mito (mitochondrial)	1,890
Eukaryota:nucl (nucleomorph)	148
Archaea	131
Eukaryota:apic (apicoplast)	22
Eukaryota:chro (chromatophore)	7

Eukaryotic primers are not eukaryote-specific: they also amplify bacterial SSU and organellar (plastid, mitochondrial) rRNA, which is exactly why PR2 carries dedicated Bacteria, plastid and mitochondrial entries in the table above rather than omitting them. Because PR2 can name those products, they are assigned to those references and not to the nearest eukaryote, and keep: [d__Eukaryota] then removes them. It removes the nucleomorph, apicoplast and chromatophore entries too, since their domain field is not plain Eukaryota either; add them to keep if they belong to your study, for example keep: [d__Eukaryota, d__Eukaryota:nucl].

PR2 files its prokaryotes into the same eukaryote-style ranks, with older phylum names: a bacterial ASV reads, for example, d__Bacteria;sg__Terrabacteria;dv__Firmicutes-Firmicutes_X;… under sintax. Do not expect these names to match those of a 16S table.

Three files are fetched on first use and cached under refdb/, and a fourth is built from one of them. The three downloads are plain files, so their fetch rules are generated from the pack — there is no 18S-specific download code in the workflow. The fourth, pr2_sintax, is PR2's UTAX file with its domain moved to the SINTAX letter d: (see below); rule convert_pr2_sintax builds it when a sintax run needs it and it is missing, and again whenever pr2_utax or the conversion script is newer than it.

Pack symbol	Cached at	Source
silva_euk	refdb/silva_euk/silva_132.18s.dada2.fa.gz	Zenodo record 1447330
pr2_dada2	refdb/pr2/pr2_SSU_dada2.fasta.gz	PR2 GitHub release v5.1.1
pr2_utax	refdb/pr2/pr2_SSU_UTAX.fasta.gz	PR2 GitHub release v5.1.1
pr2_sintax	refdb/pr2/pr2_SSU_UTAX_domain-first.fasta.gz	built from pr2_utax by rule convert_pr2_sintax (method: sintax only)

Nothing about these needs declaring in a run config. An entry under references: is only needed to relocate one — for example to point at a shared copy on a cluster filesystem instead of the project-local refdb/.

PR2 carries species names inside the trainset itself, so the DADA2 addSpecies step that 16S uses does not run here (taxonomy_species_db is null in the pack). Species assignments come straight out of the main classification.

The probe result is cached per primer pair, keyed by a hash of the two primer FASTA files, at refdb/cache/probe_18S_silva_euk_18s_v132_<hash>.json. Changing the primers changes the hash and re-probes automatically; reusing the same primers across runs reuses the cached measurement. Details in Amplicon length and truncation.

No region extractor

16S runs can pass their ASVs through Metaxa2 and ITS runs through ITSx, trimming each ASV back to the marker region before the length filter sees it. MetaFlux has no equivalent step for 18S — Metaxa2 can handle SSU, but it is deliberately not wired in here — so extractor is none in the pack.

Because no rule produces an extracted FASTA, leaving extraction.enabled: true would make the workflow ask for a file nothing creates. Instead it is forced off at start-up, with a warning on stderr:

[MetaFlux] warning: marker 18S has no target-region extractor; ignoring
amplicon.extraction.enabled: true (no extraction step will run)

The consequences are small:

• The ASV length filter still runs, on the post-DADA2 sequences rather than on extracted ones. In auto mode its window is [probe_q1 − window_margin, probe_p95 + window_margin].
• Taxonomy is assigned from 05.dada2/seqs_lenfilt.fasta, the length-filtered ASVs, as it is for every marker.
• Flanking conserved sequence that an extractor would have trimmed stays on the ASV. Since the probe measures the same primer-defined fragment, probe distribution and ASV distribution remain comparable, which is what the length filter needs.

PR2 ships two rank depths

The two PR2 release files do not describe taxonomy at the same depth. Verified against v5.1.1:

PR2 file	Used by	Ranks	Rank names
pr2_SSU_dada2.fasta.gz	method: rdp	9	Domain; Supergroup; Division; Subdivision; Class; Order; Family; Genus; Species
pr2_SSU_UTAX.fasta.gz	method: sintax, through a copy with the domain first (see below)	8	Domain; Supergroup; Division+Subdivision merged into one field; Class; Order; Family; Genus; Species

Most markers use one rank model for both classifiers. 18S declares two — tax_levels for the rdp path and sintax_tax_levels for the sintax path — so each classifier is told the true depth of the file it is reading. If the depth were left at DADA2's default of seven ranks (Kingdom to Species), assignTaxonomy would still return all nine PR2 columns — it sizes its output from the reference, not from the rank names it was handed — but it names those columns from taxLevels, so only the first seven would get names, and they would be the wrong ones. PR2's Supergroup would be labelled Phylum, its Family labelled Species, and Genus and Species would never reach the taxonomy string at all. DADA2 would report no error.

Rank prefixes are keyed to what a rank means, not to the letters the underlying file happens to use. PR2's UTAX file labels its ranks k:, d:, p:, c: …, with the Domain under k: and the Supergroup under d::

>AB353770.1.1740_U;tax=k:Eukaryota,d:TSAR,p:Alveolata-Dinoflagellata,c:Dinophyceae,…

VSEARCH does not read SINTAX letters as names. It reads them as a fixed top-down order, d, k, p, c, o, f, g, s, and counts the confidence of a rank only over the bootstrap rounds that agreed on every rank above it. On the file as published, the Supergroup would be VSEARCH's top rank and the Domain would sit below it, so an ASV whose Supergroup stayed under sintax_cutoff would lose its Domain too, even when nearly every round agreed on it. MetaFlux therefore classifies against a copy of the file with the two letters swapped, the Domain under d: and the Supergroup under k:. Rule convert_pr2_sintax builds the copy, headers only, with the sequences untouched; it stops on a file that does not have PR2 5.1.1's letter order, so a future PR2 release with a different layout is noticed rather than misread.

For ASV_71 of the 18S test set (PRJDB12231), whose best hit in PR2 v5.1.1 is Angulamoeba microcystivorans (Amoebozoa), VSEARCH's confidences read:

letters as published : d:Amoebozoa(0.66),k:Eukaryota(0.66),p:Evosea-Evosea_X(0.66),…
domain first         : d:Eukaryota(1.00),k:Amoebozoa(0.66),p:Evosea-Evosea_X(0.66),…

With sintax_cutoff: 0.8 the first gives the ASV no taxonomy at all; the second gives it d__Eukaryota. Nothing below the domain changes.

The tables render the Domain as d__ and the Supergroup as sg__, to match the rdp output. Only the raw VSEARCH output, 06.taxonomy/_sintax_raw.tsv, shows the swapped letters, as in the second line above: there d: holds the domain and k: the supergroup. The taxonomy string of a dinoflagellate ASV reads, under the two methods:

rdp    : d__Eukaryota;sg__TSAR;dv__Alveolata;sbd__Dinoflagellata;c__Dinophyceae;…
sintax : d__Eukaryota;sg__TSAR;dv__Alveolata-Dinoflagellata;c__Dinophyceae;…

Only the merged Division/Subdivision field differs. Everything above and below it is identical, which is why a d__, sg__, c__ … s__ token behaves the same under either method. A Division token does not: under rdp, dv__ holds the Division alone (dv__Opisthokonta), while under sintax it holds Division and Subdivision joined by "-" (dv__Opisthokonta-Fungi), so dv__Opisthokonta matches nothing on a sintax run, and dv__Opisthokonta-Fungi nothing on an rdp run. The start-up check of the filter tokens cannot see this, because the joined form is a property of PR2's names, not of the rank, so check a dv__ token against the taxonomy column of a run made with the same method. sbd__ exists only under rdp: under sintax an sbd__ token stops the run at start-up.

Full prefix set, in order:

Rank	Prefix	Present under rdp	Present under sintax
Domain	d__	yes	yes
Supergroup	sg__	yes	yes
Division	dv__	yes	merged into dv__
Subdivision	sbd__	yes	merged into dv__
Class	c__	yes	yes
Order	o__	yes	yes
Family	f__	yes	yes
Genus	g__	yes	yes
Species	s__	yes	yes

Rank columns in the output tables

06.taxonomy/taxon_seq_table.txt has one column per rank, so an rdp run of 18S produces nine rank columns and a sintax run eight. The taxonomy string column in asv_table.txt is built from the prefixes above and is directly comparable between the two. Scripts that assume a fixed set of seven ranks will need adjusting for this marker; pick the rank columns by name — see reading the tables in R.

PR2 already stores species as complete Genus_species binomials. On the rdp path the string builder recognises this and emits them as they are, instead of gluing the genus on a second time the way it does for a bare lowercase SILVA epithet.

Setting the taxon filter

The keep/discard lists are rank-aware token lists with no hidden default: whatever is listed is exactly what is applied. For 18S the sensible starting point is

taxonomy:
  filter:
    enabled: true
    keep:    [d__Eukaryota]
    discard: []

Each token is matched against a whole ;-delimited segment of the taxonomy string, so d__Eukaryota matches the Domain field and nothing else. keep runs first and defines what survives; discard then prunes what is left.

A 16S filter on an 18S run

A 16S keep list written for MetaFlux v2.4.0 or earlier, [k__Bacteria, k__Archaea], stops the run at start-up, before any step: 18S tables have no k__ slot. The current 16S list, keep: [d__Bacteria, d__Archaea], does pass that check, because PR2 holds bacteria and archaea too, and then keeps only the few prokaryotic ASVs, removing every eukaryotic one. Should a list remove every ASV, the workflow stops with an error naming the marker and the number of ASVs removed, rather than write a header-only table. Either message almost always means the filter lists were copied from another marker.

keep also removes the unclassified

An ASV that SINTAX or RDP left unassigned at Domain has no d__ segment, so a keep: [d__Eukaryota] list drops it along with any genuine off-target sequence. Such ASVs are few (1 of 1,268 on the 18S test set, PRJDB12231, with sintax): both classifiers judge the domain on its own (for sintax this is what the swapped letters above are for), so an ASV that cannot be placed in a supergroup usually still reads d__Eukaryota. keep: [] keeps every ASV. Which is right depends on whether unassigned ASVs are noise or the interesting part of the dataset.

To drop lab or kit contaminants that a negative control revealed, add them to discard using the exact rank-prefixed name as it appears in the taxonomy column of asv_table.txt. More on both directions in Keeping and discarding taxa.

Amplicon length by primer pair

Expected amplicon length is a property of the primers, not of the marker. The figures below were measured in silico against SILVA-Euk across the common eukaryotic pairs:

Region	Primer pair	Core length	Read budget needed to merge	expected_length
V4	TAReuk454FWD1 / TAReukREV3	~380 bp	2×300	auto
V4	Parfrey 515F / 1119r	~565 bp	beyond 2×300	auto
V4	Hadziavdic 566F / 1200R	~595 bp	beyond 2×300	auto
V4	515Y / 926R (V4–V5 in 16S nomenclature)	~551 bp	beyond 2×300	auto
V9	Euk1391F / EukBr	~130 bp	any	manual, e.g. [115, 150]

"Core length" is the amplicon between the primers, after primer trimming. The spread across pairs is wide — 380 to 595 bp for V4 alone — which is exactly why the probe measures the actual pair rather than assuming a marker-wide number.

Three of these pairs cannot be merged on a standard paired-end run

The overlap constraint is truncLen_R1 + truncLen_R2 ≥ expected_length + min_overlap. A 2×300 run leaves roughly 560 bp once the primers are trimmed off both reads, so only the shortest V4 pair clears it:

Primer pair	Constraint at p95	2×250 (~460 bp)	2×300 (~560 bp)
TAReuk454FWD1 / TAReukREV3	~506 bp	fails	clears
515Y / 926R	~679 bp	fails	fails
Parfrey 515F / 1119r	~692 bp	fails	fails
Hadziavdic 566F / 1200R	~729 bp	fails	fails

The run does not produce a silently truncated result: pick_trunclen stops with the deficit in bp, so the problem surfaces before DADA2 is invoked. The options are a longer read length, a shorter primer pair, or accepting that the pairs will not merge. probe_length_stat.18S: median relaxes the constraint by sizing it from the typical amplicon rather than the p95 tail. With it, 515Y / 926R needs ~564 bp, which is marginal on 2×300.

V9 must be set manually

Only 41% of SILVA-Euk references (24% of PR2's) are recovered in silico for Euk1391F / EukBr. EukBr sits at the 3′ terminus of the gene and most references are truncated before it, so the reverse primer simply has nothing to match on. auto would then size the window from that biased minority rather than from the real amplicon distribution. Set an explicit range instead:

expected_length: [115, 150]

V9 is also the case to watch for truncation. The insert is only ~130 bp, so on 2×250 or 2×300 reads cutadapt's 3′ pass trims away the read-through and most reads come back at about that length — while on a clean run pick_trunclen has no quality drop to cut on and takes the ceiling instead, potentially a hundred bases past the end of the amplicon. filterAndTrim discards every read shorter than truncLen, which would be all of them.

trunc_len.min_read_coverage_pct should prevent that overshoot: it caps the ceiling at the length 95% of reads reach, and after the 3′ trim that length is the amplicon length itself. V9 is not covered by the 18S test set (Test datasets), so verify rather than assume: check the ceilings reported in logs/pick_trunclen.log against the post-trimming read lengths. Setting trunc_len.mode: manual with manual_r1 / manual_r2 just below the post-trimming length remains the safe choice for a marker this short. The failure mode and its diagnosis are under Troubleshooting.

Write the TAReuk primers with the IUPAC ambiguity codes published by Stoeck et al. (2010). With the codes replaced by single bases, reads that carry the other base at those positions are trimmed only because cutadapt tolerates mismatches (amplicon.cutadapt.max_error_rate, 0.2 by default).

What the run produces

Same layout as any other amplicon run. The 18S-specific parts:

File	Contents
refdb/cache/probe_18S_silva_euk_18s_v132_<hash>.json	Probed amplicon-length distribution for this primer pair, against SILVA-Euk
05.dada2/seqs_lenfilt.fasta	Length-filtered ASVs, straight from DADA2 since no extraction runs
06.taxonomy/asv_table.txt	Counts per sample plus one PR2 taxonomy string per ASV
06.taxonomy/taxon_seq_table.txt	One column per rank — nine under rdp, eight under sintax
06.taxonomy/_sintax_raw.tsv	sintax only: the raw VSEARCH output, with the domain under d: and the supergroup under k: (written sg__ in the tables)

References

Reproduced from the MetaFlux reference list:

• Quast, C., et al. (2013). The SILVA ribosomal RNA gene database project. Nucleic Acids Research.
• Guillou, L., et al. (2013). The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucleic Acids Research. (database v5.1.1 — https://github.com/pr2database/pr2database)
• Stoeck, T., et al. (2010). Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Molecular Ecology. (18S V4 primers TAReuk454FWD1 / TAReukREV3)
• Parada, A. E., Needham, D. M. & Fuhrman, J. A. (2016). Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environmental Microbiology. (18S V4 primers 515Y / 926R — the pair known as V4–V5 in 16S nomenclature)
• Parfrey, L. W., et al. (2014). Communities of microbial eukaryotes in the mammalian gut within the context of environmental eukaryotic diversity. Frontiers in Microbiology. (18S V4 primers 515F / 1119r)
• Hadziavdic, K., et al. (2014). Characterization of the 18S rRNA gene for designing universal eukaryote specific primers. PLoS ONE. (18S V4 primers 566F / 1200R)
• Amaral-Zettler, L. A., et al. (2009). A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes. PLoS ONE. (18S V9 primers Euk1391F / EukBr)

The full list, including the tools used at every step, is on Citation and references.