A method for high precision sequencing of near full-length 16S rRNA genes on an Illumina MiSeq

Extraction of microbial DNA from foot skin

DNA was extracted from skin swabs taken from the feet of three different healthy individuals. Twelve samples were taken in total. Skin swabs were collected by swabbing either the ball or heel area of the left or right foot with a rayon swab moistened in a solution of 0.15 M NaCl and 0.1% Tween 20. The swab was rubbed firmly over the skin for approximately 30 s. Swab heads were cut into bead beating tubes, and DNA was extracted from the swabs using the BioStic DNA extraction kit (Mo-Bio), as per the manufacturers instructions. DNA was quantified on a Qubit with a HS-DNA assay (Life Technologies). Ethics approval for this study was given by the University of Technology Sydney Human Research Ethics Committee (approval number 2013000170), and participants provided written consent.

Preparation of short read 16S rRNA gene libraries for Illumina sequencing

Libraries (n = 12) of the V4 region of the 16S rRNA gene were prepared for Illumina sequencing from the microbial foot skin DNA samples using a modification of a previously published method (Caporaso et al., 2012). Briefly, samples were amplified using primers based on the Caporaso design (Caporaso et al., 2012), which were modified to include 8 nt rather than 12 nt barcodes and include a barcode on both the forward and reverse primer. The V4 region was amplified from 500 pg template DNA using 10 cycles of PCR with the modified Caporaso primers (V4_forward and V4_reverse), using different barcoded primer pairs for each sample (Table S1). After removal of excess primer via a magnetic bead clean up the samples were pooled and subjected to a further 20 cycles of PCR to enrich for amplicons containing the Illumina adapters, using primers Illumina_E_1 and Illumina_E_2 (Table S1). Pooling of samples during the enrichment PCR allows for an assessment of the putative recombination rate, by examining the rate of invalid barcode combinations that occur in the final paired end sequencing data. The method for each PCR reaction is described in detail below.

PCRs were carried out with a Taq core PCR kit (Qiagen), under the following conditions. For the initial 10 cycle PCR, reactions contained 1 × PCR buffer, 1 × Q solution (Qiagen), 250 μM dNTPs, 0.5 μM each of V4_forward and V4_reverse barcoded primers, 500 pg DNA template, and 1.25 U Taq DNA polymerase in a 50 μl reaction volume. Thermal cycling was carried out at 95 °C for 2 min, followed by 10 cycles of 95 °C for 15 s, 50 °C for 30 s and 72 °C for 90 s, followed by a final extension at 72 °C for 5 min. After a magnetic bead clean-up using 0.8 volume of SPRIselect beads (Beckman Coulter), the cleaned PCR reactions were pooled and used as input for the second PCR reaction. This PCR contained 1 × PCR buffer, 1 × Q solution (Qiagen), 250 μM dNTPs, 0.25 μM each of Illumina_E_1 and Illumina_E_2 primers (see Table S1), 31 μl pooled PCR products from the first PCR, and 1.25 U Taq DNA polymerase in a 50 μl reaction volume. Thermal cycling was carried out at 95 °C for 2 min, followed by 20 cycles of 95 °C for 15 s, 55 °C for 30 s and 72 °C for 90 s, followed by a final extension at 72 °C for 5 min. These PCR reactions were again cleaned via a magnetic bead clean-up as above and run on an Agilent Bioanalyzer using a HS-DNA chip to confirm the amplicon size and determine the concentration. Negative control PCRs were included at all stages, and all PCR products were discarded if there was any evidence of a product in the negative controls.

The short read 16S rRNA libraries were sequenced using a Nano flow cell and a 500 cycle V2 kit on an Illumina MiSeq, using custom primers as described previously (Caporaso et al., 2012). This method will be referred to as “short sequencing” and data produced with this method as “V4” data. Read pairs were merged with FLASH (Magoc & Salzberg, 2011) and de-multiplexed using a new module implemented in a previously published version of PhyloSift (Darling et al., 2014).

Sequencing of near full-length 16S rRNA gene sequences on the Illumina MiSeq

We present a method to sequence near full-length 16S rRNA gene amplicons using Illumina technology. The technique incorporates randomized molecular tags on both ends of individual 16S rRNA gene template molecules prior to PCR amplification. Copies of the templates are fragmented and sequenced and the dual tag information is used to accurately re-assemble near full-length 16S rRNA gene sequences. An overview of the method is shown in Fig. 1, and a detailed description is provided below.

Preparation of near full-length 16S rRNA gene libraries for Illumina sequencing with unique molecular tags on both ends

Primers for amplification of the 16S rRNA gene contained the 27F (Weisburg et al., 1991) or 1391R (Turner et al., 1999) bacterial primer sequences, an 8 nt barcode sequence, a 10 nt random tag and partial Illumina PE adapter sequences (Fig. 2; Table S1). The use of a 10 nt random tag on both forward and reverse primers (∼1 million possible unique tags at each end, ∼1 trillion combinations) allowed us to uniquely tag each 16S rRNA gene molecule in our pool, by modifying previously described tagging approaches (Faith et al., 2013; Lundberg et al., 2013). Template DNA was subjected to one cycle of annealing and extension with the forward primer (long_forward, Table S1), followed by a magnetic bead clean up to remove excess primer, then another cycle of annealing and extension with the reverse primer (long_reverse, Table S1), followed by another magnetic bead clean up. The first PCR carries out extension of the 16S rRNA gene from the forward primer, which uniquely tags different 16S rRNA gene templates in the reaction. The second PCR uses extension products from the first PCR as a template to produce molecules with unique tags at both ends. While the original 16S rRNA gene molecules may also act as a template in the second PCR reaction, these products will only contain an Illumina adapter at one end and will therefore not be amplified in the enrichment PCR. The enrichment PCR (34 cycles) amplifies the tagged 16S rRNA gene molecule pool, using primers that are complementary to the Illumina adapter sequences at the ends of each tagged 16S rRNA gene molecule (primers PE_1 and PE_2, Table S1).

PCRs were carried out using the Taq PCR core kit (Qiagen), and differently barcoded primers were used for each sample. Reactions contained approximately 500 pg DNA template, 0.25 μM long_forward primer, 250 μM dNTPs, 1 × PCR buffer, 1 × Q solution, and 1.25 U Taq polymerase in a 50 μl volume. Cycle conditions were 95 °C for 1 min, 50 °C for 2 min then 72 °C for 3 min. This allows extension of the 16S rRNA gene from the forward primer, which uniquely tags the forward end of each 16S rRNA gene molecule in the reaction. PCR reactions were then subjected to a magnetic bead clean up using 0.6 volumes of SPRIselect beads (Beckman Coulter) as per the manufacturer’s instructions, except that the DNA was eluted in 35 μl nuclease free water. The second PCR was set up as described above, except that 0.25 μM of the long_reverse primer was used and the template was 31 μl of the bead-cleaned first round annealing and extension reaction. Only 31 of the 35 μl of bead cleaned first round PCR was used, to prevent contamination with the magnetic beads. Cycling conditions were as in the previous step: 95 °C for 1 min, 50 °C for 2 min and 72 °C for 3 min. During this second reaction, the uniquely tagged extension products from the first reaction act as the template to produce 16S rRNA gene molecules with unique tags on the forward and reverse ends. This was followed by another magnetic bead clean up, as described above, and the output of this step was used as a template for the final PCR reaction. The final enrichment PCR reaction contained 0.5 μM of each PE_1 and PE_2 primers, 250 μM dNTPs, 1 × PCR buffer, 1 × Q solution, 31 μl template (from the bead clean up) and 1.25 U Taq polymerase in a 50 μl volume. Cycling conditions were 95 °C for 2 min, followed by 34 cycles of 95 °C for 1 min, 58 °C for 30 s, 72 °C for 2 min, and a final extension of 72 °C for 5 min. PCRs were again subjected to a magnetic bead clean up as described above, before being analysed using a high-sensitivity DNA chip on a Bioanalyser (Agilent) to determine amplicon size and concentration. Negative control PCRs were included at all stages, and all PCR products were discarded if there was any evidence of a product in the negative controls. This was assessed via Bioanalyser traces from HS-DNA chips, although we acknowledge it is possible that products below the limit of detection may still have been present.

Reconstructing full-length 16S rRNA gene sequences from tagged Illumina reads

Sequencing produced data from two kinds of fragments, those that span the entire 16S rRNA gene (end+end fragments) and those that contained either the forward or reverse end of the gene at one end with a region in the middle of the gene at the other end (end+internal fragments). Sequences from end+end fragments encoded a pairing of molecular tags and sample barcodes. Sequences were assigned to bins of original 16S rRNA gene progenitor molecules via the unique tags at either end of the molecule and were re-assembled to provide near full-length 16S rRNA gene sequences. Figure 3 shows an overview of the process.

To assign sequences to samples, the two 8 nt sample barcode regions were matched against the collection of known sample barcodes with up to one mismatch tolerated in each 8 nt barcode. Because internal regions of the 16S rRNA gene sequence might match a sample barcode, all reads with a potential sample barcode match were then screened for the presence of the proximal or distal 16S rRNA gene primer annealing sequence downstream from the sample barcode. Reads lacking a known sample barcode or the primer annealing sequence in one end were presumed to derive from an end+internal fragment.

Consensus molecular tags and elimination of recombinants

Due to sequencing error, the reads derived from the same template molecule may have had slightly different 10 nt random tagging sequences. To estimate the original 10 nt random tag sequences of tagged template molecules the UCLUST (Edgar, 2010; Edgar, 2013) algorithm was applied to identify clusters of matching random tag sequences at > 89% identity (e.g. one out of 10 bases mismatch) and to report the consensus sequences of these clusters. Clusters of molecular tags in the end+end fragments (the clustered sequences consisting of both 10 nt random tags, both 8 nt sample barcodes, and the first 14 nt of the 16S rRNA gene amplicon sequence in each read) were first identified. This was followed by the identification of the highest abundance cluster with the same combination of 10 nt random tags (one from either end) and discarding of any cluster containing one or both 10 nt random tags that were found in a different, more abundant cluster. This step aimed to identify and discard combinations of molecular tags that arose due to in vitro recombination. Recombinant forms are likely to be at lower abundance than the parental templates.

Finally, molecular tags from the entire set of reads (end+end and end+internal) were matched against the collection of consensus sequences and the reads were grouped into clusters for later assembly.

Assembly of read clusters

Each read cluster contained reads that, with high probability, originated from the same template molecule. A de novo assembly algorithm was applied on the read cluster to reconstruct as much of the original template molecule as possible. The reads were assembled using a version of the A5 pipeline called A5-miseq (Coil, Jospin & Darling, 2015) that has been modified to support assembly of reads up to 500 nt long and to trim out adapter sequence from reads instead of discarding reads containing adapter sequence. Only the first two stages of the A5-miseq pipeline were applied, involving adapter trimming, quality trimming, error correction, and contig assembly.

This method of amplifying, sequencing and assembling 16S rRNA gene sequences will be referred to as “Long-16S” and data produced with this method as “Long-16S” data, from here on in.

Analysis of V4 and Long-16S data

Both V4 and Long-16S data generated from the 12 skin samples were analysed using the software package QIIME (Caporaso et al., 2010). For comparison, the corresponding V4 region was extracted from the Long-16S sequences (which we will refer to as extracted-V4). Only those Long-16S sequences that were > 1,300 nt in length were included in the downstream analysis. V4 sequences were initially quality filtered using the default settings, with the exception of sequence length, which was altered to remove sequences less than 240 nt and longer than 260 nt. V4 sequences were additionally checked for the presence of chimeras using the UCHIME (Edgar et al., 2011) method, both against a reference database (RDP Classifier 16S training set No. 9, accessed from https://sourceforge.net/projects/rdp-classifier/files/RDP_Classifier_TrainingData/), as well as using the dataset itself as the reference. Long-16S sequences were quality filtered using default settings and excluding sequences longer than 1,400 nt. Quality filtered sequences from the V4, Long-16S and extracted-V4 datasets were then combined, and sequences were assigned to OTUs using the closed reference picking method, which assigns sequences to pre-clustered OTUs at 97% similarity from the chimera filtered Greengenes database (DeSantis et al., 2006; McDonald et al., 2012). OTUs with less than two sequences were filtered from the OTU table. Taxonomy was determined based on membership to the database of pre-clustered OTUs, and the relative abundance of taxa at different levels was generated using the summarize_taxa.py script. Differences in abundance of taxa at phylum and genera level between the V4 and Long-16S data were tested for significance using the group_significance.py script with taxa summary biom tables as input and using the Kruskal-Wallis test with the Benjamin Hochberg FDR correction for multiple testing.

Assessment of bias reduction using unique molecular tags

The use of molecular tagging has previously been shown to reduce the effect of PCR bias in RNA-seq data, for better quantitative assessment of sequences from the original samples (Islam et al., 2014). Assuming that each uniquely tagged 16S rRNA gene molecule from our skin samples was present at the same abundance as all other uniquely tagged molecules (i.e. one copy of each) and that unbiased amplification would result in an equal abundance of each cluster, we estimated the amount of biased amplification that occurred during PCR by comparing the differences in the abundance of end+end sequence clusters.

Near full-length 16S rRNA gene sequences from an Illumina MiSeq generated by molecular tagging

Sequencing of both the full-length and the tagmented amplicon pools was successful with a cluster density of 400–500 k/mm² and 832,293 read pairs. Higher clustering and sequencing output may be possible by optimising the loading concentration and ratio of end+end and end+internal pools. Clustering of end+end reads resulted in 5,085 clusters. Of these, 2,265 (44.6%) were deemed to be putative recombinant clusters, with predicted parental templates on average 29 times more abundant than putative recombinants (Fig. 4). Putative recombinant end+end sequences represented 4,378 of the total 42,715 sequences in the end+end read pool, indicating an average recombination rate of 10.2% among all samples. After binning and assembly of end+end and end+internal read clusters, 2,304 16S rRNA gene sequences were assembled from 558,053 Illumina read pairs. Sequence lengths ranged from 449 to 1,372 nt (full-length), and ∼70% of these were greater than 1,300 nt. The assembly of sequences with less than the expected length (i.e. those 400–1,300 bp) is possibly due to a lack of coverage across the internal regions for some end+end clusters. The range of sequence lengths generated is shown in Fig. 5.

Assembled sequences had consistently high quality scores across their length, with the per site average estimated PHRED quality scores at each position ranging from 54.0–89.5 (median 68.0) (Fig. 6A). This indicated estimated base-calling accuracies of greater than 99.999% at each position of the assembled 16S rRNA gene sequences. We note that errors due to base misincorporations that occur during early cycles of the enrichment PCR cannot be directly measured with this method, therefore these estimates of consensus accuracy may overestimate the true accuracy of the reconstructed 16S rRNA genes. Higher qualities for the Long-16S sequences were associated with higher coverage, which is particularly apparent at each end of the reconstructed sequences (up to 200 and beyond 1,200 nt), which were associated with one read from every read pair (end+internal or end + end) in the data set (Fig. 6B).

Short sequencing of the 16S rRNA gene V4 region

A total of 296,864 paired end 16S rRNA gene V4 sequences were generated from the 12 skin samples. Of these sequences, 11,240 could not be assigned to a sample due to invalid forward and reverse barcode combinations (e.g. combinations which were never assigned to a sample), indicating an in-vitro recombination rate of at least 3.8%, although a small proportion of this may be due to cross-contamination of barcoded primers (Kircher, Sawyer & Meyer, 2012). These sequences were removed from the dataset. We note that in-vitro recombination could also create barcode combinations that would match a valid sample and therefore be undetectable recombination events. In contrast, when attempting to detect recombination products using the chimera detection software UCHIME (as implemented in QIIME), only 0.05% of the sequences were flagged as chimeric when compared against a reference database (SILVA) and 0.2% when using the dataset itself as the reference. This highlights the difficulties of using software alone to detect recombination products from PCR in the absence of sample barcode and molecular tag information. Sequences that were flagged as chimeric using UCHIME, which had not been identified as chimeric based on sample barcode combinations (as described above) were also removed from the dataset.

Assembled near full-length 16S rRNA gene sequences produce data consistent with short read sequencing

Taxonomy, as assigned in QIIME, was similar to previous reports for skin communities, dominated by Firmicutes, Actinobacteria, and Proteobacteria. Long-16S and extracted-V4 OTUs showed the same broad taxonomic distribution as the V4 sequence data (Fig. 7). There was a small decrease in the representation of Firmicutes and an increase in the representation of Actinobacteria and Proteobacteria (Fig. 7), however these differences were not significant (Kruskal-Wallis with Benjamin-Hochberg FDR correction for multiple testing, p > 0.05). Similar taxonomic assignments between the different sequencing methods were also observed at the level of genus, with communities dominated by Staphylococcus, followed by Corynebacterium, Enhydrobacter and Acinetobacter. The Corynebacterium genus had an increased representation in the full-length data set as compared to the V4 data, which likely accounts for the observed difference in representation for the Actinobacteria phyla, but as above, this difference was not significant (Kruskal-Wallis with Benjamin-Hochberg FDR correction for multiple testing, p > 0.05).

Assessment of bias reduction using unique molecular tags

We estimated the amount of biased amplification that occurred during PCR by comparing the differences in the abundance of end+end sequence clusters. The average abundance was calculated from all clusters, and the relative mean error was 2.08 or 1.81 if singleton clusters (possible recombinants) were excluded. This indicates a standard deviation of approximately two times the average across the dataset under the particular amplification conditions used here. Figure 8 shows the distribution of the estimated amplification bias, which ranges from 0.06 to ∼32 times the average cluster abundance. This potential bias is eliminated by considering each assembled sequence cluster as having a count of one.