Deconstructing the Polymerase Chain Reaction II: an improved workflow and effects on artifact formation and primer degeneracy

DNA templates

A single microbial genomic DNA (gDNA) sample obtained from chinchilla feces was used for this study. Chinchilla feces were acquired from the cage of pet animals and the sample serves as a representative complex microbial community. No specific permissions were required for collection of the chinchilla feces. The fecal sample was extracted using the PowerSoil DNA extraction kit (Mo Bio Laboratories, Carlsbad, CA, USA).

Primer synthesis

The primers used for this study are 341F (CCTACGGGAGGCAGCAG) (Muyzer, De Waal & Uitterlinden, 1993; Caporaso et al., 2011) and 806R (GGACTACHVGGGTWTCTAAT) (Caporaso et al., 2011; Walters et al., 2016). The 806R primer pool is 18-fold degenerate, with theoretical melting temperatures ranging from 54.7 °C to 61 °C. Melting temperatures of the primers were calculated using the OligoAnalyzer3.1 tool (Owczarzy et al., 2008), assuming 250 nM primer concentration, 2 mM Mg²⁺, and 0.2 mM dNTPs. Synthesis of the primers was performed either as single degenerate primer pools (standard approach), or as individual primers without degeneracies by Integrated DNA Technologies (IDT; Coralville, IA). Primers were synthesized as LabReady and ordered at a fixed concentration of 100 micromolar. Primers contained common sequence linkers (CS1 and CS2) at the 5′ ends, as shown in Table 1. Linker sequences are required for the later incorporation of Illumina sequencing adapters and sample-specific barcodes.

Standard PCR protocol

The standard PCR protocol or targeted amplicon sequencing (TAS) protocol is a two-stage NGS library preparation protocol for generating barcoded amplicons ready for Illumina sequencing, and was performed as described previously (Naqib et al., 2018) (Fig. 1A). Briefly, gDNA was amplified by PCR with primers CS1_341F and CS2_806R. The first stage PCR reaction was conducted in a total reaction volume of 10 µl. Each reaction contained 5 µl of MyTaq HS master mix (Bioline, Taunton, MA, USA), 0.5 µl of each primer or degenerate primer at a concentration of 5 µM (e.g., CS1_341F and CS2_806R; leading to a 250 nM working concentration), 10 ng of gDNA template, and water up to 10 µl total volume. The first stage of the PCR was conducted using the following thermocycling conditions: 95 °C for 5 min, followed by 28 cycles of 95 °C for 30 s, annealing temperature (from 40 °C to 60 °C) for 30 s, 72 °C for 30 s; and a final elongation step at 72 °C for 7 min. Subsequently, a second PCR amplification was performed in 10 µl reactions in 96-well plates to incorporate Illumina sequencing adapters and a sample-specific barcode. A mastermix for the entire plate was made using the MyTaq HS 2X mastermix. Each well received a separate primer pair with a unique 10-base barcode, obtained from the Access Array Barcode Library for Illumina (Item: 100-4876; Fluidigm, South San Francisco, CA, USA). These Access Array primers contained the CS1 and CS2 linkers at the 3′ ends of the oligonucleotides. One µl of reaction mixture from the first stage amplification was used as input template for the second stage reaction, without cleanup. Cycling conditions were as follows: 95 °C for 5 min, followed by eight cycles of 95 °C for 30″, 60 °C for 30″and 72 °C for 30″. A final, 7-min elongation step was performed at 72 °C. Samples were pooled and sequenced on an Illumina MiSeq employing V2 chemistry and 2x250 base reads.

Deconstructed PCR (DePCR) Protocol

As with the TAS method, the DePCR method is also a two-stage PCR process (Fig. 1C) and is a modification of the previously described PEX PCR method (Fig. 1B). For each sample, the first stage reaction of DePCR (four total cycles) was conducted in a 96-well plate with each well containing 5 µl of MyTaq master mix, 0.5 µl of each primer or degenerate primer at a concentration of 5 µM (e.g., CS1_341F and CS2_806R; leading to a 250 nM working concentration), 10 ng of template, 1 µl Access Array Barcode Library containing a unique sample-specific barcode, and water up to 10 µl. The thermocycler conditions for first stage were composed of two cycles of denaturation at 95 °C for 5 min and annealing (40 °C–60 °C, depending on experiment) for 20 min, followed by two cycles of denaturation for 5 min at 95 °C and annealing at 60 °C for 20 min, and a final extension temperature of 72 °C for 10 min. For temperature gradient experiments, annealing temperatures of 40 °C, 45 °C, 50 °C, 55 °C, and 60 °C were tested. For single reverse primer variant (RPV) analyses, an annealing temperature of 50 °C was used for both TAS and DePCR amplification reactions. Subsequently, a pool composed of 5 µl from the first reaction of each sample was collected and processed for cleanup using AMPure XP beads (Beckman-Coulter) at 0.7X per the manufacturer’s recommendations. The cleaning step was performed twice, sequentially. A final elution volume of 20 µl was used to concentrate the sample prior to the second stage of the DePCR reaction. The second stage reactions were conducted in a final volume of 20 µl; the reaction contained 10 µl of MyTaq HS master mix, 1 µl of Illumina P5 (AATGATACGGCGACCACCGA) and P7 (CAAGCAGAAGACGGCATACGA) primers, 2 µl of purified template from pooled first stage PCR, and water up to 20 µl. The thermocycler conditions were: 95 °C for 3 min, 30 cycles at 95 °C for 15 s, 60 °C for 30 s and 72 °C for 30 s. Prior to sequencing the pooled libraries were purified using a Pippin Prep DNA Size Selection System (Sage Science), employing a 2% agarose gel cassette and selecting for fragment sizes from 450-600 bp. Sequencing of the amplified pool was performed on an Illumina MiSeq employing V2 chemistry and 2x250 base reads, and demultiplexing of sequence data was performed on instrument. Library preparation and sequencing were performed at the UIC Sequencing Core (UICSQC).

Sequence data analysis

Raw sequence FASTQ files were merged using the software package PEAR  (Zhang et al., 2013), with default parameters. For analysis of primer utilization profiles, merged sequences were trimmed using the software package trimmomatic (Edgar, 2010), and sequences shorter than 400 bases and longer than 500 bases were removed. Using Unix bash scripting, exact primer sequences were searched for within these trimmed sequences and counted. For microbial community analysis, PEAR-merged sequences were initially processed through the software package CLC genomics (v10; Qiagen, Aarhus, Denmark) to remove primer sequences, to perform quality trimming (below Q20 removed), and size trimming (below 400 bases removed). Sequences were then screened for chimeras using the USEARCH61 algorithm  (Edgar, 2010), and putative chimeric sequences were removed from the dataset. Subsequently, sequences were pooled and clustered into operational taxonomic units (OTUs) at a threshold of 97% similarity (QIIME v1.8.0)  (Caporaso et al., 2010). Representative sequences from all OTUs were annotated using the UCLUST algorithm and the Greengenes 13_8 reference database (McDonald et al., 2012b), and a biological observational matrix (BIOM) was generated by this annotation pipeline (McDonald et al., 2012a). The BIOM file was analyzed and visualized using the software package Primer7 (Clarke & Gorley, 2015) and the R environment (R Core Team, 2013). The R package ‘vegan’ (Oksanen et al., 2011) was employed to generate alpha diversity indices (Shannon, richness, and evenness indices) and to produce rarefied BIOM files. Bray–Curtis dissimilarity indices were calculated within the R package ‘vegan’ and these indices were used to evaluate differences in composition between samples. Analysis of similarity (ANOSIM) calculations were performed at the taxonomic level of genus, using square root transformed data. Initial analysis and processing of the samples was performed using QIIME (v1.8.0) package scripts. Metric multi-dimensional scaling (mMDS) plots were generated using the cmdscale and ggplot2 functions (Wickham, 2016) within the R programming environment. Ellipses, representing a 95% confidence interval around group centroids, were drawn assuming a multivariate t-distribution. Some visualizations were performed using the software package OriginPro 2018 (OriginLab, Northampton, MA, USA). Rarefaction and group-significance testing (i.e., non-parametric Kruskal-Wallis test) were performed using the QIIME software package.

Data archive

Raw sequence data files were submitted in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI). The BioProject identifier of the samples is PRJNA506229. Full metadata for each sample are provided in Table S1.

Theory

The Deconstructed PCR (DePCR) method is based on the polymerase-exonuclease (PEX) PCR method described previously (Green, Venkatramanan & Naqib, 2015). We previously noted that the first two cycles of PCR are unique in that no amplification of the template is performed. Rather, linear copying of the template nucleic acid prepares the reaction for exponential amplification, starting in the third cycle. In the prior manuscript, linear copying of the original gDNA template was separated from exponential amplification of target copies using exonuclease I (Fig. 1B). Locus-specific primers containing 5′ linker sequences anneal to genomic DNA during two cycles of amplification. Subsequently, exonuclease I was used to remove unused primers from reaction mixtures. Finally, the copied templates were exponentially amplified using primers targeting the 5′ linker sequences but not the source genomic DNA. This approach is viable, but cumbersome due to the need for exonuclease treatment of each sample, and for individual amplification of each sample with primers containing Illumina sequencing adapters and sample-specific barcodes.

We modified the original PEX PCR protocol by including both locus-specific primers containing 5′ linkers as well as primers with Illumina sequencing adapters, sample-specific barcodes, and 3′ linkers together in the first linear stage of the DePCR reaction (Figs. 1C and 2). Thus, the DePCR approach combines primer sets used in both stage A and B of the PEX PCR method in the same reaction. In addition, four cycles of linear copying are performed, instead of two as in the PEX PCR method (Figs. 1 and 2). The resulting products are target copies containing Illumina sequencing adapter sequences, sample specific barcodes, linker sequences, and the region of interest. The four cycles of copying serve to prepare the templates for exponential amplification but also (unlike PEX PCR) incorporate a sample-specific barcode so that samples can be pooled and amplified exponentially simultaneously in the second stage. As with PEX PCR, the linear amplification stage—if operating at 100% efficiency—does not increase the total number of targets from that present in the source template DNA.

After linear copying during the first four cycles, the reactions are pooled and purified to remove unincorporated primers. It is essential for the proper functioning of the method that the primers from the initial stage of the reaction are completely removed; otherwise these locus-specific primers continue to interact with template and amplicons during exponential amplification cycles. We observed that a single cleanup using AMPure XP beads (0.7X) was not sufficient to fully remove all primers; therefore, a double cleanup (i.e., two sequential AMPure XP 0.7X cleanups of the pooled reactions) is performed. The final purified DNA includes a range of DNA types, but only the fragments that contain Illumina sequencing adapters at both ends of the molecule have been generated only through linear copying steps and are available for amplification using Illumina P5 and P7 primers (Fig. 2). The entire pool is then used as input template for subsequent amplification using primers consisting of Illumina P5 and P7 sequences. Linear-copied DNA fragments from all samples within the pool, each now containing a sample-specific barcode, are thus subject to exponential amplification simultaneously. One useful feature of this approach is that hundreds of samples can be amplified simultaneously within a single reaction. The theoretical advantages of this novel workflow include: (1) the elimination of a separate exonuclease step for each sample, (2) the rapid reduction of many reactions into a single reaction for purification and exponential amplification, and (3) all associated benefits of the prior PEX PCR, in which linear and exponential amplification stages of PCR are isolated from each other and where locus-specific primers are only active for two linear cycles of copying.

Validation of the DePCR method

To assess the effects of amplification method (TAS vs DePCR) and annealing temperature on observed microbial community structure, a single genomic DNA sample was amplified across multiple annealing temperatures using both amplification strategies. Five technical replicates for each condition were performed, and amplicons were sequenced together. The data were analyzed to determine if there were significant differences in sequence metrics (chimera formation), alpha diversity (richness and Shannon index), and observed community structure (beta diversity analyses including multi-dimensional scaling, analysis of similarity (ANOSIM), and taxon-level group-significance testing). Rates of detectable chimera formation were several orders of magnitude lower with the DePCR pipeline relative to the TAS pipeline, regardless of annealing temperature (Table 2). Average chimera detection rate for TAS-processed samples range from 5.16 to 6.53%, while that for DePCR-processed samples ranged from 0.03–0.1%; this difference was significant at all annealing temperatures tested (ANOVA, P < 0.001). Low rates of detectable chimeras were found in all experiments conducted with DePCR, with averages in the range of 0.01–0.1% (Table 2). Conversely, alpha diversity metrics (genus-level richness and Shannon index), were slightly and significantly higher in TAS-based analyses relative to DePCR. Genus-level richness was on average from 1.06–1.21X higher in TAS analyses relative to DePCR, across annealing temperatures from 40 °C to 60 °C (one-way ANOVA; p values ranged from 1.9E-5 to 1.3E-1; Table 3). Shannon indices were from 1.03–1.06X higher in TAS analyses relative to DePCR across annealing temperatures from 40 °C to 60 °C (ANOVA; p<8.13E-4; Table 3).

A strong, significant effect of annealing temperature on the observed microbial community structure was seen in both TAS and DePCR amplification methods (Fig. 3A). Although the overall scale of difference between TAS and DePCR was modest (maximum Bray–Curtis dissimilarity between samples = 0.23 between a TAS sample with 60 °C annealing temperature and a DePCR sample with 40 °C), there was a significant effect of amplification method on observed microbial community at all temperatures (Table S2). Two-way ANOSIM analyses indicated significant differences by temperature across methods (R = 0.832; p = 0.0001; Fig. 3B), and by amplification method across temperatures (R = 0.988; p = 0.0001; Fig. 3C). Similar trends were observed for increases in annealing temperature in both methods, with temperature loading primarily on MDS axis 1. As previously noted (Green, Venkatramanan & Naqib, 2015), greater variability in observed microbial community structure was noted with DePCR with low annealing temperature, particularly at 40 °C (Fig. 3A).

One key feature of the DePCR methodology is the ability to determine which primers in a degenerate pool are interacting with the source genomic DNA. This is achieved as the exponential amplification of the template is performed using primers targeting Illumina sequencing adapters and not the locus-specific primers (Figs. 1 and 2). Locus-specific primers only interact with the gDNA and the first linear copies of gDNA during the first two cycles of the DePCR method. These primer sequences are retained during exponential amplification with primers targeting linker sequences. Conversely, in standard PCR, the locus-specific primers interact with both the genomic DNA template and with copies made from the genomic DNA during exponential amplification; thus, information regarding primer-gDNA template interactions are lost  (Green, Venkatramanan & Naqib, 2015). We thus examined the so-called “primer utilization profiles” (PUPs) for these reactions (Fig. 4). The relative frequency of each of the 18 unique primer variants is shown for each replicate at each PCR condition (temperature x method). Standard PCR amplification protocol (TAS) removes primer-template interaction information as primer-amplicon interactions throughout the amplification reaction tolerate mismatches; all 18 primer variants are used at similar frequencies, regardless of annealing temperature (Fig. 4A). Some patterning is observed in the TAS method, but overall diversity of primer utilization is extremely high and only small differences were observed between temperatures of 40−60 °C (Fig. 4B). The average Shannon index for PUP profiles of TAS samples across all annealing temperatures was 2.859–2.864; the maximum possible natural log Shannon index for 18 features is 2.890. This PUP diversity profiling demonstrates that for standard TAS PCR, the primers used in copying throughout the amplification reaction are not dependent on annealing temperature.

Conversely, a strong effect of annealing temperature is observed on the PUP of samples amplified using the DePCR protocol (Figs. 4A and 4B). A shift in PUP patterning is observed with increasing annealing temperature, and at 60 °C two primer variants (RPV5 and RPV15) dominate. At lower annealing temperatures, a broader range of primers are utilized in the initial stages of gDNA copying. The relationship between annealing temperature and primer utilization richness (here represented as the Shannon index) was best fit with a polynomial equation and is shown in Fig. 4C. As annealing temperature increases, fewer and fewer primer variants interact with the gDNA template. Conversely, at the lowest tested annealing temperature of 40 °C, the Shannon index of the DePCR amplicons nearly matched that of the TAS. Several primer variants, however, including RPVs 10, 12, 14 and 18, were poorly utilized in DePCR amplifications regardless of annealing temperature (Fig. 4A). These four variants included variants with high melting temperatures (57.4, 57.5, 58 and 59.8 °C), while the two most utilized RPVs at PCR annealing temperatures of 60 °C had moderate to high annealing temperatures (56.4 and 58.7 °C). Thus, the melting temperature of the primer did not directly correlate with utilization at different PCR annealing temperatures in this system. The observed primer utilization profiles represent a template-specific phenomenon, and different PUPs would be recovered with different DNA templates.

Determination of linearity in DePCR amplification

In the DePCR protocol, after four initial cycles of linear copying during the first stage of DePCR, samples are pooled prior to purification and second stage amplification with Illumina P5 and P7 primers. The pooling of samples can only be performed because of the incorporation of a sample-specific unique barcode for each sample during the first stage. During the second stage amplification, primers targeting the Illumina adapters are used for amplification, and all templates from all samples are amplified simultaneously (Fig. 1C). Since there is no opportunity for primer-template bias during the second stage (i.e., Stage B of Fig. 1C) of amplification (all amplifiable template molecules contain Illumina sequencing adapters) and primers are non-degenerate, the relative abundance of template molecules from a single sample within the pool should be maintained during amplification. To determine if the relative abundance of template molecules from each sample was maintained in the DePCR protocol, we performed an experiment in which input gDNA (feces) was varied from 1.25 ng to 20 ng per 10 µl reaction. All input levels were performed with five technical replicates. After the first stage (4 cycles) of the DePCR, all replicates from all gDNA input levels were pooled in equal volume and purified. The purified product was then amplified with P5 and P7 primers, and the final pool sequenced. We first assessed whether the input DNA concentration was correlated with the total number of reads generated using this approach (Fig. S1). Since all samples were amplified together, and low input DNA samples should theoretically provide fewer molecules to the combined pool, we hypothesized that a linear relationship should exist between input DNA in the first stage and the number of reads generated per sample. A significant positive correlation between input gDNA concentration and absolute number of reads recovered from each sample was observed, though substantial variability at each input concentration was observed (R² = 0.58, Fig. S1C). We also sought to determine if the input gDNA concentration from the same sample had a significant effect on the observed microbial community structure. Although there was a positive correlation between input gDNA and total number of sequences recovered, we observed no significant effect of input gDNA on the microbial community structure (Fig. S1A; Global ANOSIM R =  − 0.034; p = 0.79). Similarly, no significant difference in primer utilization was observed with different gDNA input concentrations (Fig. S1B). Thus, increasing input gDNA concentration alters the number of molecules passing to the second stage of the DePCR reaction, but within the measured concentration range did not affect the primer utilization profile or final observed microbial community structure.

Assessing the effect of individual primers in a degenerate primer pool

Degenerate primer pools are generally used to amplify genomic DNA, although not all primers actively interact with the source gDNA (Fig. 4A). This degenerate mixture of primers is employed to target a broad range of taxa, and the presence of additional primer variants in pools has been shown to improve detection of known microbial lineages (Hayashi, Sakamoto & Benno, 2004; Frank et al., 2008; Parada, Needham & Fuhrman, 2016; Apprill et al., 2015). In standard PCR, all primers do eventually interact with amplified copies of gDNA during the many cycles of exponential amplification; however, many primers do not interact with the source genomic DNA due to preferential annealing of other primers (Fig. 4A). We sought, therefore, to determine how much microbial diversity could be detected using each primer variant independently in PCR reactions using both the TAS and DePCR methods. In addition, we sought to determine how the observed microbial community structure differed by single primer variant usage. We hypothesized that the single primer variant PCR would better approximate degenerate primer pools when using the DePCR method relative to the TAS method, as our prior work showed that a deconstructed PCR approach was more tolerant of mismatches between primer and gDNA template than TAS (Green, Venkatramanan & Naqib, 2015). The tolerance of mismatches may lead to better capture of microbial community diversity when a greater number of mismatches between primer and template are present, as is expected in a single primer PCR. To explore this, we PCR-amplified a single gDNA template (feces) with the 18 unique reverse primer variants (RPVs) from the degenerate primer pool. Each reaction was performed in technical duplicates, and each reaction was performed using the DePCR and the TAS method. Three RPVs from the TAS method were removed from the analysis due to pipetting error, as determined by primer utilization profiles. These included one replicate of RPV5 and both replicates of RPV15 (Table S1). We compared alpha and beta diversity analyses of the PCRs employing 15–18 unique RPVs to those generated with the fully degenerate primer set. All alpha and beta diversity analyses were performed on data rarefied to a depth of 1,800 sequences/sample (Table S1—experiment 3).

When employing fully degenerate primer pools, observed alpha diversity (Shannon index) of the fecal sample was slightly, but significantly higher when analyzed using the TAS protocol relative to the DePCR protocol (average Shannon index, five replicates, 2.71 to 2.66; ANOVA P < 0.001; Table 4). We then calculated average Shannon indices for analyses of the same gDNA sample with individual RPVs, employing TAS and DePCR protocols. The average Shannon index for the TAS reactions with unique RPVs (2.40) was significantly lower than that measured for the DePCR reactions (2.58) (ANOVA P < 0.001; Table 4). Finally, all RPV data, rarefied to 1,800 sequences per sample, was pooled together for TAS and DePCR approaches, independently. These combined datasets were then randomly sub-sampled to 1,800 sequences. These rarefactions were performed five times, and the average Shannon index for the combined RPVs was calculated. In this approach, average Shannon index from TAS (2.48) was significantly lower than for DePCR (2.69) (ANOVA P < 0.001; Table 4). Across all three methods of calculating observed diversity, there was no significant difference in measured Shannon index for the DePCR method (ANOVA, P = 0.377), while a significant decrease with each RPV independently was observed with the TAS method (ANOVA, P = 3.69e−8). When each RPV is used independently in the TAS protocol, the overall captured diversity is lower than with reactions with degenerate pools (Table 4) due to the greater number of potential mismatch interactions that can occur when a complex template is amplified with a single, non-degenerate primer. As the DePCR method is more tolerant of mismatches, no significant decrease in average Shannon index was observed. However, the observed variance in Shannon index among the individual RPVs was greater for the DePCR than for the TAS method (Table 4).

We next examined the structure of the observed fecal microbial communities in standard TAS and DePCR with degenerate primer pools, and with reactions conducted using RPVs (Fig. 5). We observed high reproducibility for five replicates using TAS (i.e., ‘TAS_pool’) or DePCR (i.e., ‘DePCR_pool’) with degenerate primer pools (Figs. 5A and 5B) and observed microbial community structure was significantly different between TAS and DePCR employing the degenerate primer pools (ANOSIM, R = 0.401, p = 0.001). Compared to amplifications with degenerate pools of primers, within-group variability was much greater for the analyses of RPVs individually with either amplification protocol (Figs. 5A and 5B, ‘TAS’ and ‘DePCR’). Within-group Bray–Curtis dissimilarity (BCD) of amplicons from the 15 (TAS) to 18 (DePCR) RPVs ranged from 0.03 to 0.36 for the TAS method and from 0.04 to 0.68 for the DePCR method (ANOVA P < 0.001; Fig. 5B). Conversely, the within-group BCD for five technical replicates generated with degenerate primer pools were 0.04 to 0.07 for TAS and 0.05 to 0.11 for DePCR (ANOVA P < 0.001). Profiles of the individual RPVs from DePCR analyses could be divided into two groups: (a) RPVs with profiles highly similar to degenerate primer pool analysis with either DePCR or TAS; and (b) RPVs with profiles divergent from the degenerate pool communities, and more similar to RPVs from TAS amplification reactions. Overall, the observed microbial community structure generated using the DePCR method with RPVs and with degenerate pools was not significantly different (ANOSIM R =  − 0.306, p = 0.99). Conversely, the observed microbial community structure generated using RPVs was significantly different that that observed with degenerate primer pools for the TAS method (ANOSIM R = 0.487; p = 0.003). Average BCD between TAS_pool and TAS RPV profiles (0.211) was significantly greater than for DePCR_pool and DePCR RPV (0.154) (ANOVA P < 0.001; Fig. 5C). DePCR BCD profiles were heavily weighted toward low dissimilarity, with a long tail of high dissimilarity comparisons. The long tail is a result of some primers generating highly divergent observed microbial communities with the DePCR protocol. Many of the primers which showed the poorest utilization within the degenerate pool (e.g., RPV10, 12, 14, and 18; node with red dot in Fig. 4A), generated the most divergent single RPV profiles. This suggests that these primers do not closely match the most dominant taxa within this particular gDNA sample.