Microbiota profiling with long amplicons using Nanopore sequencing: full-length 16S rRNA gene and whole rrn operon

Samples and DNA extraction

We used two DNA mock communities as simple, well-defined microbiota samples:

We also sequenced a pure bacterial isolate of Staphylococcus pseudintermedius obtained from the ear of a dog affected with otitis.

As a complex microbial community, we used two DNA sample pools from the skin microbiota of healthy dogs targeting two different skin sites: i) dorsal back (DNA from two dorsal samples from Beagle dogs); and ii) chin (DNA from five chin samples from Golden Retriever/Labrador crossed dogs). Skin microbiota samples were collected using Sterile Catch-All™ Sample Collection Swabs (Epicentre Biotechnologies) soaked in sterile SCF-1 solution (50 mM Tris buffer (pH 8), 1 mM EDTA, and 0.5% Tween-20). DNA was extracted from the swabs using the PowerSoil™ DNA isolation kit (MO BIO) and blank samples were processed simultaneously (for further details on sample collection and DNA extraction see 27).

PCR amplification of ribosomal markers

There were two ribosomal markers evaluated in this study: full-length 16S rRNA gene (~1,500 bp) and the whole rrn operon (~4,500 bp). Before sequencing, bacterial DNA was amplified using nested PCR, with a first PCR to add the specific primer sets (Table 1) tagged with the Oxford Nanopore universal tag and a second PCR to add the barcodes from the barcoding kit (EXP-PBC001). Each PCR reaction included a no-template control sample to assess possible reagent contamination.

For the first PCR, we targeted the full 16S rRNA gene using 16S-27F and 16S-1492R primer set and the whole rrn operon (16S rRNA gene–ITS–23S rRNA gene) using 16S-27F and 23S-2241R primer set (Table 1).

PCR mixture for full-length 16S rRNA gene (25 μl total volume) contained 5 ng of DNA template, 5 μl of 5X Phusion® High Fidelity Buffer, 2.5 μl of dNTPs (2 mM), 1 μl of 16S-27F (0.4 μM), 2 μl of 16S-1492R (0.8 μM) and 0.25 μl of Phusion® Hot Start II Taq Polymerase (0.5 U) (Thermo Scientific, Vilnius, Lithuania). The PCR thermal profile consisted of an initial denaturation of 30 s at 98°C, followed by 25 cycles of 15 s at 98°C, 15 s at 51°C, 45 s at 72°C, and a final step of 7 min at 72°C.

PCR mixture for the rrn whole operon (50 μl total volume) contained 5 ng of DNA template, 10 μl 5X Phusion® High Fidelity Buffer, 5 μl dNTPs (2 mM), 5 μl each primer (1 μM) and 0.5 μl Phusion® Hot Start II Taq Polymerase (1 U). The PCR thermal profile consisted of an initial denaturation of 30 s at 98°C, followed by 25 cycles of 7 s at 98°C, 30 s at 59°C, 150 s at 72°C, and a final step of 10 min at 72°C.

The amplicons were cleaned-up with the AMPure XP beads (Beckman Coulter) using a 0.5X and 0.45X ratio for the 16S rRNA gene and the whole rrn operon, respectively. Then they were quantified using Qubit™ fluorometer (Life Technologies, Carlsbad, CA) and volume was adjusted to begin the second round of PCR with 0.5 nM of the first PCR product or the whole volume when not reaching the required concentration (mostly for samples that amplified the rrn operon).

PCR mixture for the barcoding PCR (100 μl total volume) contained 0.5 nM of first PCR product, 20 μl 5X Phusion® High Fidelity Buffer, 10 μl dNTPs (2 mM), and 1 μl Phusion® Hot Start II Taq Polymerase (2 U). Each PCR tube contained the DNA, the PCR mixture and 2 μl of the specific barcode. The PCR thermal profile consisted of an initial denaturation of 30 s at 98°C, followed by 15 cycles of 7 s at 98°C, 15 s at 62°C, 45 s (for the 16S rRNA gene) or 150 s (for rrn operon) at 72°C, and a final step of 10 min at 72°C.

Again, the amplicons were cleaned-up with the AMPure XP beads (Beckman Coulter) using a 0.5X and 0.45X ratio for the 16S rRNA gene and the whole rrn operon, respectively. For each sample, quality and quantity were assessed using Nanodrop and Qubit™ fluorometer (Life Technologies, Carlsbad, CA), respectively.

In most cases, the different barcoded samples were pooled in equimolar ratio to obtain a final pool (1000–1500 ng in 45 μl) to do the sequencing library.

Nanopore sequencing library preparation

The Ligation Sequencing Kit 1D (SQK-LSK108; Oxford Nanopore Technologies) was used to prepare the amplicon library to load into the MinION^TM (Oxford Nanopore Technologies), following the manufacturer’s protocol. Input DNA samples were composed of 1–1.5 μg of the barcoded DNA pool in a volume of 45 μl and 5 μl of DNA CS (DNA from lambda phage, used as a positive control in the sequencing). The DNA was processed for end repair and dA-tailing using the NEBNext End Repair/dA-tailing Module (New England Biolabs). A purification step using 1X Agencourt AMPure XP beads (Beckman Coulter) was performed.

For the adapter ligation step, a total of 0.2 pmol of the end-prepped DNA were added in a mix containing 50 μl of Blunt/TA ligase master mix (New England Biolabs) and 20 μl of adapter mix and then incubated at room temperature for 10 min. We performed a purification step using Adapter Bead Binding buffer (provided in the SQK-LSK108 kit) and 0.5X Agencourt AMPure XP beads (Beckman Coulter) to finally obtain the DNA library.

We prepared the pre-sequencing mix (14 μl of DNA library) to be loaded by mixing it with Library Loading beads (25.5 μl) and Running Buffer with fuel mix (35.5 μl). We used two SpotON Flow Cells Mk I (R9.4.1) (FLO-MIN106). After the quality control, we primed the flowcell with a mixture of Running Buffer with fuel mix (RBF from SQK-LSK108) and Nuclease-free water (575 μl + 625 μl). Immediately after priming, the nanopore sequencing library was loaded in a dropwise fashion using the SpotON port.

Once the library was loaded, we initiated a standard 48 h sequencing protocol using the MinKNOW™ software v1.15.

Data analysis workflow

The samples were run using the MinKNOW software. After the run, fast5 files were base-called and de-multiplexed using Albacore v2.3.1. A second de-multiplexing round was performed with Porechop v0.2.3³⁰, where only the barcodes that agreed with Albacore were kept. Porechop was also used to trim the barcodes and the adapters from the sequences (Figure 1).

Figure 1. Bioinformatics analysis workflow from raw reads to final data.

Moreover, we removed 45 extra base pairs from each end that correspond to the length of the universal tags and custom primers. After the trimming, reads were selected by size: 1,200 bp to 1,800 bp for 16S rRNA gene; and 3,500 to 5,000 bp for the rrn operon. We mapped the sequences obtained to the rrn database using Minimap2 v2.9³¹. Afterwards chimeras were detected and removed using yacrd v0.3³².

To assign taxonomy to the trimmed and filtered reads we used to strategies: 1) a mapping-based strategy using Minimap2³¹; or 2) a taxonomic classifier using What’s in my Pot (WIMP)³³, a workflow from EPI2ME in the Oxford Nanopore Technologies cloud (based on Centrifuge software³⁴).

For the mapping-based strategy, we performed Minimap2 again with the non-chimeric sequences. We applied extra filtering steps to retain the final results: we kept only those reads that aligned to the reference with a block larger than 1,000 bp (for 16S rRNA gene) and 3,000 bp (for the whole rrn operon). For reads that hit two or more references, only the alignments with the highest Smith-Waterman alignment score were kept. After filtering, the multimapping was mostly present in cases with entries that belonged to the same taxonomy.

The reference databases used in this study were:

For the taxonomic classification using the WIMP workflow, which uses the NCBI database, only those hits with a classification score >300 were kept³⁴.

An earlier version of this article can be found on bioRxiv (doi: https://doi.org/10.1101/450734)

Quality filtering results

After Albacore basecalling and Porechop processing, we lost around 5% of the initial reads (3-13%). After length trimming step, we lost more sequences (Table 2). In general, the samples amplified using 16S rRNA marker gene recovered a higher percentage of reads after the quality control when compared to the rrn operon: 74–95% vs. 32–80%. Particularly for rrn operon, the largest percentage of reads was lost during the length trimming step: some of the reads included in that barcode presented the length of the 16S rRNA gene.

After this first quality control, we performed an alignment with the mock and the rrn databases and checked for chimeras. Chimeras detected were dependent on the database used for the alignment. As a positive control, we used mock samples with their mock database. Chimera ratio was higher for 16S rRNA gene amplicons (around ~40%) than for rrn operon (~10%), suggesting that PCR conditions for the 16S rRNA gene need to be adjusted or the PCR cycles reduced.

To conclude, the final useful sequences when amplifying for either amplicon were ~40%. In 16S rRNA gene, sequences were lost in the chimera checking step. In the rrn operon, sequences were lost in the length trimming step, probably due to the underrepresentation of the amplicon in the flowcell, since we ran them together with full-length 16S rRNA amplicons in the same flow-cell.

Mock community analyses

Microbial Mock Community HM-783D contained genomic DNA from 20 bacterial strains with staggered ribosomal RNA operon counts (from 1,000 to 1,000,000 copies per organism per µl). The bacterial composition detected should be proportional to the operon counts. This mock community would allow us determining if our approach reliably represents the actual bacterial composition of the community, especially considering low-abundant species.

We analyzed HM-783D mock community against its own database, which contains only the 20 representative species. On the one hand, using 16S rRNA gene we were able to detect all the bacterial species present in the mock community, even the low-abundant ones. On the other hand, using the rrn operon we were able to detect only the most abundant species (at least 10⁴ operon copies) (Figure 2). This could be due to the lower sequencing depth obtained with rrn when compared with 16S rRNA, and probably due to the underrepresentation of the rrn amplicon in the flowcell when running together with the full-length 16S rRNA amplicons in the same flow-cell, as detailed above. Moreover, the relative abundances of rrn operon sequences were more biased than those obtained from 16S rRNA gene sequencing, when compared to those expected, which confirmed that the primers for rrn need to be improved for universality.

Figure 2. Heat map representing the HM-783D mock community when mapped to its mock database.

The darkest blue represents the bacteria that were not detected (<10⁴ copies with rrn operon), whereas the darkest red represents the most abundant bacteria.

Zymobiomics mock community presents the same amount of genomic DNA from 8 different bacterial species; the expected 16S rRNA gene content for each representative is also known, so we are able to determine if our approach represents the actual bacterial composition of the community reliably.

Both 16S rRNA gene and rrn operon sequencing were able to detect 8 out of 8 bacterial species for the Zymobiomics mock community, using Minimap2 and WIMP. The “Other taxa” group in Figure 3A can indicate: i) not expected taxa (wrongly-assigned species, or previous contamination); or ii) higher taxonomic rank taxa (sequences not assigned to species level).

Figure 3. Zymobiomics mock community taxonomic analysis and diversity.

(A) Bar plots representing the relative abundance of the Zymobiomics mock community. “REF” bar represents the theoretical composition of the mock community regarding the 16S rRNA gene content of each bacterium. (B) Alpha diversity rarefaction plot using the rrn database at the species level. (C) WIMP output represented in a taxonomic tree and percentage of reads in each taxonomic rank.

Using the mock community database (that contains only the 8 members of that community), we aimed to assess the biases regarding the actual abundance profile. 16S rRNA gene better represented the bacterial composition of the mock community, when considering the abundances. The rrn operon amplification over-represented Escherichia coli and Staphylococcus aureus and under-represented Enterococcus faecalis.

The rrn database²⁵ contains 22,351 different bacterial species, including representatives of the species in the mock community. When using the rrn database, we found that the rrn operon was a better marker than 16S rRNA: more than 98% of the sequences mapped to the corresponding species, and only <2% of the total sequences mapped to a wrong species with the rrn operon, whereas ~15% of the sequences were given the wrong taxonomy with 16S rRNA. We performed alpha diversity analyses using the same rrn database. The rrn operon hit 26 different species, whereas 16S rRNA over-estimated the actual diversity, with hits to 202 different species (Figure 3B). However, when considering abundances, the diversity values are more similar, with Shannon indices of 1.95 and 2.51, when using rrn operon and 16S rRNA, respectively (at 30,000 sequences/sample).

Using WIMP, we confirmed again the higher resolution power of rrn operon: ~70% of the sequences were assigned to the correct species compared to ~45% for the 16S rRNA gene. Among all the bacterial species included in the mock community, Bacillus subtilis presented the most trouble for the correct taxonomic classification. The theoretically expected abundance for B. subtilis is 17% using the 16S rRNA gene. When using WIMP, only 5% of the total sequences were correctly classified at the species level, another 5% was classified correctly at the genus level, and another 10% was incorrectly classified as other Bacillus species (Figure 3C).

Apart from the mock communities, we also sequenced an isolate of Staphylococcus pseudintermedius obtained from canine otitis. When using WIMP approach with rrn operon, 97.5% of the sequences were correctly assigned to the S. pseudintermedius. However, with the 16S rRNA gene, 68% of the sequences were correctly assigned at the species level and 13% at the genus (Table 3). The wrong assigned species for rrn operon was ~2.5%, compared to ~20% for the 16S rRNA gene. On the other hand, through mapping the sequences to rrn database using Minimap2, we obtained no hit to S. pseudintermedius, since there is no representative in the rrn database. Instead, they were hitting mostly to Staphylococcus schleiferi, which is a closely related species; there were also few hits to Staphylococcus hyicus and Staphylococcus agnetis. These results highlight the need of comprehensive databases that include representatives of all the microorganisms relevant to a microbiome to correctly assign taxonomy.

Complex microbial community analyses

After the first analyses with the mock communities, we were able to detect that the taxonomic resolution was higher when using rrn operon; however, the abundance profile was more reliable using 16S rRNA marker gene. If a bacterial species is not present in the database, the mapping strategy will give us the closest sequence resulting to an inaccurate taxonomic profile, such as we have seen for the Staphylococcus pseudintermedius isolate.

Here, we aimed to taxonomically profile two complex and uncharacterized microbial communities from dog skin (chin and dorsal) using the two different markers and comparing the mapping strategy (Minimap2 and rrn database) with the WIMP workflow (NCBI database).

For chin samples of healthy dogs, we found a high abundance of Pseudomonas species followed by other genus with lower abundances such as Erwinia and Pantoea. Focusing on Pseudomonas, at the species level we were able to detect that the most abundant species was Pseudomonas koreensis, followed by Pseudomonas putida and Pseudomonas fluorescens (Figure 4A and Supplementary Table 1). On the other hand, dorsal skin samples were dominated by bacteria from the genera Stenotrophomonas, Sanguibacter, and Bacillus. We reached species level for Stenotrophomonas rhizophila and Sanguibacter keddieii. It should be noted that Glutamicibacter arilaitensis is the same species as Arthrobacter arilaitensis, but is the up-to-date nomenclature³⁵ (Figure 4B and Supplementary Table 1). For both skin sample replicates, the results of the most abundant species converged and allowed for characterizing this complex low-biomass microbial community at the species level.

Figure 4. Microbiota composition of complex communities: skin samples of healthy dogs.

(A) Chin samples: upper part of the graphic, bar plot of the composition at the genus level using WIMP (left) and Minimap2 (right); lower part, heat map of the Pseudomonas species within the community (scaled at 100%). (B) Dorsal skin samples: upper part of the graphic, bar plot of the composition at the genus level using WIMP (left) and Minimap2 (right); lower part, heat map of the ten most abundant taxa within the community. N/A, taxon was not present in the database.

Finally, analyzing the dorsal skin samples, we also detected the presence of contamination from the previous nanopore run. We sequenced dorsal skin samples twice: one with a barcode previously used for sequencing the HM-783D mock community and another one with a new barcode (Table 2). We were able to detect mock community representatives within the re-used barcode (Figure 5). Some of them were found only in the sample that was using the re-used barcode (Sample_1); others were also present in the skin sample, such as Bacillus cereus or Staphylococcus aureus. In total, this contamination from the previous run was representing ~6% of the sample composition.

Figure 5. Heat map representing the HM-783D mock community contamination.

Samples Skin_1 are the ones re-using the HM-783D barcode from the previous run within the same flowcell. Samples Skin_2 are using a new barcode (not used in a prior run of the same flowcell). Values within the heat map are relative abundances in percentages. The darkest blue represents the bacteria that were not detected and the darkest red the most abundant bacterial contaminants.