Development of a genus-specific next generation sequencing approach for sensitive and quantitative determination of the Legionella microbiome in freshwater systems
© The Author(s). 2017
Received: 11 January 2017
Accepted: 21 March 2017
Published: 31 March 2017
Next Generation Sequencing (NGS) has revolutionized the analysis of natural and man-made microbial communities by using universal primers for bacteria in a PCR based approach targeting the 16S rRNA gene. In our study we narrowed primer specificity to a single, monophyletic genus because for many questions in microbiology only a specific part of the whole microbiome is of interest. We have chosen the genus Legionella, comprising more than 20 pathogenic species, due to its high relevance for water-based respiratory infections.
A new NGS-based approach was designed by sequencing 16S rRNA gene amplicons specific for the genus Legionella using the Illumina MiSeq technology. This approach was validated and applied to a set of representative freshwater samples.
Our results revealed that the generated libraries presented a low average raw error rate per base (<0.5%); and substantiated the use of high-fidelity enzymes, such as KAPA HiFi, for increased sequence accuracy and quality. The approach also showed high in situ specificity (>95%) and very good repeatability. Only in samples in which the gammabacterial clade SAR86 was present more than 1% non-Legionella sequences were observed. Next-generation sequencing read counts did not reveal considerable amplification/sequencing biases and showed a sensitive as well as precise quantification of L. pneumophila along a dilution range using a spiked-in, certified genome standard. The genome standard and a mock community consisting of six different Legionella species demonstrated that the developed NGS approach was quantitative and specific at the level of individual species, including L. pneumophila. The sensitivity of our genus-specific approach was at least one order of magnitude higher compared to the universal NGS approach. Comparison of quantification by real-time PCR showed consistency with the NGS data. Overall, our NGS approach can determine the quantitative abundances of Legionella species, i. e. the complete Legionella microbiome, without the need for species-specific primers.
The developed NGS approach provides a new molecular surveillance tool to monitor all Legionella species in qualitative and quantitative terms if a spiked-in genome standard is used to calibrate the method. Overall, the genus-specific NGS approach opens up a new avenue to massive parallel diagnostics in a quantitative, specific and sensitive way.
For many aspects of microbiology only part of the microbiome is of specific interest. For instance, in medical and environmental diagnostics the prime interest often lays in pathogenic or potentially pathogenic microorganisms. We will call this part of the microbiome “patho-microbiome” and it can consist of viruses, bacteria and lower eukaryotes. Despite the importance of assessing the total patho-microbiome, the methodologies to address all major pathogenic microorganisms are very complex . To make such an enterprise more amenable to a straight forward methodological approach, we have chosen the genus Legionella as a subject, which is not only relevant for medical and environmental microbiologists but also taxonomically well defined.
Within the genus Legionella, several described species can cause clinical disease in humans, such as L. bozemanii, L. longbeachae, L. micdadei and L. pneumophila [2, 3]. Nonetheless, the latest, i.e. L. pneumophila, is the most common causative agent of legionellosis (about 90% of the cases worldwide) or Legionnaires’ disease, the majority of them community-acquired . The infections caused by almost all non-L. pneumophila species have been associated with freshwater systems in hospital settings . Natural aquatic environments are the major reservoir for Legionella species, with aerosol-generating man-made water systems such as cooling towers, drinking water supply systems (DWSS) and recreational waters being the main sources of human exposure to Legionella species . Exceptionally, L. longbeachae infections are commonly linked with the presence of the bacteria in potting soils, instead of freshwater environments, with a high frequency of reported cases in Australia .
In freshwater systems, monitoring of pathogenic bacteria, hazard prediction and risk assessment are of high relevance to human health. To successfully achieve this, rapid and accurate detection of clinically relevant pathogens and information on bacterial diversity in contaminated systems are needed in order to improve prevention and achieve a better pathogen control. Currently, there are several methods for detection and/or quantification of Legionella species in freshwater systems including: culture, fluorescent in situ hybridization (FISH), flow cytometry, endpoint PCR and qPCR assays . Though EU governmental agencies request cultivation as the official reference standard method, several limitations have been reported, such as the need of complex culture media and relatively long incubation times for Legionella growth; the overgrowth of competing and unwanted organisms; and the ability of Legionella species to undergo a viable but non-culturable (VBNC) state [8, 9]. To overcome these limitations, molecular assays have been developed to provide rapid, highly specific and sensitive detection and quantification of Legionella species [10, 11].
Advances in molecular anaylsis in the last decade have provided new opportunities to approach and improve diagnosis of pathogens. Next-Generation Sequencing (NGS) will be more frequently applied for the assessment of microbial water quality with the expectation to understand which part of the DWSS and which pathogens are critical to human health . Especially, sequencing of PCR amplicons of 16S rRNA genes has been very successfully used to understand community dynamics of the drinking water microbiome . Thus, high-throughput sequencing of PCR amplicons has the potential to provide valuable information regarding not only the environmental distribution and diversity of Legionella species but also the temporal and spatial behaviour of the whole Legionella microbiome. To this end, we pursued a new NGS approach for the assessment of the Legionella microbiome based on the genus-specific amplification of 16S rRNA genes using the Illumina MiSeq technology. Such an approach should be validated with environmental DNA obtained from reference water samples of several freshwater environments, such as cold and hot drinking water and water from a cooling tower. Genome standards for L. pneumophila and other Legionella species should be used to calibrate the molecular assay and determine its accuracy and precision. Along these lines, the aim of our study was to develop and validate a Legionella genus-specific Illumina-based NGS approach that allows accurate identification and precise quantification of the most important Legionella species, i. e. the Legionella microbiome. Validation of the developed approach was done in silico, in vitro and in situ, with emphasis on the latter by using a set of representative freshwater samples including a comparison with real-time PCR measurements in terms of the accuracy of quantification of Legionella genus-specific abundances by the NGS-based approach.
Freshwater systems and sampling
Cold drinking water was sampled in March and April 2009. It was taken from the tap of laboratory D0.04 on the campus of the Helmholtz Centre for Infection Research (HZI), Braunschweig, Germany, with five minute flushing to prevent stagnant water to be sampled. This drinking water originates from two surface water reservoirs situated in the Harz Mountains, 40 km south of Braunschweig. More details on the respective drinking water supply system are given elsewhere . Hot drinking water was sampled, in March and April 2009, from a shower next to laboratory D0.04 also with several minutes flushing. The protocol is described in more detail in a preceding publication . Cooling tower water was monthly sampled between January 2013 and December 2014 from a cooling tower on the campus of the HZI. The studied cooling tower is used for discharge of waste heat generated by air conditioning and other heat producing processes at the site. It is an open circuit cooling tower with antifouling coating (Gohl GmbH, Singen, Germany). Regular drinking water is used as make up water, whereas blow down is an automated process based on the conductivity of the water. Disinfection of the system was achieved by continuous silver plus hydrogen peroxide treatment.
Bulk water microorganisms were sampled by filtration according to Eichler et al. . Briefly, 5 l of drinking water and 3 l of cooling tower water were filtered through a filter sandwich consisting of a 0.2 μm pore size polycarbonate filter (90 mm diameter; Nucleopore; Whatman, Maidstone, UK) with a precombusted glass fiber filter on top (90 mm diameter; GF/F; Whatman, Maidstone, UK). Biomass harvested on filter sandwiches was stored at−70 °C for the molecular analyses.
DNA extraction of water samples
DNA was extracted from the filter sandwiches according to a modified DNeasy (Qiagen, Hilden, Germany) protocol by Henne et al. . This DNA extraction protocol consists of cutting sandwich filters into pieces and an incubation step with lysis buffer containing 10 mg/ml lysozyme (Sigma) for 60 min at 37 °C. Next, samples were digested with proteinase K digestion according to the manufacturer instructions and heated to 70 °C in a water bath for 20 min. After filtration through a polyamide mesh with 250 μm mesh size, absolute ethanol was added to the filtrate (ratio filtrate/ethanol 2:1) and the mixture was applied to the adequate spin-column of the Qiagen kit. From this point forward, the washing and elution steps were performed according to the manufacturer instructions. DNA extracts were stored at−20 °C until further molecular analyses.
Illumina MiSeq library design and preparation
NGS libraries construction and preparation were adapted from previously described protocols [15–18]. A dual-indexing strategy was chosen and primer sequences are described in the supplementary material (Additional file 1: Table S1). More specifically, the amplification of a 421 bp region of the 16S rRNA gene of Legionella including the variable regions V3 and V4 was performed using primers Lgsp17F 5’-GGCCTACCAAGGCGACGATCG-3’ (Escherichia coli 16S rRNA gene position 264–284) and Lgsp28R 5’-CACCGGAAATTCCACTACCCTCTC-3’ (E. coli 16S rRNA gene starting position 661–684), previously developed by Kahlisch et al. . The Legionella genus-specific primer pair Lgsp17F/Lgsp28R was then combined with the indexing primers for multiplexing and sequencing on the Illumina MiSeq platform. All primers were synthesized and HPSF purified by Eurofins MWG Operon, Ebersberg, Germany. The 16S rRNA gene amplicons (details of PCR conditions see next paragraph) were size-selected by gel electrophoresis on a 2% agarose gel pre-stained with GelRed Nucleic Acid Gel Stain (Biotium, Hayward, CA, USA). The bands with the expected target length were extracted over an UV light box and then recovered using the QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). QIAquick Gel Extraction Kit manufacturer’s specifications were slightly modified. Briefly, 1.5 ml of Buffer QG was added to the gel band followed by incubation at room temperature. Next, 250 μl of isopropanol (AppliChem, Darmstadt, Germany) were added, the mixture vortexed and placed in a QIAquick spin column and centrifuged for 1 min at 13,000 rpm. The remaining steps were performed according to the manufacturer instructions to obtain clean library amplicons, concentrated in 30 μl of EB buffer. DNA concentration of the extracted amplicons was assessed on a VICTOR X3 2030 Multilabel Plate Reader (Perkin Elmer, Germany) using a Quant-iT Picogreen dsDNA assay kit (Life Technologies, Oregon, USA) to ensure, when required, equimolar amounts of DNA with unique indices were pooled. Subsequently, pooled samples were purified by the MinElute PCR Purification Kit according to manufacturer protocol (Qiagen, Hilden, Germany). Molarity was quantified and library fragment size confirmed with Agilent Bioanalyzer. The resulting libraries were sequenced, in the Genome Analysis Department of the HZI, using Illumina MiSeq platform, generating paired-end 250 bp sequences.
Amplification and Sequencing of Legionella-specific 16S rRNA fragment using KAPA HiFi DNA polymerase and HotStarTaq DNA polymerase
We compared the proof reading DNA polymerase KAPA HiFI with the most widely used HotStarTaq DNA polymerase to understand the effect of the DNA polymerase on the sequencing error rate. The Legionella-specific PCR reaction, using KAPA HiFi DNA polymerase, had the following components (50 μl total volume): KAPA HiFi HotStart ReadyMix 1X (Peqlab, Erlangen, Germany), 0.3 μM of each primer and 5 ng of the extracted environmental DNA. The PCR was carried out on a Biorad Thermo cycler 96-well iCycler under the following conditions: an initial denaturation step of 95° for 5 min, followed by 30 cycles of 98 °C for 20 s, 68 °C for 20 s, 72 °C for 20 s, with a final extension at 72 °C for 5 min. A second PCR (Multiplexing) had the following components (50 μl total volume): KAPA HiFi HotStart ReadyMix 1X (Peqlab, Erlangen, Germany), 0.3 μM of each primer and 2 μl of the previous PCR amplicons. This PCR was carried out on a Biorad Thermo cycler 96-well iCycler under the following conditions: an initial denaturation step of 95 °C for 5 min, followed by 10 cycles of 98 °C for 20 s, 59 °C for 20 s, 72 °C for 20 s, with a final extension at 72 °C for 7 min.
For HotStarTaq DNA polymerase-based amplification, the PCR mix comprised 0.1 mM concentration of each dNTP (Bioline, Luckenwalde, Germany), 0.4 mM MgCl2, 1X PCR reaction buffer, 0.03 U of HotStarTaq DNA polymerase (Qiagen, Hilden, Germany), 0.4 μM of each primer and 5 ng of the extracted environmental DNA in a total volume of 50 μl. The PCR included an initial denaturation step of 95 °C for 15 min; 30 cycles of 1 min at 95 °C, 30 s at 66.5 °C and 30 s at 72 °C followed by a final extension at 72 °C for 10 min. 2 μL of the target-specific PCR product from the previous PCR reaction were used as template for a second PCR reaction (Multiplexing PCR). Additionally, the mix contained 0.1 mM concentration of each dNTP (Bioline, Luckenwalde, Germany), 0.75 mM MgCl2, 1X PCR reaction buffer, 0.03 U of HotStarTaq DNA polymerase (Qiagen, Hilden, Germany), 0.4 μM of each primer. The PCR conditions included a cycle at 95 °C for 15 min, 10 cycles of 95 °C for 45 s, 57 °C for 30 s and 72 °C for 30 s, ending with an extension at 72 °C for 10 min.
Pan-bacterial amplification and sequencing of 16S rRNA gene
We applied a commonly used PCR protocol with universal primers  for a comparison of the Legionella genus-specific sequences and their read abundances with the ones obtained from all Bacteria. In this universal PCR the variable regions V4 and V5 of the 16S rRNA gene of Bacteria were amplified using the modified universal primers 519 F (5’-CAGCAGCCGCGGTAATAC-3’) and 907R (5’-CCGTCAATTCCTTTGAGTTT-3’) . All reactions were carried out with 50 μl reaction mixtures. The target-specific PCR mix comprised 0.1 mM concentration of each dNTP (Bioline, Luckenwalde, Germany), 0.4 mM MgCl2, 1X PCR reaction buffer, 0.03 U of HotStarTaq Polymerase (Qiagen, Hilden, Germany), 0.4 μM of each primer and 1 ng of environmental DNA. Cycling was performed on a Biorad Thermo cycler 96-well iCycler with an initial denaturation step at 95 °C for 15 min; 30 cycles of 1 min at 95 °C, 40 s at 55 °C, 40 s at 72 °C followed by a final extension at 72 °C for 10 min. The multiplexing PCR was performed as described above for the HotStarTaq DNA polymerase-based amplification and the primers were analogous to the ones provided in Additional file 1: Table S1, i.e. the primer pair Lgsp17F/Lgsp28R was replaced by primers 519 F/907R.
Assessment of sensitivity and quantitative precision of NGS
For the assessment of sensitivity of NGS and its potential for absolute quantification of Legionella species a certified standard of genomic DNA was used. For these spike-in experiments ten-fold dilutions of a molecular standard based on genomic DNA of L. pneumophila ATCC 33152T (GPS-Genetic PCR Solutions, Alicante, Spain) ranging from 105 to 101 genome copies were added to triplicate 5 ng environmental DNA aliquots of a drinking water sample. The sample without spiked L. pneumophila was also analysed. In addition, the same range of L. pneumophila was spiked into nuclease free water to verify the existence of amplification inhibition and interference of bacterial community in the quantification of L. pneumophila in the drinking water sample. All samples were amplified by the Legionella genus-specific assay and after pooling of equal volumes of NGS libraries, sequenced with Illumina MiSeq platform.
For the determination of specificity and precision of the relative abundance assessment, a mock community was constituted by adding equivalent concentrations of 16S rDNA of the following nine strains: L. anisa ATCC 35292T, L. bozemanii ATCC 33217T, L. feelei ATCC 35072T, L. longbeachae ATCC 33462T, L. micdadei ATCC 33218T, L. pneumophila ATCC 33152T, Acinetobacter baumannii DSM 30008, Klebsiella pneumoniae DSM 30104 and Pseudomonas aeruginosa DSM 50071T. Technical replicates (triplicates) of the defined synthetic community were amplified by the described Legionella genus-specific protocol and sequenced by the Illumina MiSeq platform.
Determination of sequence accuracy and quality
The error rate of sequence reads was obtained for NGS assays with HotStarTaq and KAPA HiFi DNA polymerases by using internal controls consisting of nuclease-free water samples spiked in with genomic DNA of L. pneumophila ATCC 33152T (GPS-Genetic PCR Solutions, Alicante, Spain). Technical triplicates were used. 5,000 paired-end merged reads of each triplicate were statistically sampled, for a total of 15,000 paired-end merged reads, and an alignment algorithm was applied against the L. pneumophila reference sequence, allowing detection of nucleotide substitutions, indels and unknown nucleotides, and subsequent calculation of the respective error rates.
Data processing, clustering and taxonomic assignment
The pre-processing of Illumina MiSeq reads, i.e. merging of paired-end reads and barcode-based de-multiplexing, was done using Mothur software package (v.1.34.0) . Reads were quality controlled and aligned with the bioinformatics pipeline of the SILVA project . Each read was aligned using SINA (v1.2.10) . Reads shorter than 200 aligned nucleotides and reads with more than 2% of ambiguities or 2% of homopolymers were excluded for downstream analysis. Reads with a low alignment quality (50 alignment identity, 40 alignment score reported by SINA), suggesting putative contaminations and artefacts, were also removed. Next, filtered reads were dereplicated, clustered and binned/classified. Dereplication (identification of identical reads ignoring overhangs) was done using cd-hit-est  using identity criterion level of 1.00. Clustering was also completed with cd-hit-est using identity criterion of 0.98. The representative sequence of each OTU was then classified using BLAST+ (version 2.2.28+)  against the non-redundant version of the SILVA SSU 115 NR dataset  with standard settings. Reads presenting a best blast hit values below the value 0.93 for the function (% sequence identity + % alignment coverage)/2 were not classified. For species-level classification, additional data processing steps were performed. A conservative approach with removal of singletons, doubletons, tripletons and quadrupletons, was pursued because application of chimera checking on low diversity genus-specific libraries is not very effective . Filtered OTU reference sequences after trimming of primer sequence were blasted using BLASTN+ (version 2.2.28+)  against a database including curated and truncated (only V3-V4 region amplified) 16S rRNA sequences of Legionella species that were retrieved from the non-redundant version of the SILVA SSU 115 NR dataset. OTU reference sequences were assigned to the species level when sequence identity was higher or equal to 97%. This threshold was defined after amplification and sequencing of L. pneumophila ATCC 33152T with HotStarTaq and KAPA HiFi DNA polymerases (Additional file 1: Figure S1). Our results revealed that at a 97% threshold more than 98% of the sequences generated were identified as L. pneumophila. For both acknowledged species and sequence clusters not assigned to a known species, the term phylotype will be used throughout the study.
Quantification of Bacteria and Legionella microbiome by real-time PCR
For quantitative real-time PCR of the Bacteria the same 519 F/907R primers were used as for the NGS approach. For this real-time PCR the same calibration standard and DNA extracts were applied as above. The PCR mixture (20 μl total) contained 10 μl of LC480 SYBR Green Master Mix (2x, Roche Diagnostics, Mannheim, Germany), 1.6 μl pan-bacterial primers (final concentration: 0.8 μM each) and 2 to 4 ng of environmental DNA in 2 μl of PCR grade water. An initial denaturing step at 95 °C for 5 min was followed by 45 cycles of 95 °C for 30 s, 60 °C for 30 s and 72 °C for 30 s.
Legionella species abundance was quantified using 16S rRNA gene primers Lgsp17F/Lgsp28R as for the NGS library preparation. The qPCR reaction consisted of a total volume of 20 μl reaction mixture containing 10 μl of Master Mix as above, 2 μl of each primer (final concentration: 0.5 μM each, Eurofins MWG Operon, Ebersberg, Germany) and 1 to 5 ng of template DNA. Amplification cycling conditions consisted of an initial 5 min cycle at 95 °C followed by 45 cycles of 30 s at 95 °C, 60 °C and 72 °C each. All real-time PCR assays were run on a Light Cycler 480 (Roche Diagnostics, Mannheim, Germany) and potential contaminations and development of primer dimers were determined by melting curve analyses.
The statistical comparison of two datasets was performed using parametric paired t-test, or Welch’s t-test in case of unequal variances. The null-hypothesis was rejected when P-value <0.05. The correlation between two variables was determined by Spearman rank of correlation (rs), or Pearson correlation (r) in case of linear relationship. T-tests and correlations were performed with Past 3 (version 3.09) . Rarefaction curves to evaluate sequencing depth were generated using Analytic Rarefaction 1.3 .
Alpha-diversity metrics and Good’s coverage were calculated with software package Explicet . Multivariate analyses were performed using PRIMER (Version 7.0.7) . Bray-Curtis similarity  matrices were generated, by comparing the standardised, untransformed abundances of phylotypes, and represented by non-metric multidimensional scaling (MDS) plots. Sample cluster analysis was performed with SIMPROF , whereas potentially significant differences between groups of samples were determined using ANOSIM .
Sequence accuracy and quality of generated NGS libraries
In order to minimise error rate, validation of our deep sequencing approach was performed with proofreading DNA polymerase KAPA HiFi and compared with the extensively used HotStarTaq DNA polymerase. With this purpose, we investigated the type, frequency and distribution of errors that occur during the amplification and sequencing steps. After NGS analysis of DNA from the L. pneumophila genome standard in nuclease-free water 6,315,000 bases of sequence information was compared. These sequences had an average error rate of 0.383% ± 0.28% that was significantly lower when compared to HotStarTaq enzyme (0.446% ± 0.283%) (T-test, P < 0.05). Erroneous sequencing was observed in all the 421 bases, with both enzymes, indicating no error-free base position in the sequences amplified. The minimum and maximum error frequencies per base with KAPA HiFi were lower than with HotStarTaq (0.033% and 1.971% against 0.067% and 2.159%) Additionally to a lower minimum, mean and maximum error frequency values per position, 65.1% of the sequences analysed were identified as error-free with KAPA HiFi, while, with HotStarTaq, this number was reduced to only 39.3% (Additional file 1: Figure S1).
When breaking down the error rate by error type, substitutions were the foremost, representing 84.1% of the total errors detected in the libraries generated by KAPA HiFi (Additional file 1: Figure S2a). The 16S rRNA sequences had an average substitution rate per base of 0.322% ± 0.277%, which represents a significant accuracy improvement to libraries amplified by HotStarTaq DNA polymerase (0.383% ± 0.275%) (T-test, P < 0.05). Unknown bases (Ns) were the next most frequent errors and had an average rate per base of 0.032% with KAPA HiFi and 0.035% with HotStarTaq. Deletions and insertions made a smaller contribution to the global error rate. Insertions with KAPA HiFi and HotStarTaq had an average rate per base of 0.011% and 0.008%, respectively. A comparison of the substitution error profiles between enzymes showed different patterns in relation to individual bases (Additional file 1: Figure S2b). With HotStarTaq, transitions and transversions were about equally represented. However, with KAPA HiFi, transversions were clearly more frequent accounting for 67.1% of all substitutions. This imbalance between transitions and transversions is mostly due to a significant reduction of 46.3% in the number of transitions with KAPA HiFi (T-test, P < 0.05).
In silico and in situ evaluation of primer specificity
The primer pair Lgsp17F/Lgsp28R was evaluated in silico using TestPrime tool against SILVA rRNA gene database (Parc SSU release 123) , while allowing none or one mismatch per primer, in order to assess its specificity. The in silico analysis showed a good coverage of 84% of the species in the targeted genus. Considerable amplification of bacteria belonging to related genus Methylococcus (51%) could be expected. Other taxonomic groups, i.e. Steroidobacter, Methylocaldum, Aquicella, Ricketsiella, were covered by the primer studied but at a substantially smaller extent (<25%). The primer pair targets the regions V3 and V4 of the 16S rRNA gene of genus Legionella, which are two of the most hypervariable regions within the genus studied (Additional file 1: Figure S3). The fragment amplified comprising 421 bp allows a very good resolution of most of the Legionella species, including the most clinically relevant species L. pneumophila (Additional file 1: Figure S4). The amplicon presents low intragenomic and intergenomic heterogeneity within species L. pneumophila (Additional file 1: Figure S5) allowing confident identification.
In situ determination of specificity of the primer pair Lgsp17F/Lgsp28R was performed using NGS libraries of triplicates of 7 water samples (two cold drinking water (samples 1 and 2), two hot drinking water [3, 4] and three cooling tower water samples [5–7]). Across these samples, we obtained a total of 1,738,388 sequences after quality filtering and removal of genera with a single read, with the median number of quality sequences in the sample set being 60,455. The mean relative distribution of Legionella sequences per sample was 95.2% ± 10.8%. Additionally, 47 other taxonomic groups, such as genera Methylococcus and Ricketsiella, were amplified but only 4 accounted for 0.1% or more of the sequences classified. These included unclassified Gammaproteobacteria SAR86 clade (4.01% ± 9.53%), genera Methylocaldum (0.29% ± 0.63%), Aquicella (0.10% ± 0.15%) and Massilia (0.10% ± 0.12%). The SAR86 clade was solely identified in the cooling tower water samples while the remaining taxa were detected in more than 80% of the total water samples analysed. In the sample set studied, the highest unspecific amplification was observed in sample 7 and corresponded to 29.14% of the classified sequences, with 25.59% of the sequences affiliated to the SAR86 clade. In all other water samples, the percentage of sequences classified as Legionella surpassed the 95% threshold, with the highest value being 99.85% (sample 4). Overall, the in situ evaluation of the primer pair Lgsp17F/Lgsp28R demonstrated its high specificity for the assessment of Legionella community composition.
Sensitivity and quantification of Legionella species using NGS
Repeatability of NGS approach for determination of Legionella microbiome
Triplicates of the 7 water samples were assessed to determine the technical variation generated by the PCR steps of the Legionella libraries preparation using primer pair Lgsp17F/Lgs28R and their subsequent amplification and sequencing. From this dataset, 1,613,470 16S rRNA Legionella gene sequences were retrieved with a median number of 76,832 quality sequences per sample classified as Legionella. Rarefaction curves for the 21 water samples were plotted (Additional file 1: Figure S6) showing enough sequencing depth as indicated by an average Good’s coverage of 95.3% ± 1.0%.
Alpha-diversity represented by the richness observed and the Shannon’s diversity index (H’) of the individual Legionella microbiomes was calculated based on both OTU and phylotype information for all samples after normalisation (Additional file 1: Table S2). Substantial variation in observed OTU richness (OTUobs) among replicates was detected with coefficient of variation (CV) values reaching values as high as 36.8% (Sample 4), with an overall mean CV value of 20.8%. In contrast, phylotype richness (PTobs) showed a substantially lower variation with an overall CV value of 5.9%. Shannon’s diversity index, based on OTU and phylotype numbers, showed a lower variation among replicates than richness metrics, with an overall mean CV value of 6.5% and 2.9%, respectively. These results seem to be in part due to a significant positive correlation of OTU numbers with sequencing depth (rs = 0.68, P < 0.05), independently of the sample. As OTU-based richness metrics showed lower consistency between replicates, its use is not recommended for Legionella microbiome studies with this approach.
Profiling of the Legionella microbiome in freshwater systems using NGS
Comparison between NGS and real-time PCR quantification of Legionella microbiome within the bacterial microbiome
Detection and quantification of individual Legionella species in freshwater by NGS
A comparison between the Legionella genus-specific NGS approach (using primers Lgsp17F/Lgsp28R) and the pan-bacterial NGS approach (using primers 519 F/907R) was performed to determine the sensitivity and accuracy of the approaches to detect individual Legionella species, using L. pneumophila as an example, in the 7 reference samples. When using the Legionella genus-specific and the pan-bacterial approaches, the average number of sequences per sample (assessed in triplicates) classified as Legionella was 76,832 ± 44,270 and 254 ± 127, respectively.
Accuracy of Legionella genus-specific NGS approach
NGS technologies have emerged as valuable molecular techniques and are becoming more frequently applied in community studies and pathogen discovery and detection . The accuracy of molecular techniques strongly depends on the agreement between the microorganisms present and their sequences derived from the assay and can be strongly influenced by chimeric artefacts, and amplification and sequencing errors [35, 36]. In NGS, the frequency and type of errors is influenced by the enrichment PCR step and the sequencing platform used. Illumina platforms have been reported to account for the lowest number of sequencing errors and to generate mostly substitutions, whereas in other platforms, such as IonTorrent PGM and Roche 454, indels dominate [37, 38]. To maximise genome-coverage and minimise the error-rate during the library preparation step, the use of high-fidelity DNA polymerases with strong proofreading activity, such as KAPA HiFi, has been recommended [39, 40]. The sequence errors can be measured by two parameters, accuracy that is derived from error rate and quality that is derived from the error-free reads. Our results showed that KAPA HiFi outperformed HotStarTaq in both parameters (Fig. 1). Yet, sequences amplified with both enzymes showed overall a low error rate (<0.5%). Additionally, substitutions were the most common errors supporting and being in the range of previous findings with distinct approaches [37, 38, 41]. With KAPA HiFi DNA polymerase, transition substitutions, which are mostly associated with PCR amplification , showed a significant reduction of 46%. This was particularly noticed in substitutions A to G and T to C, which tend to predominate with low fidelity enzymes [43, 44]. Nonetheless, transversions represented a substantial part of the sequencing errors detected which have been reported to be caused mainly by Illumina deep sequencing [41, 45]. These results associated with other findings, such as higher error rate at the end of the reads and in reverse reads confirm that the Illumina sequencing step is responsible for an important portion of the errors observed. Additionally, common error hotspots and systematic errors were found with both enzymes supporting the described existence of DNA sequence and motif-dependent errors [41, 46]. Overall, our results demonstrated the beneficial effect of choosing a high-proofreading enzyme during NGS library preparation.
Specificity and repeatability of the Legionella genus-specific NGS approach
When designing a molecular approach, in silico evaluation comprising the target gene to address and the primer pair to choose is of unequivocal importance. The in silico analysis of primer pair Lgsp17F/Lgsp28R showed not only good coverage and resolution to the species-level of the bacteria englobed in the genus Legionella, but also good specificity. Though in silico predictive analysis is of great help, empirical in situ analysis should be performed as results may differ considerably . The in situ analysis revealed high specificity of the primer pair Lgsp17F/Lgsp28R supporting previous findings . Moreover, the primer pair showed good agreement with in silico data, with most of the expected taxa being amplified. The main exception was amplification of bacteria in clade SAR86, which is constituted of uncultured bacteria mostly described in marine environments . Our results suggest that these bacteria can also be abundant in man-made freshwater environments, especially in cooling towers. In general, the primer pair tested showed high specificity (mostly >95%) and provided adequate coverage and sequencing depth of the targeted Legionella species.
In NGS studies, technical variability is often introduced during sample preparation, library preparation and sequencing . Although resources tend to be more directed to the use of biological replicates with technical variability being occasionally assessed, measurement and proper estimation of technical variability is of critical importance for reliable data interpretation [50, 51]. Our results for the developed NGS approach demonstrated that technical variability as a result of PCR and/or sequencing stochasticity is undisputable in all water samples but it did not significantly modify the Legionella microbiome structure and composition (Fig. 4). Overall, we would draw the same beta-diversity conclusions if no technical replicates were used. These results are in accordance with previous studies using distinct approaches with Illumina technology [52, 53].
Precision and sensitivity of the Legionella genus-specific NGS approach
In PCR-based studies the abundance of microbial taxa is often misrepresented and under-detected due to poor primer coverage and reduced amplification efficiency, inadequate sequencing depth, the presence of spurious amplification and sequencing artefacts, plus incorrect sequence clustering and classification [36, 54]. Our results of the mock community revealed that the approach does not seem to strongly favour the amplification of a particular Legionella species (Fig. 3). Yet, our findings give us hints that Legionella diversity might be slightly overestimated and that taxonomic assignment might be to a small degree wrongly inferred due to possibly incorrect sequences and/or wrong clustering/classification. The spike-in experiments suggest not only high sensitivity of the approach for detection of L. pneumophila but also that sequence abundance is a good indicator of the concentration of the targeted bacteria, allowing precise absolute quantification of L. pneumophila abundance in freshwater environments (Fig. 2). The relative abundance of the Legionella microbiome among the whole bacterial microbiome was determined independently with two real-time measurements using the same primer pairs as for the NGS methods (Fig. 6). This comparison showed good correspondence and correlation between both independent abundance determinations and demonstrated the reliability of the quantification of the genus-specific NGS approach.
In addition, when comparing the developed Legionella genus-specific approach with a broader pan-bacterial approach for the detection of L. pneumophila in freshwater samples, a higher detection and more precise quantitative measurement of L. pneumophila with the Legionella genus-specific approach was shown (Fig. 7 and Additional file 1: Figure S9). A study by Vierheilig et al. , using 454 pyrosequencing, failed to detect faecal indicator E. coli in wastewater, though it was identified by cultivation. These concerns over the capability of a pan-bacterial approach to detect and exactly quantify low-abundant bacteria seem to be confirmed by our results, despite applying a substantially higher sequencing effort with Illumina MiSeq platform. These findings demonstrated that the developed genus-specific assay is a more adequate tool for sensitive detection and quantification of the targeted Legionella species.
Characterisation of Legionella microbiome to the species level by NGS
We applied the Legionella genus-specific NGS approach to characterise the Legionella microbiome diversity and composition of seven representative freshwater samples. Our results revealed that a major percentage of the 16S rRNA sequences analysed (about 71%) were not affiliated to any described species, demonstrating the importance of molecular tools such as NGS to more accurately measure true Legionella diversity and abundances. The analysis of the composition of 7 Legionella communities revealed the presence of 75 phylotypes, demonstrating the great Legionella richness in freshwater environments, which should be more thoroughly studied (Fig. 5). The breakdown of Legionella community composition showed that L. pneumophila was the most abundant species and the most abundant potentially pathogenic species in the community, especially in hot water. These results corroborate the single-strand conformation polymorphism (SSCP) fingerprint analyses of cold and hot drinking water by Lesnik et al.  but revealed a substantially lower taxonomic resolution of this community profiling technique. We obtained 75 Legionella phylotypes whereas the SSCP fingerprints using the same primers only found half the number of phylotypes in almost five times more samples. These findings underline the public health relevance of L. pneumophila and justify its attentive monitoring in man-made freshwater systems. In addition, the data presented are indicative of a considerable spatial and temporal variation of the Legionella microbiome, highlighting the need of a broader and more detailed study of Legionella species at distinct sites and throughout time in order to better comprehend their dynamics, distribution and the factors that might trigger their growth, persistence and dominance.
In summary, the developed Legionella genus-specific Illumina-based approach provides an accurate and reliable detection, quantification and unambiguous taxonomic classification of the Legionella microbiome comprising all major Legionella species. This assay represents an improvement, in terms of sensitivity and precision, to the widely used pan-bacterial approaches and has the potential to be applied to water quality monitoring of Legionella species, most importantly L. pneumophila. In addition, it will be a valuable molecular tool in water research to augment the scientific understanding of the Legionella species dynamics and the biotic and abiotic drivers impacting Legionella community structure and diversity. In general, the developed genus-specific NGS approach should be applicable to other environmentally or medically relevant genera and opens up a new window to look in greater detail into specific parts of the microbiome, especially the patho-biome.
We thank Verena Maiberg, Josefin Koch, Rene Lesnik, Axel Mohr, Monika Pape and Marina Pecellin for their skilful technical support with respect to field work, DNA extraction and real-time PCR analyses.
This work was supported by funds from the Deutsche Forschungsgemeinschaft (DFG-project No. HO-930/5-1/2) to MGH and the EU project AQUAVALENS (No. 311846) jointly supporting JP and MGH. The DFG project provided the design of the study including sample collection and environmental background data. The EU project enabled the development and application of the NGS approach and the data analysis.
Availability of data and materials
All data generated or analysed during this study are included in this published article and its supplementary information files.
RP designed and performed the experiments, did the data analyses and wrote most the manuscript. JP performed the bioinformatics of the raw sequence data. IB provided the real-time analyses and detailed data analyses. MH designed the study and wrote parts of the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Uyaguari-Diaz MI, Chan M, Chaban BL, Croxen MA, Finke JF, Hill JE, et al. A comprehensive method for amplicon-based and metagenomic characterization of viruses, bacteria, and eukaryotes in freshwater samples. Microbiome. 2016;4:20.View ArticlePubMedPubMed CentralGoogle Scholar
- Fields BS, Benson RF, Besser RE. Legionella and Legionnaires’ disease: 25 years of investigation. Clin Microbiol Rev. 2002;15(3):506–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Muder RR, Yu VL. Infection due to Legionella species other than L. pneumophila. Clin Infect Dis. 2002;35(8):990–8.View ArticlePubMedGoogle Scholar
- Beauté J, Zucs P, de Jong B. European Legionnaires’ Disease Surveillance Network. Legionnaires disease in Europe, 2009–2010. Euro Surveill. 2013;18(10):20417.PubMedGoogle Scholar
- van Heijnsbergen E, Schalk JAC, Euser SM, Brandsema PS, den Boer JW, de Roda Husman AM. Confirmed and potential sources of legionella reviewed. Environ Sci Technol. 2015;49(8):4797–815.View ArticlePubMedGoogle Scholar
- Whiley H, Bentham R. Legionella longbeachae and legionellosis. Emerg Infect Dis. 2011;17(4):579–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Buse HY, Schoen ME, Ashbolt NJ. Legionellae in engineered systems and use of quantitative microbial risk assessment to predict exposure. Water Res. 2012;46(4):921–33.View ArticlePubMedGoogle Scholar
- Ditommaso S, Gentile M, Giacomuzzi M, Zotti CM. Recovery of Legionella species from water samples using an internal method based on ISO 11731: suggestions for revision and implementation. Diagn Microbiol Infect Dis. 2011;70(2):200–6.View ArticlePubMedGoogle Scholar
- Whiley H, Taylor M. Legionella detection by culture and qPCR: Comparing apples and oranges. Crit Rev Microbiol. 2014;0(0):1–10.View ArticleGoogle Scholar
- Joly P, Falconnet P-A, André J, Weill N, Reyrolle M, Vandenesch F, et al. Quantitative real-time Legionella PCR for environmental water samples: data interpretation. Appl Environ Microbiol. 2006;72(4):2801–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Collins S, Jorgensen F, Willis C, Walker J. Real-time PCR to supplement gold-standard culture-based detection of Legionella in environmental samples. J Appl Microbiol. 2015;119(4):1158–69.View ArticlePubMedGoogle Scholar
- Tan B, Ng C, Nshimyimana JP, Loh LL, Gin KY-H, Thompson JR. Next-generation sequencing (NGS) for assessment of microbial water quality: current progress, challenges, and future opportunities. Front Microbiol. 2015;6:1027.PubMedPubMed CentralGoogle Scholar
- Eichler S, Christen R, Höltje C, Westphal P, Bötel J, Brettar I, et al. Composition and dynamics of bacterial communities of a drinking water supply system as assessed by RNA- and DNA-based 16S rRNA gene fingerprinting. Appl Environ Microbiol. 2006;72(3):1858–72.View ArticlePubMedPubMed CentralGoogle Scholar
- Henne K, Kahlisch L, Höfle MG, Brettar I. Seasonal dynamics of bacterial community structure and composition in cold and hot drinking water derived from surface water reservoirs. Water Res. 2013;47(15):5614–30.View ArticlePubMedGoogle Scholar
- Hamady M, Walker JJ, Harris JK, Gold NJ, Knight R. Error-correcting barcoded primers for pyrosequencing hundreds of samples in multiplex. Nat Methods. 2008;5(3):235–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Bartram AK, Lynch MDJ, Stearns JC, Moreno-Hagelsieb G, Neufeld JD. Generation of multimillion-sequence 16S rRNA gene libraries from complex microbial communities by assembling paired-end illumina reads. Appl Environ Microbiol. 2011;77(11):3846–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Camarinha-Silva A, Jáuregui R, Pieper DH, Wos-Oxley ML. The temporal dynamics of bacterial communities across human anterior nares. Environ Microbiol Rep. 2012;4(1):126–32.View ArticlePubMedGoogle Scholar
- Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl Environ Microbiol. 2013;79(17):5112–20.View ArticlePubMedPubMed CentralGoogle Scholar
- Kahlisch L, Henne K, Draheim J, Brettar I, Höfle MG. High-resolution in situ genotyping of Legionella pneumophila populations in drinking water by multiple-locus variable-number tandem-repeat analysis using environmental DNA. Appl Environ Microbiol. 2010;76(18):6186–95.View ArticlePubMedPubMed CentralGoogle Scholar
- Schwieger F, Tebbe CC. A new approach to utilize PCR–single-strand-conformation polymorphism for 16S rRNA gene-based microbial community analysis. Appl Environ Microbiol. 1998;64(12):4870–6.PubMedPubMed CentralGoogle Scholar
- Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75(23):7537–41.View ArticlePubMedPubMed CentralGoogle Scholar
- Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 2007;35(21):7188–96.View ArticlePubMedPubMed CentralGoogle Scholar
- Pruesse E, Peplies J, Glöckner FO. SINA: Accurate high-throughput multiple sequence alignment of ribosomal RNA genes. Bioinformatics. 2012;28(14):1823–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.View ArticlePubMedGoogle Scholar
- Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10(1):421.View ArticlePubMedPubMed CentralGoogle Scholar
- Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41(D1):D590–6.View ArticlePubMedGoogle Scholar
- Haas BJ, Gevers D, Earl AM, Feldgarden M, Ward DV, Giannoukos G, et al. Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons. Genome Res. 2011;21(3):494–504.View ArticlePubMedPubMed CentralGoogle Scholar
- PAST: Paleontological Statistics Software Package for Education and Data Analysis. Available from: http://palaeo-electronica.org/2001_1/past/issue1_01.htm.
- UGA Stratigraphy Lab. Available from: http://strata.uga.edu/software/anRareReadme.html.
- Robertson CE, Harris JK, Wagner BD, Granger D, Browne K, Tatem B, et al. Explicet: graphical user interface software for metadata-driven management, analysis and visualization of microbiome data. Bioinformatics. 2013;29(23):3100–1.View ArticlePubMedPubMed CentralGoogle Scholar
- Clarke K, Warwick R. Change in marine communities: an approach to statistical analysis and interpretation. 3rd ed. Plymouth, UK: Primer-E Ltd; 2014.Google Scholar
- Bray JR, Curtis JT. An ordination of the upland forest communities of southern Wisconsin. Ecol Monogr. 1957;27(4):326–49.View ArticleGoogle Scholar
- Clarke KR. Non-parametric multivariate analyses of changes in community structure. Aust J Ecol. 1993;18(1):117–43.View ArticleGoogle Scholar
- Pallen MJ. Diagnostic metagenomics: potential applications to bacterial, viral and parasitic infections. Parasitology. 2014;141(14):1856–62.View ArticlePubMedPubMed CentralGoogle Scholar
- Kircher M, Heyn P, Kelso J. Addressing challenges in the production and analysis of illumina sequencing data. BMC Genomics. 2011;12:382.View ArticlePubMedPubMed CentralGoogle Scholar
- Pinto AJ, Raskin L. PCR biases distort bacterial and archaeal community structure in pyrosequencing datasets. PLoS One. 2012;7(8):e43093.View ArticlePubMedPubMed CentralGoogle Scholar
- Loman NJ, Misra RV, Dallman TJ, Constantinidou C, Gharbia SE, Wain J, et al. Performance comparison of benchtop high-throughput sequencing platforms. Nat Biotech. 2012;30(5):434–9.View ArticleGoogle Scholar
- Jünemann S, Sedlazeck FJ, Prior K, Albersmeier A, John U, Kalinowski J, et al. Updating benchtop sequencing performance comparison. Nat Biotech. 2013;31(4):294–6.View ArticleGoogle Scholar
- Quail MA, Otto TD, Gu Y, Harris SR, Skelly TF, McQuillan JA, et al. Optimal enzymes for amplifying sequencing libraries. Nat Methods. 2012;9(1):10–1.View ArticleGoogle Scholar
- Oyola SO, Otto TD, Gu Y, Maslen G, Manske M, Campino S, et al. Optimizing Illumina next-generation sequencing library preparation for extremely AT-biased genomes. BMC Genomics. 2012;13(1):1.View ArticlePubMedPubMed CentralGoogle Scholar
- Schirmer M, Ijaz UZ, D’Amore R, Hall N, Sloan WT, Quince C. Insight into biases and sequencing errors for amplicon sequencing with the Illumina MiSeq platform. Nucleic Acids Res. 2015;43(6):e37.View ArticlePubMedPubMed CentralGoogle Scholar
- Brodin J, Mild M, Hedskog C, Sherwood E, Leitner T, Andersson B, et al. PCR-induced transitions are the major source of error in cleaned ultra-deep pyrosequencing data. PLoS One. 2013;8(7):e70388.View ArticlePubMedPubMed CentralGoogle Scholar
- McInerney P, Adams P, Hadi MZ, McInerney P, Adams P, Hadi MZ. Error rate comparison during polymerase chain reaction by DNA polymerase. Mol Biol Int. 2014;2014:e287430.View ArticleGoogle Scholar
- Bracho MA, Moya A, Barrio E. Contribution of Taq polymerase-induced errors to the estimation of RNA virus diversity. J Gen Virol. 1998;79(12):2921–8.View ArticlePubMedGoogle Scholar
- Dohm JC, Lottaz C, Borodina T, Himmelbauer H. Substantial biases in ultra-short read data sets from high-throughput DNA sequencing. Nucleic Acids Res. 2008;36(16):e105.View ArticlePubMedPubMed CentralGoogle Scholar
- Nakamura K, Oshima T, Morimoto T, Ikeda S, Yoshikawa H, Shiwa Y, et al. Sequence-specific error profile of Illumina sequencers. Nucleic Acids Res. 2011;39(13):e90.View ArticlePubMedPubMed CentralGoogle Scholar
- Morales SE, Holben WE. Empirical testing of 16S rRNA gene PCR primer pairs reveals variance in target specificity and efficacy not suggested by in silico analysis. Appl Environ Microbiol. 2009;75(9):2677–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Lesnik R, Brettar I, Höfle MG. Legionella species diversity and dynamics from surface reservoir to tap water: from cold adaptation to thermophily. ISME J. 2016;5:1064–80.View ArticleGoogle Scholar
- Dupont CL, Rusch DB, Yooseph S, Lombardo M-J, Alexander Richter R, Valas R, et al. Genomic insights to SAR86, an abundant and uncultivated marine bacterial lineage. ISME J. 2012;6(6):1186–99.View ArticlePubMedGoogle Scholar
- Robasky K, Lewis NE, Church GM. The role of replicates for error mitigation in next-generation sequencing. Nat Rev Genet. 2014;15(1):56–62.View ArticlePubMedGoogle Scholar
- Ficetola GF, Pansu J, Bonin A, Coissac E, Giguet-Covex C, De Barba M, et al. Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data. Mol Ecol Resour. 2015;15(3):543–56.View ArticlePubMedGoogle Scholar
- Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Huntley J, Fierer N, et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 2012;6(8):1621–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Kennedy K, Hall MW, Lynch MDJ, Moreno-Hagelsieb G, Neufeld JD. Evaluating bias of Illumina-based bacterial 16S rRNA gene profiles. Appl Environ Microbiol. 2014;80(18):5717–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Lee CK, Herbold CW, Polson SW, Wommack KE, Williamson SJ, McDonald IR, et al. Groundtruthing next-gen sequencing for microbial ecology-biases and errors in community structure estimates from PCR amplicon pyrosequencing. PLoS One. 2012;7(9):e44224.View ArticlePubMedPubMed CentralGoogle Scholar
- Vierheilig J, Savio D, Ley RE, Mach RL, Farnleitner AH, Reischer GH. Potential applications of next generation DNA sequencing of 16S rRNA gene amplicons in microbial water quality monitoring. Water Sci Technol. 2015;72(11):1962–72.View ArticlePubMedPubMed CentralGoogle Scholar