Microbiota profiling with long amplicons using Nanopore sequencing: full-length 16S rRNA gene and the 16S-ITS-23S of the rrn operon

Samples and DNA extraction

We first sequenced a pure bacterial isolate of S. pseudintermedius obtained from the ear of a dog affected by otitis.

Then, we used two DNA mock communities as simple and well-defined microbiota samples:

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

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

For the first PCR, we targeted: i) the full-length 16S rRNA gene using 16S-27F²⁸ and 16S-1492R²⁹ primer set and ii) the 16S-ITS-23S of the rrn operon using 16S-27F and 23S-2241R²⁸ primer set (Supplementary Table 1). All the three primers contained the Oxford Nanopore tag, which is an overhang that allows barcoding the samples during the second PCR.

PCR mixture for the full-length 16S rRNA gene (25 μl total volume) contained 5 ng of DNA template (or 2.5 μl of unquantifiable initial DNA), 1X Phusion® High Fidelity Buffer, 0.2 mM of dNTPs, 0.4 μM of 16S-27F, 0.8 μM of 16S-1492R and 0.5 U of Phusion® Hot Start II Taq Polymerase (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 16S-ITS-23S of the rrn operon (50 μl total volume) contained 5 ng of DNA template (or 2.5 μl of unquantifiable initial DNA), 1X Phusion® High Fidelity Buffer, 0.2 mM μl dNTPs 1 μM each primer and 1 U Phusion® Hot Start II Taq Polymerase. 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 16-ITS-23S of the rrn operon, respectively. Then, they were quantified using Qubit™ fluorometer (Life Technologies, Carlsbad, CA) and the volume was adjusted to begin the second round of PCR with 0.5 nM of the first PCR product or the complete volume when not reaching the required DNA mass (mostly in the samples that amplified with the 16S-ITS-23S genetic marker).

PCR mixture for the barcoding PCR (100 μl total volume) contained 0.5 nM of the first PCR product (50 ng for the 16S rRNA gene and 142 ng for the 16S-ITS-23S), 1X Phusion® High Fidelity Buffer, 0.2 mM μl dNTPs, and 2 U Phusion® Hot Start II Taq Polymerase. 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. The samples with higher DNA concentrations were checked by agarose gel to see the size profile of the PCR products (Supplementary Figure 1).

The different barcoded samples were pooled in equimolar ratio to obtain a final pool (1,000–1,500 ng in 45 μl) to do the sequencing library. In few cases, 16S-ITS-23S amplicons did not reach the initial amount of required DNA and we proceeded with lower input material.

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, as well as 45 extra base pairs from each end that correspond to the length of the universal tags and custom primers (See Supplementary Figure 2 for a schematic overview of the process and Supplementary File 1 for the bioinformatics workflow of the mapping approach).

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 16S-ITS-23S of the rrn operon. Afterwards, we removed chimeras with the following approach: i) we mapped each mock community to its mock database and the complex samples to the complete rrn database using Minimap2 v2.16 (with base-level alignment and z-score set to 70)³¹; ii) chimeras were detected and removed using yacrd v0.5³².

To assign taxonomy to the trimmed and filtered reads we used to strategies: 1) a mapping-based strategy using Minimap2 v2.16³¹ (with base-level alignment and z-score set to 70); 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 equal or larger than 1,000 bp (for 16S rRNA gene) and 3,000 bp (for the 16S-ITS-23S of the rrn operon). For reads that hit two or more references, only the alignments with the highest Smith-Waterman alignment score (AS score) were kept.

The reference databases used in this study were:

For assessing the mapping-based strategy, we have made a subset with the rrn DB to exclude all the operons that were representatives of Gammaproteobacteria class as an example to see how the alignment-based approach performs when missing main references within the database. These operons were identified by introducing a list of all the genera of the rrn DB as a batch in NCBI Taxonomy browser. The sequences belonging to Gammaproteobacteria (code 1236) were removed from the rrn DB.

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

Ampvis2 package in R was used to plot the heatmaps³⁵ and the Phyloseq package, to plot the alpha rarefaction curves³⁶.

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

Bacterial isolate analysis

We first sequenced an isolate of S. pseudintermedius obtained from a canine otitis. When using WIMP approach with the 16S-ITS-23S of the rrn operon, 97.5% of the sequences were correctly assigned at the species level as S. pseudintermedius. However, with the full-length 16S rRNA gene, 68% of the sequences were correctly assigned at the species level as S. pseudintermedius, while 13% at the genus one and ~20% were wrongly assigned, either by not reaching the species level or by giving an incorrect species (Table 1).

When using the mapping approach with the rrn DB, we obtained no hit to S. pseudintermedius. 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. This result was due to the rrn DB did not contain any representative of S. pseudintermedius. When including S. pseudintermedius sequence to the rrn DB (rrn + S. pseudintermedius) both markers retrieved the correct result with more than 97% of the assignments hitting the correct reference.

When comparing the alignment results obtained with Minimap2 and rrn DB vs rrn DB + S. pseudinteremedius, we found that the Smith-Waterman alignment score (AS) presented higher values in the correct alignments, especially when using 16S-ITS-23S marker gene (Figure 1). Thus, the AS score could be a filter to identify a wrong taxonomic assignment due to the lack of a reference in the database.

Figure 1. Violin plot representing the distribution of AS score for S. pseudintermedius isolate.

Alignment scores for each genetic marker were obtained using Minimap2 and either rrn DB or rrn DB with a reference added for S. pseudintermedius. The first two plots are for 16S-ITS-23S, whereas the second ones are for 16S rRNA.

Assessment of the mapping-based strategy

As we have already seen for S. pseudintermedius, in the mapping strategy a lack of a reference in the database led to an incorrect taxonomic assignment. Since both marker genes chosen for the microbiota profiling are highly conserved among bacteria, the mapping strategy (through Minimap2) will always align to some reference.

To check the behavior of the mapping approach when using an incomplete database, we have performed an example test using the mock communities. We have mapped the mock communities both against the complete rrn DB and against a subset of the rrn DB without any representative of the Gammaproteobacteria class. The rrn DB²⁵ contains 22,351 different bacterial species, including representatives of the species in both mock communities. We have chosen Gammaproteobacteria because each mock community contains three Gammaproteobacteria species, representing around 24% of the total microbial composition.

We checked the alignment score values and the alignment block length to detect any differences on the alignment performance when using complete or an incomplete database. We plotted two histograms: i) read counts distributed by the alignment block length; and ii) read counts distributed by the alignment score (AS) (Figure 2). For the 16S-ITS-23S genetic marker, we detected a clear pattern: when aligning to the rrn DB without Gammaproteobacteria, both histograms changed from a left-skewed distribution to a bimodal distribution with two peaks (Figure 2). A new peak appeared at the lower values that included the wrong taxonomic assignments, which are species not present in the mock community or the non-concordant hits when compared to the complete rrn DB results. Thus, for the 16S-ITS-23S marker, the initial filtering step by alignment block length will get rid of most of the incorrect taxonomic assignments. However, this pattern was not observed with the full-length 16S rRNA gene (Figure 2) or when closely-related references were present in the database, as seen above for the S. pseudintermedius isolate. So, to further confirm taxonomic results (especially for the 16S rRNA gene), we assigned taxonomy using two different bioinformatics approaches that work with different databases.

Figure 2. Assessment of the mapping-based strategy using a sample of Zymobiomics mock community.

Top part, histograms of the read counts distributed by alignment block depending on the database used (complete rrn DB vs rrn DB without Gammaproteobacteria). Alignment threshold at 3,000 bp (for 16S-ITS-23S) and 1,000 (for 16S rRNA gene) are marked with a red dashed line. Below part, histograms of the read counts distributed by AS scores depending on the database used (complete rrn DB vs rrn DB without Gammaproteobacteria).

Mock community analyses

We analyzed two microbial mock communities to validate the ability of the presented approach: i) to quantify what is expected and detect biases of the technique; and ii) to reach a reliable taxonomic assignment at the species level.

For the first aim, we used the HM-783D mock community that contained genomic DNA from 20 bacterial strains with staggered ribosomal RNA operon counts (from 10³ to 10⁶ copies per organism per µl). This mock community would allow us determining if our approach reliably represents the actual bacterial composition of the community, especially considering the low-abundant species. We assigned taxonomy with Minimap2 and a database containing only the 20 representative bacterial genomes in the HM-783D mock community (Minimap2-mock DB).

On the one hand, using the full-length 16S rRNA gene we detected all the bacterial species present in the mock community, even the low-abundant ones. On the other hand, using the 16S-ITS-23S of the rrn operon we detected only the most abundant species (at least 10⁴ operon copies) (Figure 3A). This could be due to the lower sequencing depth obtained with 16S-ITS-23S when compared with 16S rRNA (Supplementary Table 2). Moreover, the relative abundances of 16S-ITS-23S sequences were more biased than those obtained from 16S rRNA gene sequencing, which confirmed that the primers for 16S-ITS-23S of the rrn operon need to be improved for a better representation of the actual abundances (Figure 3B and 3C).

Figure 3. HM-783D mock community analysis.

(A) Heat map representing the HM-783D mock community composition when mapped to its mock database. Grey colour represents the bacteria that were not detected (<10⁴ copies with rrn operon). (B) Linear regression analysis of relative read proportions obtained using full-length 16S rRNA gene for all bacterial species present in HM-783D mock community and the actual operon copies (in log scale). (C) Linear regression analysis of relative read proportions obtained using the 16S-ITS-23S genetic marker for all bacterial species present in HM-783D mock community and the actual operon copies (in log scale).

To assess if the technique and the analyses would give a reliable taxonomy at the species level we used Zymobiomics mock community, which contains equal quantities of 8 bacterial species. The expected 16S rRNA gene content for each representative is also known, so we were able to determine if the different analysis approaches reliably represented the actual bacterial composition of the community. We sequenced the Zymobiomics mock community twice per marker gene and found that the replicates presented equivalent results.

We assigned taxonomy with three different approaches: i) Minimap2 and a database containing only the 8 bacterial species of the correspondent mock community (mock DB); ii) Minimap2 and a database containing sequences for the rrn operon of 22,351 different bacterial species (rrn DB (25)) and iii) WIMP from EPI2ME and NCBI RefSeq database.

Similarly to what we have seen for the HM-783D mock community, when using the mapping strategy with Minimap2 and the mock database, we detected that the full-length 16S rRNA gene retrieved better the actual abundances of the mock community. The 16S-ITS-23S genetic marker over- and underrepresented most of the bacterial species in the mock community. When using larger databases such as rrn DB and NCBI RefSeq, both the full-length 16S rRNA gene and the 16S-ITS-23S were able to detect 8 out of 8 bacterial species of the Zymobiomics mock community (Figure 4A). However, we also detected other taxa that included mostly higher taxonomic rank taxa (sequences not assigned to species level), but also not expected taxa (wrongly-assigned species), especially with WIMP and NCBI RefSeq database (see Supplementary Table 3 for complete taxonomic assignments).

Figure 4. Zymobiomics mock community taxonomic analysis and diversity.

(A) Heat map representing the relative abundance of the Zymobiomics mock community. “REF” column represents the theoretical composition of the mock community regarding the 16S rRNA gene content of each bacterium. (B) Alpha diversity rarefaction plot using observed taxa metrics.

On one hand, the mapping approach using the rrn DB provided highly similar results to the reference, despite the larger size of this database (Figure 4A), especially with 16S-ITS-23S marker (99% of the reads were correctly assigned at the species level with 16S-ITS-23S; near 90% for 16S). On the other hand, with the WIMP workflow and NCBI RefSeq database, a larger number of sequences are classified as “Other taxa”. Again, this is especially remarkable when using the full-length 16S rRNA gene, with > 50% of the taxonomic assignments not hitting the expected bacterial species. The results were also confirmed by alpha diversity analyses: WIMP strategy overestimated the actual bacterial diversity, when compared to rrn DB and the reference (Figure 4B).

Complex microbial community analyses

We profiled two complex and uncharacterized microbial communities from dog skin (chin and dorsal). We used both long-amplicon markers and the two bioinformatics approaches –Minimap2 and rrn DB and WIMP with NCBI RefSeq database– to corroborate the results.

For chin samples of healthy dogs, we found a high abundance of Pseudomonas species (>40% of total relative abundance using 16S rRNA and >60% using 16S-ITS-23S) 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 5A and Supplementary Table 3). 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, with newer nomenclature (Figure 5B and Supplementary Table 3). For both skin sample replicates, the results of the most abundant species converged using the two different methods and allowed for characterizing this complex low-biomass microbial community at the species level.

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

(A) Chin samples: heat map representing the relative abundance of the main bacterial species in chin samples using WIMP and Minimap2. (B) Dorsal skin samples: heat map representing the relative abundance of the main bacterial species using WIMP and Minimap2. The lower heat map represents the remaining contamination from the previous run using HM-783D mock community within the same flowcell. Samples marked with a red line shared barcode with the mock community.

Finally, analyzing the dorsal skin samples, we also detected the presence of contamination from the previous nanopore run (Supplementary Table 2). 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. We were able to detect mock community representatives within the re-used barcode (Figure 5B). 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.

De novo assembly of the 16S-ITS-23S genetic marker

To further confirm the results obtained with the complex samples using directly the raw reads, we performed the assembly and consensus of the 16S-ITS-23S genetic marker using canu and we assigned taxonomy of the consensus sequences using BLAST.

For the HM-783D mock community, we were able to retrieve some of the most abundant bacterial species blasting with >99% of identity to their reference (Escherichia coli, Staphylococcus epidermidis, Streptococcus mutans and Clostridium beijerinnkii). For the Zymobiomics mock community, we found a consensus sequence for all the bacterial species with >99% of identity. The only two exceptions were Salmonella enterica that presented a consensus sequence with an identity of 98.8% and Escherichia coli that presented no consensus sequence (Table 2 and Supplementary Table 4).

For the complex microbial communities, the de novo consensus sequences usually presented lower identity than those of the mock community. For the dorsal skin, the ones with higher consensus accuracy were Stenotrophomonas rhizophila (99.3 %), Streptococcus mutans (99.1 %) and Arthrobacter arilaitensis (98.6 %), thus confirming the results retrieved using directly the raw reads. For the chin, the three contigs with higher consensus accuracy hit Pseudomonas fluorescens (>98.5%). We detected also other contigs with lower consensus accuracy values, with previously not seen taxonomic assignments (Supplementary Table 4).

Underlying data

The datasets analyzed during the current study are available in the NCBI Sequence Read Archive, under the Bioproject accession number PRJNA495486.

Extended data

All the supplementary data has been added in an OSF repository (doi: http://doi.org/10.17605/OSF.IO/8MYKV)²². We provide here a complete list:

-  Supplementary Table 1. Primer sequences for amplifying the full-length 16S rRNA gene and 16S-ITS-23S of the rrn operon.

-  Supplementary Table 2. Samples included in the study, run summary and quality control results. *For mock communities, the mock DB was used. For complex communities, the rrn DB.

-  Supplementary Table 3. Taxonomic assignments table of each sample with the different approaches.

-  Supplementary Table 4. De novo results obtained with canu and their taxonomic assignment using BLAST.

-  Supplementary Figure 1. Photo of the agarose gel electrophoresis of some of the samples.

-  Supplementary Figure 2. Main workflow overview. Detailed bioinformatics workflow can be found in Supplementary File 1.

-  Supplementary File 1. Bioinformatics workflow used for the mapping approach.