Thermal-bias PCR: generation of amplicon libraries without degenerate primer interference

Reaction modeling

Modeled PCR data were generated using Excel (Microsoft). Input cells referenced for the model included: “seed”, “max”, and “K_D”. A “Maximum efficiency” reference cell converted a desired reaction efficiency to a yield relative to a 2-fold amplification in the model (=(efficiency * 2) −1 ), such that an efficiency input of 1.0 generated a value of 1.0, and an efficiency input of 0.98 generated a value of 0.96 (98% of 2-fold). The yield of the first cycle referenced the seed input cell, and the yields of subsequent cycles referenced the value in the cycle preceding them. An example of a model cell is “=B6*($L$9+(($G$6-B6)/$G$6)-(B6/($J$6+B6)))”, where B6 is the preceding cell, $L$9 is the efficiency cell, $G$6 is the max cell, and $J$6 is the K_D cell.

PCR template preparations

Escherichia coli genomic DNA was purchased from Fisher Scientific (cat. AAJ14380MA), resuspended in DNA buffer (10 mM Tris-Cl, 0.1 mM EDTA, pH 8.0) and quantified using UV absorbance (260 nm). Serial dilutions were prepared in DNA buffer. A mock community bacterial genome mixture ‘ABRF-MGRG 10 Strain Even Mix Genomic Material’ was purchased from the American Type Culture Collection (ATCC, cat. MSA-3001) (Foox et al., 2021). It contains a mixture of Bacillus subitilis (Bs), Chromobacter violaceum (Cv), Escherichia coli (Ec), Entercoccus faecalis (Ef), Halobacillus halophilus (Hh), Haloferax volcanii (Hv), Micrococcus luteus (Ml), Pseudomonas fluorescens (Pf), Pseudoalteromonas haloplanktis (Ph), and Staphylococcus epidermidis (Se) genomes. This mixture is listed by the ATCC as having an even mixture of each genome and a company representative indicated that it contains an equal weight of each. However, the creators of this mixture prepared it as having an even concentration of each genome and they subsequently determined that the genomes may not be evenly represented (Foox et al., 2021). The match and mismatch reporter templates were generated using PCR with NEBNext Ultra II Q5 polymerase (New England Biolabs, cat. M0544S) by amplifying segments of the E. coli V3-V4 region of the 16S rRNA gene using primers that introduced restriction sites as insertions and also modified the targeting primer-binding sites in the mismatch version. Overlap extension PCR (Q5 polymerase) was then used to fuse the two segments. The final products were gel purified and the concentrations were determined by UV absorbance in a SynergyMX spectrophotometer (Biotek Instruments) with a sample diluted 10-fold in TE buffer in a 1 cm cuvette. Yields were ~0.2 μg/μL and had a 260/280 ratio of 1.85. Templates were then normalized to 0.1 μg/μL (0.31 μM), and stored in DNA buffer.

Primer and protocol performance assays

E. coli genomic DNA was diluted to 1.95 μg/mL and 1 μL was used per 20 μL of reaction mixture (97.5 pg/μL final). Platinum II Taq and Ultra II Q5 reactions were supplemented with EvaGreen dye (Jena Bioscience, cat. # PCR-379); iTaq Universal SYBR Green (Bio-Rad, cat. # 1725120) and SsoFast EvaGreen (Bio-Rad, cat.172-5201) were used without modification. The Earth Microbiome Project (EMP) ‘16S Illumina Amplicon Protocol’ uses a 45 s melt at 94 °C, 60 s anneal at 50 °C, and a 90 s extension at 72 °C (https://earthmicrobiome.org/protocols-and-standards/16s/) (Ul-Hasan et al., 2019; Gilbert, Jansson & Knight, 2014; Thompson et al., 2017). The evaluated thermal-bias protocols included an enzyme activation step (2 min at 98 °C), two targeting cycles (10 s melting at 98 °C and annealing/extension ranging from 2.5–10 min at 46–55 °C), and 30 or 40 amplification cycles (10 s melting at 98 °C and annealing/extension for 1 min ranging from 74–84 °C). Fluorescence measurements were taken at the end of each amplification cycle. These protocols included a melt-curve analysis at the end that progressed from 75–95 °C in 0.5 °C steps.

Primer sequences are listed in Data SD1 and were purchased from Eurofins Genomics as ‘NGS Grade’ and stored as 100 μM stocks in DNA buffer. Primer pairs were mixed in High-Performance Liquid Chromatography (HPLC)-grade water at 1.11 μM each and added to reactions at 45 % (0.5 μM final). Standard PCRs were performed in a Bio-Rad MJ Mini thermal cycler and qPCRs were performed in a Bio-Rad CFX96 Touch Real-Time system. For global fitting, fluorescence data were exported from the real-time system and saved in .csv format before being analyzed using an online qPCR fitting program (http://www.bioinformatics.org/ucfqpcr/) (Carr & Moore, 2012). The seed, max, and K_D values were then transferred to Excel (Microsoft) for subsequent processing. Figure data were plotted using Prism 9 (Graphpad Software) and annotated using Illustrator (Adobe).

Reporter template amplification and restriction digestion

Conventional PCRs were performed using Platinum II Taq polymerase and qPCRs were performed using SsoFast EvaGreen unless otherwise indicated. The reporter templates were diluted 10,000-fold in DNA buffer and 1 μL was used per 20 μL of reaction mixtures. DNA from those PCRs was purified using an in-house PEG precipitation protocol: mixing with an equal volume of Polyethylene Glycol (PEG) precipitation reagent (20% PEG-8000, 1 M NaCl); incubated at room-temperature for 10 min; centrifuged at 14,000 RCF for 15 min; liquid removed; pellet washed once with 75% ethanol, centrifuged 10 min, liquid removed; washed once with 95% ethanol, centrifuged 5 min, liquid removed; air dried. DNA was resuspended in a volume of DNA buffer matching the original PCR reaction volumes. Separate gel electrophoresis evaluations of this PEG procedure indicated that a single round reduced a 2 μM spike of primers to a barely detectable level with a complete recovery of DNAs larger than ~300 bp. As with magnetic bead-adhesion protocols, the size cut-off of recovered DNAs can be increased by lowering the added fraction of PEG reagent.

Restriction digestions were performed at 37 °C for 120 min using either BamHI-HF or SpeI-HF (NEB, cat. R3136S and R3133S) according to the manufacturer’s instructions. In preliminary restriction digestion studies, it was noted that neither the BamHI or SpeI digestions went to completion. Suspecting there may have been contamination by templates lacking the restriction sites, BamHI- and SpeI-cleaved amplicon fragments were independently gel purified and then ligated back together to ensure the final templates contained the restriction sites prior to Q5 amplification. However, these preparations were also unable to be fully digested. Digestion kinetics were also evaluated over 24 h and using increased enzyme concentration and using a different lot of enzyme. Although the digestion progression was completed within 30 min, a small amount of residual undigested material persisted. Amplicon reannealing was performed by heating samples to 95 °C for 5 min in a thermal cycler and subsequent cooling to 25 °C.

Sequencing library preparation and analysis

Each 100 μL PCR mixture contained 2.5 μL of ATCC MSA-3001 genome mixture (~0.175 ng/μL final) and either 2.5 μL of DNA buffer (Library_1), a 1:1 reporter template mix at a 1/100,000 dilution (Library_2), or a 1/200,000 reporter dilution (Library_3). The thermal-bias primers TB_F_Tm55_v2 and TB_R_Tm55_v2 contained Illumina TruSeq adapter sequences and had calculated targeting Tms of 55 °C and amplification Tms of 80 °C. SsoFast EvaGreen mixtures were distributed into four 20 μL reaction wells and the thermal-bias protocol used two targeting cycles (annealed 5 min at 48 °C) and 30 amplification cycles (annealed 1 min at 78 °C). After cycling, reaction wells were pooled and the amplicons were purified using two rounds of PEG precipitation. Dried amplicons were resuspended in 60 μL of DNA buffer, quantified using a Quant-iT PicoGreen assay (Invitrogen), and normalized to 20 ng/μL. The normalized samples were evaluated using gel electrophoresis to confirm amplicon size and primer removal. Bar-coding and 250 nt Illumina sequencing was performed by a commercial service (Azenta Life Sciences).

Sequence reads were processed and mapped on a public Galaxy server (usegalazy.org) (Afgan et al., 2016): Cutadapt v5.0 (Martin, 2011) was used to remove reads having a quality scores less than 20 and lengths less than 200 nt, Bowtie2 v2.5.3 (Langmead & Salzberg, 2012) was used to map reads onto a reference FASTA file containing concatenated V3-V4 regions of each bacterium and the match and mismatch reporter sequences. Alignment maps and their indices were downloaded and visualized using JBrowse2 (Diesh et al., 2023), which was also used to determine read depths at each V3-V4 region.

Primer degeneracy reduces PCR quality

We have used global fitting of qPCR data during other projects to establish relative template abundances (seed values) (Carr & Moore, 2012; Smith et al., 2019; Bose et al., 2021). Inspection of the associated max and K_D values from those fitting operations revealed that the max/K_D ratios for a given reaction were reflective of the amplification quality, with more ‘upright’ reaction curves (indicating more productive reactions) having lower ratios. With an interest in using max/K_D ratios as a metric to assess reaction quality, modelling was used to establish the impact of varying max and K_D terms on PCR profile shapes (Fig. S1). That exercise revealed that reaction ’poisoning’ is primarily responsible for inefficient amplification and that end-point yield measurements are unreliable indicators of reaction quality.

The amplification performance of degenerate primer pairs was evaluated by monitoring the accumulation V3-V4, V4-V5, or V4 amplicons from Escherichia coli genomic DNA under reaction conditions recommended by the Earth Microbiome Project (EMP). Primer names, sequences, and calculated Tms are listed in Data SD1. In each target case, the performance of the degenerate primer pair was considerably poorer than the non-degenerate pair, yielding higher max/K_D ratios, lower overall yield, and lower seed values (Fig. 1). It is notable that the greatest performance disparity occurred during amplification of the V3-V4 region, in which the forward primer contained degeneracy at the −4 position relative to the 3′ end. It is well documented that polymerases are less tolerant of mismatches near the 3′ end (Bru, Martin-Laurent & Philippot, 2008; Wu, Hong & Liu, 2009; Gohl et al., 2021), and structural studies revealed that positions −1 through −4 interact with the DNA duplex binding region of polymerases and mismatches in those locations disrupt the active site (Johnson & Beese, 2004). Therefore, the location of this degenerate position likely exacerbated poor performance.

The relative performance of degenerate primers was also evaluated using other polymerases and reaction conditions. In each case, the degenerate pair performed poorly compared to the non-degenerate pair (Fig. S2). Importantly, even though the degenerate sets contained primers that perfectly matched the initial targets, the reaction performances were compromised before the amplicons substantially accumulated above the baselines, when less than 1% of the final amplicon pools had been generated. Therefore, limiting cycle numbers would not avoid this performance disparity.

Templates containing mismatches at priming sites

To evaluate primer performance on target templates containing mismatches to non-degenerate primers, an engineered V3-V4 reporter template was constructed that contained two ‘mismatch’ nucleotide changes in each of the forward and reverse primer binding sites of the target (Fig. 2A). These template changes are in positions predicted to be compensated by these degenerate primers (Fig. S3); yet, they represent a rather exotic case for the amplification of V3-V4 regions, as non-canonical eubacterial 16S rRNA genes from environmental samples usually contain only one or two mismatches to the consensus sequences at those positions, with rarer examples having highly deleterious mismatches farther into the 3’ primer-binding region (Gohl et al., 2021). A V3-V4 control template was also constructed that contained consensus ’match’ primer-binding sites. With an eye toward evaluating the amplification of mixed template samples, each engineered template additionally contained a unique restriction enzyme site (SpeI in the ’mismatch’ template, BamHI in the ’match’ control) so that amplicon pools could be evaluated for the presence of products derived from either template using gel electrophoresis.

When monitored by qPCR, the non-degenerate primers again performed better than the degenerate pair on either template (Fig. 2B). In reactions containing a 1:1 mixture of each template, amplicons derived from the match template dominated using either primer pair and products derived from the mismatched template were undetectable by restriction digestion (Fig. 2C). Therefore, the degenerate primers did not compensate for the template mismatches and they impeded amplicon production from both templates. Taken together, these data suggested that using non-degenerate primers to amplify heterogeneous targets could be a fruitful strategy if the performance on mismatched targets could be improved.

Thermal-bias PCR for co-amplification of a mismatched template

Commonly-employed V3-V4 degenerate primer pairs have predicted annealing Tms that are substantially different (Data SD1). In pilot experiments, it was found that using non-degenerate primers with matched Tms improved the co-amplification of the mismatched reporter to detectable levels (Fig. S4). Following these insights, a new amplification strategy termed “thermal-bias PCR” was developed that uses non-degenerate consensus primers having matched targeting Tms and containing long 5′ tails harboring sequencing adapters and additional nucleotides to substantially raise and match their full-length Tms (Fig. 3A). The thermal-bias protocol employs a low targeting annealing/extension temperature for only the first two cycles (during which the requisite complementary amplicon strands are generated). The cycling protocol then switches to use a very high annealing/extension temperature for the remaining amplification cycles, which forces the reactions to only use amplicons as templates instead of continuously re-targeting the initial templates. Thus, any bias against priming a particular target would only be present for two cycles, with the remainder of the cycles ‘locking in’ relative abundance.

An additional protocol improvement came from the use of an alternative polymerase. SsoFast is a Pfu/Vent hybrid polymerase fused to the Sso7d DNA binding protein. The Sso7d domain not only improves the processivity (extension rate) of polymerases to which it is attached, but it also increases primer annealing Tms (Wu et al., 2017; Wang et al., 2004). Although Taq polymerase functioned in trial thermal-bias protocols, it was found that that SsoFast polymerase tolerated degenerate primers and mismatched templates better (Fig. S5), so it was selected for further studies.

The amplification stages could be performed as high as 80 °C without substantially inhibiting the reactions (the max/K_D ratios were similar from 74 °C to 80 °C,) (Fig. S6A). The thermal-bias against re-priming on the initial targets during amplification was confirmed using protocols that lacked the initial targeting steps, which failed to yield amplicons (Fig. S6B). Altering the temperatures or times of the targeting cycles had a negligible impact on performance as long as that temperature was at least ~4 °C lower than the calculated targeting Tm. Thus, a thermal-bias PCR protocol can be employed that has an annealing temperature difference of ~25–30 °C between the targeting and amplification stages.

Using a thermal-bias protocol, amplification of the mismatched V3-V4 template appeared nearly as efficient as the match template and amplicons derived from an even mixture of the two templates contained both versions after 40 PCR cycles (Figs. 3B and 3C). It was noted that qPCR data from reactions using mixed reporter templates exhibited an apparent lag in amplification (e.g., Fig. 3B ‘Mix’), which was puzzling because amplicon strands were presumably being independently primed. The signals generated in those experiments were dependent on a fluorescent dye binding to dsDNA, so one explanation might be that heterogeneously-annealed amplicons (containing a match and mismatch strand) might not have bound as much dye as homogeneous amplicons. Heterogeneous amplicons would also not be digestible by SpeI, which would lower the apparent abundance of the mismatch product. To evaluate this idea, homogeneous match and mismatch amplicons were mixed in equal proportions before digestion. An aliquot of that mixture was then thermally melted and reannealed before each sample was subjected to SpeI digestion (Fig. 3C, lanes 5 and 6). It was evident that reannealing the mixture reduced the abundance of SpeI-digestible amplicons, indicating that heterogeneity in the digestion assays had lowered the abundance estimates.

Thermal-bias PCR accurately reports changes in template abundance

A common goal of amplicon deep-sequencing is to measure changes in the relative abundance of organisms, which requires that invariant populations be consistently reported and that changes in other populations be accurately reported (Morton et al., 2019). To determine the utility of thermal-bias PCR for generating sequencing libraries, the V3-V4 reporter match and mismatch templates were spiked into a commercial mixture of ten bacterial genomes at two different levels (one having half as much spiked as the other). These mixtures were then subjected to a 30-cycle thermal-bias PCR using primers that appended Illumina adapters (Fig. 4A). Replicate samples were pooled, cleaned, and normalized prior bar coding and paired-end sequencing. Sequence reads were then mapped onto a reference sequence that contained each bacterial V3-V4 region as well as the match and mismatch reporter sequences. In this experiment, the ‘reverse’ reads covered the restriction sites in the reporters and allowed them to be unambiguously mapped relative to amplicons derived from E. coli genomes. Read depths were divided by the 16S rRNA gene copy number for each organism and plotted as a percentage of total (Fig. 4B).

The relative abundances of the genome-derived V3-V4 amplicons was highly similar among the libraries and, for the most part, aligned with several independently measured genome compositions of this mixture (Foox et al., 2021). Nine of these genomes had consensus primer-binding sites and the relative abundances were not correlated with the G/C content of their V3-V4 regions. One notable exception was a substantial underrepresentation of the Archaeon H. volcanii, which has a high G/C content in its genome, but not in the V3-V4 region. This underrepresentation was unexpected because the primer-binding sites of H. volcanii presented three mismatches: two were the same as those in the engineered mismatch template (G to C at −8 in the forward site, and T to G at −14 in the reverse), with the third H. volcanii mismatch being outside of a 3′ region (A to G at −9 in the forward site). This V3-V4 region has as other properties that likely contributed to an inefficient targeting stage (discussed below).

Both the match and mismatch reporter amplicons were present in the spiked libraries, but the mismatch product was reduced by ~37% compared to the match product in each case. Importantly, each reporter version was reduced in abundance by ~50% in Library_3 compared to Library_2, accurately reflecting the difference in spiking level. Taken together, any biases in the targeting stages were preserved through the amplification stages along with apparent differences in relative template abundance. Therefore, the thermal-bias protocol is a promising alternative method for generating sequencing libraries and it can tolerate non-consensus targets.