PCR effects of melting temperature adjustment of individual primers in degenerate primer pools

Primer design and synthesis

The most recent Earth Microbiome Project (EMP) primers, 515F (Parada) [“515F”] and 806R (Apprill) [“806R”], were utilized as the default primer set (Parada, Needham & Fuhrman, 2016; Apprill et al., 2015; Walters et al., 2016; Quince et al., 2011). These primers are 4-fold degenerate (515F) and 24-fold degenerate (806R), with theoretical melting temperatures ranging from 66.9 to 71.8 °C (515F) and from 54.7 to 61.7 °C (806R) (Table S1). Primer theoretical melting temperatures were calculated using the OligoAnalyzer3.1 calculator (Owczarzy et al., 2008), assuming 250 nM primer concentration, 2 mM Mg²⁺, and 0.2 mM dNTPs. Maximum Delta G values for each sequence were calculated using the self-dimer option in the OligoAnalyzer software. Primers were synthesized either as single degenerate primer pools (standard approach), or as individual primers without degeneracies by Integrated DNA Technologies (IDT, Coralville, IA, USA). All primers were synthesized as LabReady and delivered at a fixed concentration of 100 micromolar. Most primers contained common sequence linkers (CS1 and CS2) at the 5′ ends, as shown in Table S2. These linker sequences are necessary for the later incorporation of Illumina sequencing adapters and sample-specific barcodes. For new primer pools containing shortened primers, a combination of either four (“515F”) or 24 (“806R”) non-degenerate primers were combined in equal volume to generate degenerate pools at 100 micromolar concentration. These primer pools were named “515F_pool” and “806R_pool”, and when these primers were used together, the primer set was named “ShortEMP”. One experiment, described below, employed ShortEMP primer pools without common sequences (“NoLinker_ShortEMP”) to assess the effect of the linker sequences on analysis of microbial community structure. A final set of primers, employing variable length spacers, was generated and named the “LongEMP” primer set. These primers included variable length spacers between the common sequence linkers and the gene-specific EMP primer regions, as described previously (Lundberg et al., 2013). Sequencing with these primers was performed on both MiniSeq and MiSeq instruments, but data from the MiniSeq using the “LongEMP” primer set did not properly merge due to read-length limitations, and was not further analyzed.

Genomic DNA templates

Four microbial genomic DNA (gDNA) samples were employed in this study. These include the ZymoBIOMICS Microbial Community DNA standard (D6306; Zymo Research, Irvine, CA, USA; ‘Zymo’), as well as three environmental samples derived from Lake Michigan sediment (‘LMC’), garden soil (‘Soil’) and mammalian (rat) feces (‘Feces’). Sediment, soil and fecal samples were extracted using a PowerSoil DNA extraction kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol.

Amplicon library preparation and sequencing

A two-stage PCR amplification strategy was used to generate sequencer-ready amplicons (Naqib et al., 2018). Genomic DNA was first PCR amplified with primer set EMP (CS1_806R and CS2_515F), ShortEMP (CS1_806R_pool and CS2_515F_pool), or LongEMP (CS1_806R_long_pool and CS2_806R_long_pool) (Table S2). All primers contained 5′ common sequence tags (known as common sequence 1 and 2, CS1 and CS2) as described previously (Moonsamy et al., 2013). First stage PCR amplifications were performed in 10 microliter reactions in 96-well plates, using MyTaq HS 2X mastermix (Bioline, Taunton, MA, USA). PCR conditions were 95 °C for 5 min, followed by 28 cycles of 95 °C for 30″, variable annealing temperature for 60″and 72 °C for 90″. For temperature gradients, annealing temperatures of 40°, 45° and 50 °C were employed. Subsequently, annealing temperatures of 45 °C were used for ShortEMP reactions and 50 °C for EMP reactions.

Second stage PCR amplifications were performed in 10 microliter reactions in 96-well plates. 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 (Fluidigm, South San Francisco, CA; Item# 100-4876). These AccessArray primers contained the CS1 and CS2 linkers at the 3′ ends of the oligonucleotides. Cycling conditions were as follows: 95 °C for 5 min, followed by 8 cycles of 95 °C for 30″, 60 °C for 30″and 72 °C for 30″.

In an experiment using the NoLinker_ShortEMP primers, the above protocol was modified slightly. To assess the effect of common sequence linkers on the observed microbial community structure, we performed the first stage PCR amplification with NoLinker_ShortEMP primers. PCR conditions were identical to those described above. During the second stage PCR amplification, both Fluidigm Access Array Barcode and the ShortEMP primers were included, and 12 cycles of amplification were performed in place of eight. In this approach, common sequences are incorporated during the second stage PCR as the ShortEMP primers with linkers amplify the NoLinker_ShortEMP amplicons; subsequently, the Fluidigm Access Array barcode primers amplify amplicons containing the common sequence linkers.

In all experiments, samples were pooled using an EpMotion5075 liquid handling robot (Eppendorf, Hamburg, Germany). Pooled libraries were purified using an AMPure XP cleanup protocol (0.6X, vol/vol; Agencourt, Beckmann-Coulter) to remove fragments smaller than 300 bp. The pooled libraries, with either a 1% or 20% phiX spike-in, were loaded onto either an Illumina MiniSeq mid-output kit (2 × 153 paired-end reads) or an Illumina MiSeq V2 kit (2 × 250 paired-end reads), as indicated in Table S3. Fluidigm sequencing primers, targeting the CS1 and CS2 linker regions, were used to initiate sequencing. De-multiplexing of reads was performed on instrument. Library preparation, pooling, and sequencing were performed at the University of Illinois at Chicago Sequencing Core (UICSQC).

Bioinformatic analysis of sequence data

Raw FASTQ files were downloaded from Illumina Basespace. Sequence reads were merged using PEAR (Paired-End Read Merger) (Zhang et al., 2013) with default parameters. Merged reads were quality trimmed (<Q20 discarded) and length trimmed (<250 bases were removed), and primer sequences were removed using the ‘Trim reads’ algorithm within the software package CLC Genomics Workbench (v11; Qiagen; Hilden, Germany). Trimmed sequences were reverse complemented using a QIIME script. Chimeras were removed using the USEARCH81 algorithm (Edgar, 2010). Subsequently, for environmental samples (LMC, Soil, Feces), sequences were pooled, renamed and clustered into operational taxonomic units (OTUs) at a threshold of 97% similarity (QIIME v1.8.0; Caporaso et al., 2010). Each OTU was annotated taxonomically based on the representative sequences using the UCLUST algorithm and the greengenes 13_8 reference database (McDonald et al., 2012b). A biological observational matrix (BIOM) was generated from the clustered OTU data and the taxonomy data (McDonald et al., 2012a). For subsequent analyses, data were rarefied using the vegan package within the R programming language. Rarefaction depths were adjusted by analysis; depth of rarefaction for each figure is shown in Table S3. The BIOMs were analyzed and visualized using the software package Primer7 (Clarke, 1993) and in the R programming environment (R Core Team, 2013). Dendrogram creation and SIMPROF tests were conducted within Primer7. The vegan package (Oksanen et al., 2011) was used to generate alpha diversity indices and to calculate pairwise Bray–Curtis dissimilarity scores. Analysis of similarity (ANOSIM) calculations were performed at the taxonomic level of genus, using square root transformed data. Metric multi-dimensional scaling (mMDS) plots were created using the cmdscale and ggplot2 functions within R. Ellipses, representing a 95% confidence interval around group centroids, were created assuming a multivariate t-distribution (Fox & Weisberg, 2011). Taxon-level differential abundances between sample groups were identified using the software package STAMP  (Parks et al., 2014) employing White’s non-parametric t-test (White, Nagarajan & Pop, 2009). P values were adjusted using the Benjamini–Hochberg False Discovery Rate correction  (Hochberg & Benjamini, 1990).

Sequence data from ‘Zymo’ samples were processed using the same pipeline, but data were not clustered. Instead, merged, trimmed and chimera-cleaned data were mapped against reference gene sequences for the eight bacterial reference organisms using the software package CLC genomics workbench v10 (Qiagen, Aarhus, Denmark). The eight reference organisms included: Lactobacillus fermentum (AJ575812), Bacillus subtilis (DQ993674), Escherichia coli (J01859), Enterococcus faecalis (EU887827), Salmonella enterica (JQ694167), Listeria monocytogenes (M58822), Pseudomonas aeruginosa (LN874213), and Staphylococcus aureus (L37597). Mapping data were then converted to biological observation matrices for use in visualization and statistical analyses, as described above. Ideal scores were calculated according to the formula described previously (Green, Venkatramanan & Naqib, 2015). Briefly, the Ideal score is a summation of the absolute difference between the expected relative abundance and the observed relative abundance for each feature in a multi-feature dataset.

Data sharing

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

Primer design

The design of new primers with lower variation in melting temperature was performed with two constraints, including: (a) no primer shorter than 16 bases, and (b) all degenerate positions are retained. Outside of these constraints, we aimed to minimize theoretical melting temperature variation within each degenerate primer pool. The high overall GC content and 3′ run of A or T bases in the 515F primer limited our ability to adjust Tm of this primer while meeting the above constraints. An underlying assumption of the primer modification was that shortening the primers could only increase the range of potential targets (Isenbarger et al., 2008). To adjust primer Tm, bases were sequentially removed from either or both 5′ and 3′ ends of each oligonucleotide (Fig. 1A; Table S1). The final primer design (“ShortEMP”) yielded pools of primers with overall lower average theoretical melting temperatures, and with smaller range in Tm between all primers within the pool compared to the original pool (“EMP”) (Fig. 1B; Table S1). The 515F primer pool, comprised of four primers, was only modestly affected by the primer alterations. The Tm range prior to modification was 66.9 to 71.8 °C; after modification, the Tm range was 66.9 to 69.1 °C. The Tm range of the modified primer pool was not significantly different from that of the original pool [Wilcoxon–Mann–Whitney (WMW) test; p = 0.37]. The 806R primer pool, comprised of 24 primers was more strongly affected by the alterations. The Tm range prior to modification was 54.7 to 61.7 °C; after modification, the Tm range was 53.6 to 56.6 °C. The Tm range of the modified primer pool was significantly different from that of the original pool (WMW test; p < 0.00001).

A secondary aim of design was to introduce nucleotide diversity at the 5′ end of the combined oligonucleotides by preferential removal of bases as the 5′ end of the gene-specific portion of each primer (Fig. 1A; Tables S1, S4 and S5). Illumina’s technical notes indicate that having relatively even proportions of all four nucleotides during each sequencing cycle is necessary for proper sequencing. This is particularly true during the first cycles of sequencing, when the sequencer is still training itself to identify clusters. To improve sequence quality of so-called “low diversity” libraries (e.g., 16S rRNA gene amplicon libraries), exogenous spike-ins of shotgun DNA libraries derived from the virus phiX174 are typically employed. Elsewhere, Lundberg et al. (2013) introduced nucleotide diversity through the synthesis of primers containing frameshift nucleotides. Briefly, a mixture of six different forward and six different reverse primers were pooled together (Lundberg et al., 2013). Each primer variant contained an identical gene-specific degenerate primer at the 3′ end of the oligonucleotides, but a variable number of nucleotides upstream of the gene-specific portion of the primers. This variable number of upstream nucleotides introduces artificial diversity into the sequencing reaction due to an effective frame shift by offsetting the start of highly conserved regions of the amplicon by 1–5 bases. We aimed to achieve the same effect by removing bases from the 5′ end of the gene-specific portion of the primers, thereby introducing the frame-shift effect (Fig. 1A). To assess the effect of 5′ nucleotide removals on the overall nucleotide diversity at each position, we generated two in silico calculations: (a) number of different nucleotides present at any given position across the first 16 nucleotides of each primer pool, and (b) the Shannon index for nucleotide diversity at each of the first 16 nucleotides of each primer pool (Tables S4 and S5). The original primer pool design has very low nucleotide diversity, as the primers are perfectly synchronized and have only two or three degenerate positions. Thus, for the original 515F primer pool, all primer positions other than the 4th and 9th positions have only a single nucleotide represented. The average number of nucleotides present at any position is 1.125, and the average Shannon index (natural log) across the first 16 positions is 0.09 (Table S4). Conversely, after modifying the individual primers, the average nucleotide diversity across the 16 positions was 2.44, with an average Shannon index of 0.77. The maximum possible Shannon index value for 4 different nucleotides is 1.39. A similar effect was observed for the 806R primer pool (Table S5). Prior to redesign, the average nucleotide diversity across the first 16 positions was 1.38, compared to 3.19 after redesign. Likewise, Shannon index increased from 0.2 to 1.0 after redesign.

Microbiome profiling using EMP and ShortEMP primers

To assess the effects of primer modifications, initially a single complex genomic DNA template (LMC) was profiled across an annealing temperature gradient using EMP and ShortEMP primer sets. We previously observed a strong and consistent shift in microbial community structure associated with increasing annealing temperatures from 40 °C to 55 °C, with specific taxa such as Prevotellaceae increasing in relative abundance with increasing annealing temperature in mammalian fecal samples (Green, Venkatramanan & Naqib, 2015). Therefore, annealing temperatures of 40 °C, 45 °C and 50 °C were tested in this study. The single gDNA template was PCR amplified at each of the temperatures and primer sets with six technical replicates. The observed microbial community from the source gDNA template was significantly affected by primer set (EMP or ShortEMP) and by annealing temperature (Fig. 2). Variation in microbial community structure by temperature was described primarily by the MDS axis 1, while that of primer set was described primarily by MDS axis 2. The magnitude of the shift was small, with overall Bray–Curtis similarity of all comparisons between EMP and ShortEMP replicates >0.84 (Fig. 3A). Within each primer set, replicates had Bray–Curtis similarity values >0.89, with ShortEMP replicates having slightly and significantly greater similarity (Fig. S1, LMC). No significant effect of primer set and annealing temperature on alpha diversity indices (e.g., richness and Shannon index) at the genus- or OTU-level were observed for LMC (Table 1). At each annealing temperature, relatively few taxa (16–28 genus-level taxa) were significantly differently abundant between the EMP and ShortEMP analyses (Fig. 2B, Fig. S3).

Based on the results of the temperature gradient analysis (above) and testing of a mock community standard (below), an annealing temperature of 45 °C was chosen for additional analyses employing the ShortEMP primers and 50 °C for EMP primers. Two additional complex microbial samples were analyzed at these temperatures, with six technical replicates for each. These samples include ‘feces’ and ‘soil’. In each of these analyses, a significant effect of primer set was observed (Figs. 3B and 3C; ANOSIM R = 1, p < 0.008), and in each case, the magnitude of the effect was similar, with overall Bray–Curtis similarity between all replicates of the same sample >0.82 (Figs. 3B and 3C). Small effects of primer set on alpha diversity indices were observed, and only the OTU-level Shannon index for the soil sample was significantly different between ShortEMP primers relative to EMP primers (7.09 vs 7.16) (Table 1).

Interrogation of a mock community with EMP, ShortEMP, LongEMP and NoLinker_ShortEMP primer sets

The ZymoBIOMICS Microbial Community DNA standard (“Zymo”) was used as a gDNA template to assess the capability of the various primer sets to characterize microbial communities. The Zymo standard is composed of 8 bacterial taxa and 2 fungal taxa at varying levels of abundance (fungi are not amplified by the 515F/806R primer set). When analyzed using 16S rRNA gene amplicon sequencing, the relative abundance of each strain should range from 4.2% to 18.4%. Both primer sets (EMP and ShortEMP) generated highly similar but significantly distinct results (Fig. 3D; ANOSIM R = 1, p = 0.0019). For each replicate, an Ideal Score (Green, Venkatramanan & Naqib, 2015) was calculated based on the expected relative abundance of each taxon. The Ideal Score represents a summation of the difference in relative abundance of each taxon in a mock community. A perfect representation of the expected relative abundance of each taxon will yield an ideal score of zero, with values greater than zero representing increasing discordance from the expected results. Ideal score results were generated for each of the six technical replicates and compared at annealing temperatures of 40 °C, 45 °C and 50 °C for EMP and ShortEMP primers (Table 2). When employing the EMP primers for analysis of the Zymo standard, annealing temperatures below 50 °C produced progressively worse results (i.e., higher ideal scores). When employing the ShortEMP primers for analysis of the Zymo standard, annealing temperatures of 45 °C and 50 °C yielded similar results, with low ideal scores relative to 40 °C (Table 2).

Ideal scores indicated that when used at their optimal annealing temperature EMP primers were slightly, but significantly, better than ShortEMP primers in recovering the expected distribution of the Zymo bacterial taxa. Although the EMP primers slightly outperformed the ShortEMP when compared to the expected distribution of eight bacterial taxa, the ShortEMP primers better tolerated a broader range of annealing temperatures. ShortEMP primers produced nearly identical results at 45 °C and 50 °C (Table 2), demonstrating that the lower overall Tm of the primer pool can contribute to a shift in optimum annealing temperature.

We sought to assess whether the addition of spacer regions, as described previously by Lundberg et al. (2013), would alter the efficiency of amplification in this system. We synthesized primers containing common sequence linkers, a 0–5 base frameshift sequence, a two-base ‘linker’, and the EMP primer sequence (Table S2; ‘LongEMP’). A pool of six forward primers and six reverse primers were used in the standard two-stage PCR protocol to amplify the mock community gDNA. The observed mock community generated from LongEMP analyses was most similar, and not significantly different, from that generated using the EMP primer set (Fig. 4; ANOSIM R =  − 0.179, p = 0.965). EMP and LongEMP primer sets generated more similar observed microbial communities, relative to the ShortEMP primer set (ANOSIM R = 0.715 − 0.722, p = 0.002). Although significant, the difference was not large, and EMP, ShortEMP and LongEMP generated highly similar results when applied to the mock community DNA standard (Fig. 4).

We also sought to assess whether the common sequencer linkers themselves contributed to distortion of the underlying microbial community structure. Conceptually, the 3′ ends of the linker sequences could interact with genomic DNA templates, leading to preferential amplification of some templates. To avoid initial interaction of linker sequences with gDNA templates, we performed the first stage PCR amplification with ShortEMP primers without common sequence linkers (NoLinker_ShortEMP). After 28 cycles of amplification, the generated PCR amplicons were transferred to the second stage PCR amplification, as performed for all other reactions. To allow these amplicons without linkers to be prepared for Illumina sequencing, both ShortEMP primers (with common sequence linkers) and Fluidigm Access Array barcoding primers were added to the second stage reaction. Final amplicons were sequenced and the resulting sequence data analyzed together with EMP, ShortEMP and LongEMP primer sets (Fig. 4). The observed community generated by the NoLinker_ShortEMP primer set was significantly different than those generated by EMP, ShortEMP and LongEMP primer sets (ANOSIM R = 1, p = 0.002), and the observed community structure was also more divergent from the expected structure than those generated with the other primer sets (Fig. 4). Thus, the common sequence linkers do not appear to substantially alter the observed mock microbial community structure, and removing the linkers leads to a more complex workflow and a poorer representation of the mock community.

We should note that the Zymo standard is not an ideal mock community for assessing the action of a degenerate primer pool. We sought to identify the sequences of the eight bacterial DNA templates in the region of the 515F and 806R primers, and despite the moderately broad range of bacterial taxa included (Proteobacteria, Actinobacteria, and Firmicutes), all sequences were identical at the primer annealing locations. This mock community is, therefore, appropriate for determining that an overall PCR and sequencing workflow is successful but is not a good approximation of a highly complex natural microbial community, where even conserved sites such as the 515F and 806R primer sites contain considerable sequence heterogeneity.

We also used the amplicons generated from the Zymo standard to assess variability introduced by sequencing on the Illumina MiSeq or MiniSeq platforms. For the EMP, ShortEMP, and NoLinker_ShortEMP, identical amplicons were sequenced and analyzed on both platforms. The LongEMP amplicon reads generated on the Illumina MiniSeq did not consistently merge, and were not included in the analysis. This is due to the additional ‘spacer’ bases (up to six on each read) and the short overlap of paired-end reads generated with 2 × 153 bases reads on the MiniSeq. MiSeq and MiniSeq results from the same amplicons clustered together, regardless of primer set (Fig. 4). For the EMP and ShortEMP primer sets, the effect of sequencing platform was not significant (ANOSIM R = 0.166, p = 0.091 and R = 0.05, p = 0.281, respectively). For the NoLinker_ShortEMP primer set, a significant effect of sequencing platform was observed (ANOSIM R = 0.416, p = 0.002). Thus, a slight effect of sequencing platform was observed, but the choice of platform did not alter biological conclusions.

Assessing the need for phiX spike-in with ShortEMP primers

We sequenced ShortEMP amplicons on an Illumina MiniSeq run without substantial phiX spike-in to determine if ShortEMP amplicons by themselves produced sufficient nucleotide diversity to allow proper clustering. These amplicons, generated from Lake Michigan sediment, at annealing temperatures of 40, 45 and 50, were analyzed together with the same amplicons generated on a MiniSeq run with a 20% phiX spike-in. The overall quality of the run was extremely high (>96% pass-filter, >97% Q30 with approximately 1.75% phiX spike-in), and results were similar to those generated on a 20% phiX run (Fig. 5). Slight trends towards small, significant differences were observed between the two sequencing runs, but the overall magnitude of the difference was small.