Two algae that were in the reads all along

A 2016 comparison of 18S primers, rebuilt on AlkhemYst from its seven mock communities.

---

The question

Most of what is known about microscopic algae in lakes, coasts and treatment plants now comes from sequencing one gene, the 18S rRNA gene, in short stretches. Which stretch to read is not a detail. A pair of primers that binds one group of organisms well and another badly will report a community that never existed, and nobody looking only at the result can tell.

In 2016 Ian Bradley, Ameet Pinto and Jeremy Guest, at the University of Illinois at Urbana-Champaign and Northeastern University, designed and tested primers for two stretches on the Illumina MiSeq [1]: V4, the region most 18S studies read, and V8-V9, near the end of the gene.

The mock communities

What makes their study worth rebuilding is how they tested the primers. They took twelve algae, six from freshwater and six marine, and cloned each one's 18S gene into a plasmid, so that every alga contributed the same number of gene copies and no difference in copy number could blur the answer. From those they mixed seven communities, from freshwater only (MC1) to marine only (MC7), and sequenced each one three times with both pairs of primers, beside samples of freshwater, wastewater and coastal seawater. Then they deposited every read.

Their findings: V4 failed to capture two of the twelve, the haptophytes Prymnesium parvum and Isochrysis galbana. V8-V9 recovered all twelve. Both regions separated the three kinds of water.

What we did

AlkhemYst's components for amplicon sequencing had been tested on 16S, the bacterial gene, and this was their first 18S study. Two canvases of eleven steps each, one per region, the same steps in both: a read quality report, primer trimming, quality filtering, denoising into exact sequence variants with DADA2 [2], taxonomy against the PR2 reference [3] with the classifier PR2 publishes for IDTAXA [4], a phylogenetic tree, the mock community check, composition, alpha and beta diversity, and a BIOM file for other tools. Each canvas became 1,028 lines of R, run on our servers in the EU with the network switched off.

The mock community check compares every sequence variant in every mock with the expected sequences of the twelve algae, which the authors deposited in GenBank, and says, sample by sample, which algae were found exactly and which were not.

The paper's result, again

V4 found ten of the twelve. No sequence variant in any of the eighteen mock samples that hold P. parvum and I. galbana matched either of them; in three, a variant one base away from I. galbana held 0.2 percent of the reads. Everything else V4 found was exact: 0.3 percent of the mock reads matched no expected sequence.

V8-V9 found all twelve. P. parvum was in sixteen of the eighteen mocks that hold it and I. galbana in fifteen; the ones that missed them were the replicates of MC2, the community with the most freshwater in it besides the first, where the marine algae are a small share.

Both regions separated the three kinds of water. PERMANOVA [5] puts the habitat behind 17 to 33 percent of the variation in V4, depending on the distance, and 20 to 38 percent in V8-V9, p = 0.001 in every case. The groups also differ in how spread out they are, so part of that separation is spread rather than position. V8-V9 found more sequence variants than V4 in 22 of the 37 samples sequenced both ways, a small difference: medians of 375 against 344 in freshwater, 231 against 214 in coastal seawater, 46 against 42 in wastewater. The paper reported more taxa with V8-V9; the direction holds, and the size is modest.

Where the two algae went

We expected the cause to be a primer that does not bind them. It is not that, as far as the sequences can say: both V4 primers match the reference sequences of both haptophytes base for base.

The run itself gave the first clue. In V4, the quality filter kept 69 percent of the reads of the freshwater-only mocks and 53 percent of the marine-only ones, the mocks where the two haptophytes are. In V8-V9 the same filter kept 83 and 87 percent. Something in the marine mocks' V4 reads was failing the filter, and it was not the primers.

The haptophytes' V4 region opens, about fifty bases in, with an unusually GC-rich run: CGGGGCGGGGCCGCCGG in P. parvum's reference. Stretches like that are hard for the sequencer. The forward read crosses it within its first hundred bases; the reverse read, which starts from the other end of a 384-base region, stops before it gets there. If that was the cause, the reverse reads alone should hold the two algae.

So we read V4 again, on a third canvas of five steps, taking only the reverse reads as single reads. All twelve algae came back. P. parvum was in all eighteen of the mocks that hold it, I. galbana in seventeen, each 17 to 23 percent of the reads in the three marine-heavy communities. The filter kept 90 to 98 percent of the reads in every mock, and 0.43 percent of the mock reads matched no expected sequence.

The paper had seen this from the other side. It reports that the haptophytes' reads carried an error rate an order of magnitude higher than the rest of the community, 2.82 percent against 0.04. The errors were real; the algae were there.

There is one more thing in the two figures worth a reader's attention. With the two haptophytes lost, the marine algae that remain look larger than they are: Heterocapsa niei is 42 to 44 percent of the marine-only mocks in V4's read pairs and 24 to 27 percent in its reverse reads. A missing member does not leave a gap; it inflates everything else.

What this means for anyone reading V4

V4 is still the region most 18S studies sequence, and it is a good one. But a V4 dataset with no haptophytes in it is not yet evidence that there were none. Reading the reverse reads alone took one canvas of five steps and ten minutes here, and it answers the question.

The part that was actually hard

The first 18S run stopped at the taxonomy step. The classifier our 16S runs use holds a profile for every name at the reference's lowest rank, and PR2 has 53,594 species: about 14 gigabytes, more than the server has. The step now uses IDTAXA with the classifier PR2 itself trains, which needs 2.

The mock community check was written for one mock, as 16S studies usually carry. This study has twenty-one, of seven different make-ups. It now checks any number, each against the algae its community should hold. And one of the twelve reference sequences, Thalassiosira pseudonana's, stops before the end of the V8-V9 region; the check now counts a variant as exact over the part a reference covers, and says so.

And uploading the data showed that an upload was held whole in memory before it was stored, which a 657 MB archive is too large for. Uploads now stream.

None of these is interesting on its own. They are here because an account in which everything worked the first time would be a worse guide to what bringing a new kind of data to a platform actually costs.

Why this story

There is nothing to correct. V4, as it is usually processed, loses these two algae, as the paper says, and V8-V9 does not. What the rebuild adds is where they go: into the forward reads, not out of the sample.

That is only knowable because of how the study was built. Seven communities of known make-up, each in triplicate, sequenced with both pairs of primers, the expected sequences in GenBank and every read public. A primer comparison without the mocks would have reported a difference; with them, the difference can be followed to the read where it starts, ten years later, by someone else.

---

Evidence

Paper. Bradley IM, Pinto AJ, Guest JS (2016). Design and evaluation of Illumina MiSeq-compatible, 18S rRNA gene-specific primers for improved characterization of mixed phototrophic communities. Applied and Environmental Microbiology 82(19):5878-5891. doi 10.1128/AEM.01630-16.

Data. SRA SRP071862 (BioProject PRJNA314977), 124 runs from the European Nucleotide Archive, every file checked against its published checksum: 63 samples for V4 and 61 for V8-V9, of which 21 each are the seven mock communities in triplicate, 2 x 250 bases. The twelve algae's 18S sequences, GenBank KU900218 to KU900229. Which algae each community holds, from the paper: freshwater only in MC1, marine only in MC7, all twelve in MC2 to MC6. The paper's exact freshwater-to-marine ratios were not used; the check tests membership.

Pipelines and runs. pipe_2BvukJ7IAiEJ, V4, eleven library components, run run_5UNkk1xckDWt, 1,028 lines of R. pipe_cYgfEazvPnJ9, V8-V9, the same eleven, run run_UIMHHptucXqf, 1,028 lines. pipe_N3rB6qevxVNN, V4 reverse reads only, five components, run run_BaCAcO72eM14, 725 lines. Each is a deployed run of a tested program, with the version of every library it ran on recorded.

Measured. V4: 2,445,599 read pairs, 75.9 percent passed filtering, 2,900 sequence variants, 70.3 percent of reads in the final table. V8-V9: 2,604,575 read pairs, 86.1 percent passed filtering, 2,989 variants, 78.4 percent in the final table. Mock communities: V4 ten of twelve algae, V8-V9 twelve, V4 reverse reads twelve.

Credit. The study, the mock communities and the data are their authors', in work led by Jeremy Guest's group at the University of Illinois at Urbana-Champaign, with Ameet Pinto at Northeastern University. This article describes a second analysis of their data and claims nothing further.

References

[1] Bradley IM, Pinto AJ, Guest JS. Design and evaluation of Illumina MiSeq-compatible, 18S rRNA gene-specific primers for improved characterization of mixed phototrophic communities. Applied and Environmental Microbiology 2016;82(19):5878-5891. doi:10.1128/AEM.01630-16

[2] Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods 2016;13(7):581-583. doi:10.1038/nmeth.3869

[3] Guillou L, Bachar D, Audic S, Bass D, Berney C, Bittner L, et al. The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucleic Acids Research 2013;41(D1):D597-D604. doi:10.1093/nar/gks1160

[4] Murali A, Bhargava A, Wright ES. IDTAXA: a novel approach for accurate taxonomic classification of microbiome sequences. Microbiome 2018;6:140. doi:10.1186/s40168-018-0521-5

[5] Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecology 2001;26(1):32-46. doi:10.1111/j.1442-9993.2001.01070.pp.x