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 the honest framing is compositional comparison across a consistently processed set of samples.

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. Worth knowing what `auto` resolves to: the statistic named in `probe_length_stat`, `p95` by default, and not the typical amplicon. For TAReuk V4 that is ~494 bp against a 381 bp mode, 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 parse time 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). 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
Archaea	131

This is a feature rather than contamination of the reference. 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 rather than silently forced onto the nearest eukaryotic reference. keep: [d__Eukaryota] then removes them cleanly — so the filter is doing real work, not applying a formality.

Three files are fetched on first use and cached under refdb/. All three are plain downloads, so the fetch rules are generated from the pack — there is no 18S-specific download code in the workflow.

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

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 parse time, 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 and worth knowing:

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 5.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`	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 Linnaean-seven default, 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. No error, just a silently mislabelled lineage.

Rank prefixes are keyed to what a rank means, not to the letters the underlying file happens to use. The UTAX file labels its ranks k:, d:, p:, c: … but its k: field holds the Domain and its d: field the Supergroup, so MetaFlux renders them as d__ and sg__ to match the rdp output. The same taxonomy string for the same ASV therefore reads:

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 one keep/discard token behaves the same under either method.

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

6.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 seven Linnaean ranks will need adjusting for this marker — see Output files.

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 empties the table

keep: [k__Bacteria, k__Archaea] would match no eukaryotic ASV at all. Rather than write a header-only table, the workflow stops with an error naming the marker and the number of ASVs removed. Encountering that 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. On a 120-sample 18S production run this trade was measured on a test subset: the keep list removed no off-target ASVs but discarded 385 merely unclassified ones, 12.5% of all ASVs. That run therefore went ahead with keep: [] on purpose, to retain unclassified diversity. 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	~563 bp	fails	marginal
Parfrey 515F / 1119r	~577 bp	fails	fails
Hadziavdic 566F / 1200R	~607 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, which can be enough to clear a marginal case.

V9 must be set manually

Only 24–41% of references 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 in the test set, though, 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.

Degenerate primer forms are worth checking before a run. In a 120-sample V4 production run, the same TAReuk primers with their IUPAC ambiguity codes collapsed to single bases — the form some secondary sources list — matched only ~21% of forward reads, while the degenerate forms as published by Stoeck et al. matched ~94% (forward) and ~97% (reverse) across all 120 samples. A low cutadapt trim rate is the symptom to look for.

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
`5.dada2/seqs_lenfilt.fasta`	Length-filtered ASVs, straight from DADA2 since no extraction runs
`6.taxonomy/asv_table.txt`	Counts per sample plus one PR2 taxonomy string per ASV
`6.taxonomy/taxon_seq_table.txt`	One column per rank — nine under `rdp`, eight under `sintax`

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.