Comparing DNA Extraction and 16S rRNA Gene Amplification Methods for Plant-Associated Bacterial Communities

Plant-associated microbial communities (the plant “microbiome”) play a significant role in global ecology and agriculture. These microbes affect how plants grow, acquire nutrients, defend against pests and disease, and otherwise carry out many essential activities (Compant et al. 2019; Müller et al. 2016). There is currently great interest in understanding how plants and microbes interact and how these interactions could be harnessed to improve human activities (Bell et al. 2019; Reid and Greene 2013).

Although the plant microbiome consists of bacteria, archaea, fungi, oomycetes, protists, and other microscopic organisms, most current research focuses on bacteria, followed by fungi. The most common method to profile the bacteria in plant-associated samples is targeted sequencing of the 16S ribosomal RNA (rRNA) gene. (An alternative method, whole-metagenome sequencing, is rarely used when plant tissue is involved because of how hard it is to separate microbial DNA from the much more abundant host DNA.) This method consists of extracting DNA from plant-associated communities, using PCR to amplify a specific gene sequence (in this case, 16S rRNA) from the bacterial community, sequencing the resulting amplicon, and performing bioinformatic analyses to identify patterns of community composition. The two major processing steps in amplicon analysis—DNA extraction and PCR amplification—are known to be subject to bias based on the specific methods and conditions used (Pollock et al. 2018). Although no method is completely free of bias, researchers understandably want to use the methods with as little bias as possible.

Bias during DNA extraction comes from incomplete recovery of DNA from the sample. This can be due to microbial cells that are not lysed at equal efficiency or due to chemical inhibitors present in the sample (Pollock et al. 2018). Although phase-separation methods (such as cetyltrimethylammonium bromide extraction) can be used to recover plant-associated microbial DNA, many researchers choose to use commercial column- or bead-based purification kits due to their greater convenience and consistency. However, these kits have generally been optimized either for the plant itself or for microbes in other environments (soil, feces, and so on). Extracting plant-associated microbes includes the combined challenges of lysing a diverse array of microbes while also dealing with potential chemical inhibitors from the plant itself, and it is unknown how well different kits perform under these circumstances.

Bias during PCR amplification comes primarily from PCR primers and bias from the DNA polymerase itself (meaning, the enzyme does not amplify all targets equally well). Plant-associated microbiome samples pose an additional challenge in that the most common primers to amplify microbial DNA also amplify plant organelles (chloroplasts and mitochondria). Because organelles are generally much more abundant than microbes in plant tissue, most sequencing reads come from host DNA unless something is done to specifically exclude them. Several methods have been used to selectively discriminate against host DNA in plant-associated samples.

Discriminating primers.

Discriminating primers are designed to selectively bind to bacterial DNA but not plant organellar DNA. Such primers are extremely simple and have been successfully used in Arabidopsis (Bodenhausen et al. 2013), poplar (Beckers et al. 2016), maize (Chelius and Triplett 2001; Wallace et al. 2018; Wasimuddin et al. 2019), apple (Shade et al. 2013), algae (Thomas et al. 2020), various nonmodel plants (Massoni et al. 2020; Vannier et al. 2018), and even for the gut microbes of plant-eating insects (Hanshew et al. 2013). The major drawback of discriminating primers is that they cause bias by also discriminating against legitimate bacterial targets (Thomas et al. 2020). Beckers et al. (2016) compared the efficiency of different discriminating primer pairs in poplar and found that the 799F-1391R pair was most suitable for broad taxonomic coverage, followed by 799F-1193R. (The pair used in the current study—799F-1115R—was not tested but is likely similar to these two because the chloroplast discrimination comes from the 799F primer.)

PCR clamps.

PCR clamps consist of modified oligonucleotides that bind to template DNA during PCR and physically prevent amplification of some products (in this case, organellar DNA). This inhibition usually works by blocking primer sites or preventing the DNA polymerase from physically moving along the template. The most common PCR clamps are peptide nucleic acids (PNAs) (Lundberg et al. 2013; Sakai and Ikenaga 2013), which have been used in a wide variety of plant species (Blaustein et al. 2017; Fitzpatrick et al. 2018; Tkacz et al. 2020). The efficiency of PNA clamps can vary widely across plant species even within similar clades (Fitzpatrick et al. 2018). A related technology uses locked nucleic acids to similar effect (Ikenaga and Sakai 2014; Ikenaga et al. 2016, 2018) but these have had limited use in the plant microbiome field and, therefore, were not included in the current analysis.

Blocking oligos.

Blocking oligos are a recently proposed method that uses a nested PCR to “poison” organellar amplicons so that they cannot be used for library preparation (Agler et al. 2016). Because the nested PCR targets DNA specific to the plant of interest, blocking oligos must be custom designed for each target species. In theory, they should be compatible among closely related species but, to our knowledge, this has not been tested.

To provide guidelines for choosing among these various options, we compared several common methods for DNA extraction and PCR amplification of plant-associated bacteria to identify the potential biases present in them. Although these methods have been used for several years, there has not yet been a head-to-head comparison among them to help researchers make informed choices as to which to use for plants. (Similar comparisons have been done for other species, especially humans) (Fiedorová et al. 2019; Teng et al. 2018). Although we did not include analysis of fungal (and other) communities in this study, we acknowledge that they are also important to the plant microbiome and encourage similar studies to be done for them.

MATERIALS AND METHODS

Sample collection.

Samples for extraction methods (experiment I) consisted of leaf tissue from Arabidopsis thaliana, corn (Zea mays), and soybean (Glycine max) growing at the Center for Applied Genetic Technology (CAGT) at the University of Georgia. All plants were grown under standard conditions (growth chamber or greenhouse), with no special treatments or sterilization. Arabidopsis samples consisted of one mature leaf. For maize and soybean, either 1 or 10 leaf discs were collected using a standard hole punch (approximately 8 mm in diameter); the punch was wiped clean with alcohol between collections to remove any carryover. Soil samples were collected from areas with low (soil 1) or high (soil 2) organic matter. All samples were done in four replicates per extraction method, with soil being mixed and subdivided among the replicates.

Samples for PCR primer tests (experiment II) were handled similarly but came from different locations. Soybean and corn samples came from the Iron Horse Plant Sciences Farm in Watkinsville, GA, and two different soil samples (a high-organic landscaped soil and a low-organic native clay soil) were collected around the CAGT. Defined community DNA (ATCC 20-Strain Even Mix Genomic Material; ATCC MSA-1002) was also tested as a positive control to check for community distortion. A single DNA extraction was used for each sample for primer tests, with four replicates aliquoted during the PCR step.

All samples were frozen at −80 C until extraction. Plant tissue was frozen in tubes preloaded with two 4.5-mm steel beads for later grinding. Soil was weighed into 200- to 250-mg aliquots before DNA extraction.

DNA extraction.

DNA extraction protocols followed the manufacturer’s guidelines, with minor modifications (e.g., adding small amounts of lysis buffer before grinding because it improved grinding efficiency; see below for details). Tissue was ground with two 4.5-mm steel balls unless otherwise stated, and all grinding steps consisted of 90 s of grinding at 850 rpm in a SPEX SamplePrep Geno Grinder 2010. Samples for the DNA extraction tests (experiment I) were extracted according to the indicated protocol. All samples for PCR primer tests (experiment II) were extracted using the Quick-DNA protocol. Details of each specific protocol follow.

Extract-N-Amp Plant PCR Kit.

Extraction buffer (200 μl; Extract-N-Amp Plant PCR Kit, Sigma-Aldrich) was added to the frozen samples, which were then ground as described above. Another 300 μl of extraction buffer was then added (total = 500 μl) and mixed, and the solution was incubated at 95°C for 10 min. Dilution buffer (500 μl) was then added before centrifuging the sample and transferring the supernatant to a clean tube.

PowerSoil DNA Isolation Kit.

Soil samples were handled according to the manufacturer’s protocol (Molecular Biosciences; now Qiagen DNEasy PowerSoil Kit), with two of the samples including 500 μl of lysis buffer in the grinding step to test its effect. (It did not appear to alter extraction efficiency; data not shown.) Initial tests showed that the beads provided with the PowerSoil kit did not do a good job grinding plant leaf disks; therefore, the grinding step of this protocol was modified as follows. Lysis buffer (200 μl) from the PowerBead tubes was added to each sample tube (with two steel beads) and the tissue ground as above. Buffer C1 (25 μl) was added and mixed, and tubes centrifuged to pellet the debris. The supernatant was transferred to a clean tube and 100 μl of buffer C2 was added. The remaining steps were carried out as per the manufacturer’s protocol, and DNA eluted in 100 μl.

DNeasy Plant kit.

Buffer AP1 (100 μl for plants or 200 μl for soils) (Qiagen) was added to the sample tubes prior to grinding, as above. After grinding, additional buffer AP1 was added to total 400 μl. For soil samples only, 4 μl of RNase A was then added. The rest of the protocol followed the manufacturer’s instructions.

Quick-DNA Fungal/Bacterial Miniprep Kit.

BashingBead buffer (200 μl; Zymo Research) was added to both plant and soil samples before grinding as described above. After grinding, 550 μl of BashingBead buffer was added (total = 750 μl) and the tube was vortexed and centrifuged to pellet debris. The rest of the protocol was according to the manufacturer’s instructions.

Design of maize-specific blocking oligos.

Blocking oligos for maize were designed following the procedure of Agler et al. (2016). In brief, BLAST+ v2.2.31+ (Camacho et al. 2009) was used to align 30-bp sections of the maize mitochondria and chloroplast sequences to the Greengenes 97% operational taxonomic unit (OUT) database v13.8 (DeSantis et al. 2006) (distributed as part of the qiime-default-reference v0.1.3 package for python) (Caporaso et al. 2010). Sequences with the fewest Greengenes hits and a low probability of forming hairpins and homodimers (as determined by Primer3) (Rozen et al. 2017) were chosen as new blocking oligos (Supplementary Table S1; primers P10, P16, and P17).

PCR amplification.

All samples had portions of the 16S rRNA gene amplified by PCR for sequencing. The exact protocol depended on the primer set being used. All primer and oligo sequences are shown in Supplementary Table S1. (Note that many of the primers include Illumina Nextera linkers used to prepare sequencing libraries.)

For the DNA extraction tests (experiment I), soil samples were amplified with universal primers (P01 and P02) while the plant samples used the discriminating primer set (P03 and P04). PCR amplification tests (experiment II) used the primer set marked for each sample, and each reaction was run in triplicate and then pooled for sequencing. A portion of each reaction was run on a 1% agarose gel to confirm amplification.

Universal primers.

Universal primer reactions targeted the V4 region of the 16S rRNA gene and consisted of 1 µl of template DNA, 0.5 µl of forward primer P01 (515f; 0.2 µM final concentration), 0.5 µl of reverse primer P02 (806rB; 0.2 µM final concentration), 12.5 µl of 2× Hot Start Taq master mix (New England BioLabs), and 10.5 µl of water (= 25 µl total volume). PCR amplification followed Earth Microbiome Project recommendations (Earth Microbiome Project 2020): 94°C for 3 min; 30 cycles of 94°C for 45 s, 55°C for 60 s, and 72°C for 90 s; and a final 72°C extension for 10 min before holding at 4°C.

Discriminating primers.

Discriminating primer reactions targeted the V5-V7 region of the 16S rDNA gene and consisted of 1 µl of DNA template, 0.5 µl of forward primer P03 (799F; 0.2 µM final concentration), 0.5 µl of reverse primer P04 (1115R; 0.2 µM final concentration), 12.5 µl of 2× Hot Start Taq master mix (New England BioLabs) and 10.5 µl of water (= 25 µl total volume). The PCR program was 95°C for 5 min; 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s; and a final 72°C extension for 5 min before holding at 4°C.

PNA clamps.

PNA clamps bind to the V4 region of the 16S rRNA gene to prevent PCR amplification. PNAs were stored at 100 µM and −20°C until use, during which time a small aliquot was diluted to 5 µM and stored at 4°C for up to 1 week. PCR assays with PNAs were performed according to the manufacturer’s recommendations; namely, 1 µl of template DNA, 1 µl of forward primer P01 (515f; 0.2 µM final concentration), 1 µl of reverse primer P02 (806rB; 0.2 µM final concentration), 1.25 µl of mPNA (0.25 µM final concentration), 1.25 µl of pPNA (0.25 µM final concentration), 12.5 µl of 2× Hot Start Taq master mix (New England BioLabs), and 7 µl of water (= 25 µl total volume). The reaction mixture was placed in the thermocycler for PCR settings prescribed by PNA Bio: 94°C for 3 min; 30 cycles of 95°C for 15 s, 75°C for 10 s, 55°C for 10 s, and 72°C for 60 s; and a final 72°C extension for 10 min before holding at 4°C.

Blocking oligos.

Blocking oligos “poison” undesired PCR amplifications by creating short amplicons that prevent further PCR assays for unwanted sequences. The blocking oligos in this study target either the V3-V4 region of the 16S rDNA gene (BO_3/4) or the V5-V7 region (BO_5/7). The blocking oligo reaction protocol was adapted from Agler et al. (2016) and consisted of two PCR steps. The first reaction uses the blocking oligos to block unwanted transcripts. Enzymatic cleanup is then followed by a second PCR step to further amplify successful targets (which are presumably now much more abundant than the unwanted sequences). Primer names are given at the end because they varied by species and target region; however, the methodology was the same.

First-step amplification consisted of 0.5 µl of template DNA, 0.16 µl of forward primer (0.08 µM final concentration), 0.16 µl of reverse primer (0.08 µM final concentration), 0.5 µl of forward blocking oligo (0.25 µM final concentration), 0.5 µl of reverse blocking oligo (0.25 µM final concentration), 10 µl of 2× Hot Start Taq master mix (New England BioLabs) and 8.18 µl of water (= 20 µl total). The first amplification reaction mixture was placed in the thermocycler for 95°C for 40 s; 10 cycles of 95°C for 35 s, 55°C for 45 s, and 72°C for 15 s; and a final 72°C extension for 2 min before holding at 4°C.

First-round PCR product was subjected to an enzymatic cleanup by adding 1 µl of Antarctic phosphatase (New England Biolabs), 1 µl of exonuclease I (New England Biolabs), and 2.44 µl of 10× Antarctic phosphatase buffer. The reaction was incubated at 37°C for 30 min followed by 85°C for 15 min to inactivate the enzymes. Reactions were centrifuged at 7,000 rpm until a white pellet formed on the wall of the tube, and 10 µl of the supernatant was transferred to a fresh PCR tube.

Second-step PCR amplification consisted 2.0 µl of DNA from step 1, 0.83 µl of forward primer (0.166 µM final concentration), 0.83 µl of reverse primer (0.166 µM final concentration), 25 µl of 2× Hot Start Taq master mix (New England BioLabs), and 21.34 µl of water (= 50 µl total volume), with the amplification program of 95°C for 40 s; 25 cycles of 95°C for 35 s, 55°C for 45 s, and 72°C for 15 s; and a final 72°C extension for 2 min before holding at 4°C. Each reaction was run in triplicate, then pooled for sequencing.

PCR primer varied by species and target region. BO_3/4 refers to primers targeting the V3-V4 region of the 16S rRNA gene, while BO_5/7 targets the V5-V7 region. For all samples, BO_3/4 used amplification primers P05 (B341F) and P06 (B806R) in both PCR steps. The soybean samples used the default chloroplast-blocking oligos from Agler et al. (2016) (P07 [3C30-F] and P08 [c11BV3-R]) in the first PCR, while all other samples used a maize-specific forward blocking oligo (P09) with the default reverse (P08 [c11BV3-R]). For technical reasons, some of the BO_5/7 reactions used amplification primers with linkers in the first PCR (P10 [799F] and P11 [1192R]), while others used ones without linkers (P12 [799F] and P13 [1192R]); all used primers P10 and P11 in the second step. This information is recorded in the sample keys (part of the analysis pipeline at https://github.com/wallacelab/paper-giangacomo-16s-methods) and did not cause any significant differences in the results (data not shown). (Also note that, even though both are called 799F in the literature, primers P10 and P12 have a slightly different sequence from primer P03.) For BO_5/7, mitochondria-blocking oligos for soybean in the first amplification step were from Agler et al. (2016) (P14 [5M30-F] and P15 [cl1BV5-R]), while all other samples used maize-specific ones (P16 and P17).

Library preparation.

For the DNA extraction tests (experiment I), PCR products were purified with a QIAQuick PCR Purification kit (Qiagen) according to the manufacturer’s instructions; those for PCR amplification tests (experiment II) were purified with AGENCOURT AMPure XP magnetic beads. Illumina sample-specific barcodes and sequencing adapters were added in a final PCR consisting of 10 μl of 2× Hot Start Taq master mix (New England BioLabs), 2 μl each of Nextera N7 and S5 barcoding primers (Illumina), 1 μl of purified template, and 5 μl of water (20 μl total) per sample. The PCR cycle was 95°C for 5 min; eight cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s; and a final 72°C extension for 5 min before holding at 4°C.

Barcoded samples were sent to the Georgia Genomics and Bioinformatics Core at the University of Georgia for quantification, normalization, pooling, and sequencing on an Illumina MiSeq. Paired-end 2 × 250 sequencing was performed for DNA extraction tests (experiment I). Because the blocking oligos protocol resulted in longer amplicons, paired-end 2 × 300 was used for the PCR amplification tests (experiment II). All read data are available at the NCBI Sequence Read Archive, Bioproject PRJNA646931.

Bioinformatic processing.

Raw sequencing reads were processed down to amplicon sequence variants (ASVs) using the QIIME 2 pipeline (version 2019-7) (Bolyen et al. 2019). DNA extraction tests were processed separately from PCR amplification tests, and different primer sets were processed separately within each test set, though all followed a common pipeline. All bioinformatic scripts are available at https://github.com/wallacelab/paper-giangacomo-16s-methods.

First, the final 100 (extraction tests) or 150 (primer tests) base pairs of reverse reads were removed with cutadapt (command ‘cutadapt-cut’) (Martin 2011) because initial checks found that the high error rate in these regions interfered with read joining (>90% of reads could be merged after trimming, whereas <50% could be merged without trimming; data not shown). Primers were removed with cutadapt (command ‘qiime cutadapt trim-paired’, with forward and reverse primer sequences supplied), followed by joining with vsearch (‘qiime vsearch join-pairs’) (Rognes et al. 2016). Joined pairs were quality filtered with ‘qiime quality-filter q-score-joined’.

Amplicon sequence variants were called with Deblur (‘qiime deblue denoise-16s’) (Amir et al. 2017), with trim lengths based on the primer set used (250 bp for universal and PNAs, 295 for discriminating, 370 for blocking oligos V5-V7, and 400 for BO V3-V4) and keeping all ASVs with at least two reads across any samples (‘--p-min-reads 2 --p-min-size 1’). Taxonomy was assigned by extracting the targeted regions from the Silva 132 0.99 OTU database (Quast et al. 2013) (‘qiime feature-classifier extract-reads’), training a custom naive Bayes classifier (‘qiime feature-classifier fit-classifiers-naive-bayes’), and using this classifier to assign taxonomy to the ASVs (‘qiime feature-classifier classify-sklearn’) (Bokulich et al. 2018; Pedregosa et al. 2011).

Because different primer sets targeted different parts of the 16S rRNA gene, a phylogenetic tree containing all samples in a test set could not be made directly from the ASVs. Instead, we used BLAST+ v2.2.31 (Camacho et al. 2009) to identify the best hit for each ASV in the Silva 132 0.99 OTU database (Quast et al. 2013), and that hit was used to anchor each ASV in the corresponding Silva phylogenetic tree using the R package ape v5.3 (Paradis and Schliep 2019).

Data unification and filtering.

Treatment comparisons were performed by exporting the Qiime2 artifacts to standard file formats: BIOM (McDonald et al. 2012) for ASV tables, FASTA for sequences, and so on. Phyloseq v1.28.0 (McMurdie and Holmes 2013) was then used to combine ASVs from each primer set into a single unified data table containing ASV counts, sample metadata, assigned taxonomy, and phylogeny, which was then collapsed to the genus level to allow comparison across primer sets. This reduced ASVs to the operational phylogenetic units (OPUs) of Agler et al. (2016). Samples were then filtered to remove any with fewer than 1,000 (experiment I) or 500 (experiment II) total reads, and OPUs were filtered to remove any with fewer than 3 total reads or that were not present in at least two samples.

Statistical analyses.

Organellar contamination of primer sets was tabulated by summing the read counts in each sample with the taxonomic assignment of Chloroplast (order) or Mitochondria (family). Organelles were then removed from the data before performing other analyses.

Alpha diversity of DNA extraction tests was calculated in phyloseq using the plot_richness() function for observed OTUs and Shannon diversity.

Shared and unique OPUs were calculated by converting reads counts to presence/absence data. We then tabulated the OPUs unique to a given sample set or shared with at least one other sample set.

Principal coordinate plots were generated by rarefying samples to 2,000 (experiment I) or 500 (experiment II) reads and calculating the weighted UniFrac distance with rbiom 1.0.0 (Smith 2019) (Phyloseq v1.28.0 appears to have a bug in its UniFrac calculation; see https://github.com/joey711/phyloseq/issues/936). The cmdscale() function in R was then used to convert the distance matrices to principal coordinates.

Taxonomic distortion was calculated by using the phyloseq tax_glom() function to collapse the data at various taxonomic levels (phylum, class, and so on) and the phyloseq_to_deseq2() function to convert it to a DESeq2 (Love et al. 2014) object, with both sample type and treatment as design variables. Because DESeq2 ignores any taxa with 0 reads in any samples, if any collapsed taxon had 0 counts, a pseudocount of 1 was then added to all counts. Data were rarefied to the lowest read count among samples, and the DESeq() function was used to find significantly different taxa. Rarefaction is usually not recommended for DESeq data but we found that, without it, the program tended to find most taxa in a sample changing in a single direction that correlated with the total read depth of that sample. DESeq significance was determined with a Wald test statistic and parametric fit for all levels except domain, which was fit with a mean-dispersion estimate. (This difference was due to a consistent error when trying to fit domain with a parametric dispersion estimate, possibly due to there being only two domains [Bacteria and Archaea] and Archaea having no reads in many samples.) Individual taxa were considered significantly distorted if their adjusted P value was < 0.01 when compared with the reference treatment.

Verification of sample identity was performed using both BLAST (Camacho et al. 2009) and Kraken2 (Wood et al. 2019). BLAST analysis was done by taking the first 1,000 reads of each raw fastq file, BLASTing them against the Silva 132 0.99 OTU database (Quast et al. 2013), and tallying up the species of all hits corresponding to “Chloroplast” or “Mitochondria”. Kraken2 analysis involved using Kraken2 to assign all reads in each sample to a taxon based on a custom database consisting of the NCBI plants database. (Reads that did not match the plant database—such as bacterial reads—were simply “unclassified.”)

Plots were generated with ggplot2 (Wickham 2016) and minor cosmetic adjustments made in Inkscape (Inkscape 2020).

Software packages.

The following software packages were used in this analysis:

• R packages ape v5.3 (Paradis and Schliep 2019), argparse v2.0.1 (Davis 2019), DESeq2 v1.24.0 (Love et al. 2014), dplyr v0.8.3 (Wickham et al. 2019), ggplot2 v3.2.1(Wickham 2016), gridExtra v2.3 (Auguie 2017), igraph v1.2.5 (Csardi and Nepusz 2006), phyloseq v1.28.0 (McMurdie and Holmes 2013), rbiom v1.0.0 (Smith 2019), tidyr v1.0.0 (Wickham and Henry 2019), and vegan v2.5.5 (Oksanen et al. 2019).

• Python packages argparse v1.1, biopython v1.74 (Cock et al. 2009), matplotlib 3.1.0 (Hunter 2007), and primer3-py v0.6.0 (Pruitt and Conway 2020).

• Command-line tools biom v2.1.7 (McDonald et al. 2012), blast+ v2.2.31 (Camacho et al. 2009), Conda v4.8.2, Clustal Omega v1.2.1 (Sievers et al. 2011), cutadapt v2.4 (Martin 2011), libprimer3 v2.5.0 (Rozen et al. 2017), QIIME2 v2019-7 (Bolyen et al. 2019), and vsearch v2.7.0 (Rognes et al. 2016).

RESULTS

Experiment I: DNA extraction.

We tested four different commercial extraction methods covering a range of manufacturers and target tissues (soil, microbial, and plant) to determine how well each captures the microbial community (Fig. 1A). Three of these—MoBio PowerSoil, Qiagen DNeasy Plant, and Zymo Quick DNA—consisted of column-based affinity purification, and one phase-separation method (Extract-N-Amp) was included for comparison. Because the true composition of the community is unknown, the MoBio (now Qiagen) PowerSoil kit recommended by the Earth Microbiome Project was used as the baseline. The recommendation was later updated to the magnetic bead-based version of this kit. (Marotz et al. 2017). Samples consisted of leaves from corn, soybean, and Arabidopsis, and also a soil sample. Because the results of the primer tests (experiment II) were not known at the time, plant samples were amplified with chloroplast-discriminating primers and the soil samples were amplified with Earth Microbiome Project universal primers (see Materials and Methods).

Download as PowerPoint

Because the samples were the same across all extraction methods, we expected that the “best” methods would recover more taxa (= higher alpha diversity). This is admittedly a simplification, because other factors such as sequencing errors and PCR conditions can affect α diversity–though, because all other conditions were the shared across treatments, we expected the effect of such things to be minor, or at least uncorrelated with the treatments. However, no single method followed this pattern for all sample types (Fig. 2). The PowerSoil and Quick-DNA kits recovered the highest α diversity of soil samples, as might be expected because the other two kits had not been designed with soil in mind. However, in plant samples, we found that the DNeasy Plant kit performed best for Arabidopsis, Quick-DNA worked best for maize, and everything performed roughly evenly for soybean. The number of unique taxa found with each method (Supplementary Fig. S1) reinforces this, because the best method for a sample usually has the largest number of unique taxa, while others share most of their taxa with at least one other method.

Download as PowerPoint

Overall community composition was assessed with weighted UniFrac analysis. As expected, samples generally clustered by sample type (Fig. 3A), except that the maize samples extracted with PowerSoil clustered with the Arabidopsis samples. Analysis of the raw data indicates that this was not due to a labeling error (Supplementary Fig. S2); thus, it appears that the PowerSoil kit had a strong effect on the recovered maize microbiome but not on soybean and Arabidopsis. Within each sample type, soil showed the most distinct clustering by extraction method (Fig. 3B). The three plant samples showed less distinct clustering of all methods from each other, although maize leaves showed a very strong separation of PowerSoil extractions from the other three, as mentioned earlier. If we take the PowerSoil kit to be the current standard, none of the other methods consistently overlapped with it, although the Quick-DNA kit came closest in all samples but maize; and, in that case, the PowerSoil kit was likely the one with stronger biases.

Download as PowerPoint

To determine just how different methods compared at recovering different bacterial taxa, we used DESeq2 (Love et al. 2014) to determine which groups at various taxonomic levels were distorted relative to the PowerSoil method (the current community standard) (Fig. 4). This clearly showed that the Quick-DNA method had the least distortion relative to PowerSoil, with most distortion limited to a single genus (Aureimonas) within the phylum Proteobacteria. The DNEasy Plant kit showed significant distortion, with the enrichment of several Proteobacteria clades (especially the genus Sphingomonas) and depletion of Actinobacteria (genera Cutibacterium and Corynebacterium) and the Firmicutes family Staphylococcaceae. A complete table of all significantly distorted clades is provided in Supplementary Table S2.

Download as PowerPoint

Experiment II: Amplification.

We compared three different methods to exclude organellar DNA from 16S rRNA gene amplification of plant samples: discriminating primers, PNA clamps, and blocking oligos (Fig. 1B to D). Each was compared with generic universal primers from the Earth Microbiome Project as the control. The universal primers targeted the V4 region of the 16S rRNA gene, which is also where PNA clamps bind to prevent amplification of plastid sequences. The discriminating primers targeted the V5-V7 region and the blocking oligos targeted either the V3-V4 region (BO_3/4) or the V5-V7 region (BO_5/7) (Fig. 1C). Samples tested were two soil samples, a defined 20-member community (ATCC 20-Strain Even Mix Genomic Material), and two plant leaf samples (maize and soybean). The soil and defined community samples were used to detect bias in the primer sets because they should not contain organelles.

As expected, the universal primers amplified almost pure organellar DNA from plant samples (Fig. 5). Discriminating primers and the blocking oligos targeting the V5-V7 region (BO_5/7) had the lowest organellar sequence (nearly zero), while the PNA clamps and V3-V4 blocking oligos (BO_3/4) had intermediate levels.

Download as PowerPoint

To determine how much each method distorts the view of the underlying community, we used them to amplify bacterial communities from two soil samples and a defined microbial community (Fig. 6). In all three samples, the PNAs and universal primer set clustered together (though not quite on top of each other with soil). Interestingly, the V5-V7 blocking oligos clustered with the universal primers when using a defined community, as reported in the original publication (Agler et al. 2016), but they were strongly distinct in soil samples, indicating that there is still significant distortion occurring. The discriminating primers and V3-V4 blocking oligos always clustered separately, indicating significant distortion.

Download as PowerPoint

To identify specific bacterial clades that are affected by the different primer sets, we determined the significantly different taxa with DESeq2 (Love et al. 2014). The universal primer was set as the reference, and we limited the analysis to just the soil samples (because the universal set has almost no bacterial reads on leaf samples) (Fig. 7). As expected, PNAs showed no distortion relative to the universal primer set. The other three sets showed large numbers of taxa being distorted, including the Archaea and most clades within the Bacteria. All distorted bacterial clades are too numerous to mention (see Supplementary Table S3 for the full list) though, in broad strokes, Proteobacteria are enriched; Verrucomicrobia, Chloroflexi, and Acidobacteria are depleted; and clades within the Actinobacteria can go either way. The distortion pattern is similar between the BO_5/7 and discriminating sets, probably because they share a primer site (799F), though it is surprising that the BO_3/4 set also shows a broadly similar pattern of distortion. A full list of distorted taxa is provided in Supplementary Table S3.

Download as PowerPoint

When looking at the total number of genera recovered by each method (Supplementary Fig. S3) only the V3-V4 blocking oligos showed significantly lower genus counts with the defined community. In soil samples, however, the PNA and universal primer sets recovered dramatically more genera than the other methods. On leaves, the discriminating primers and PNA recovered the most genera; the low number of genera for the universal primer set was likely due to almost all of its reads coming from chloroplasts instead of bacteria.

DISCUSSION

All current microbiome methods are known to result in some degree of bias (Pollock et al. 2018). Because removing all bias is impossible, current best practices are to treat everything within an experiment as identically as possible so that all samples at least suffer the same biases. Nonetheless, researchers generally want to reduce bias as much as possible to get the clearest results they can. With that in mind, our results suggest the following recommendations regarding DNA extraction and amplification of plant-associated samples.

Choice of extraction method (experiment I).

The extraction methods we tested were chosen based on their prevalence in the field or ease of use. Principal coordinate analysis showed a distinct clustering for each method, confirming that each has unique biases (Fig. 4). For soil samples (and by extension, rhizosphere samples), the best extraction methods were unsurprisingly those designed for soil extraction: the PowerSoil and Quick-DNA kits. These kits had the highest recovered α diversity from soil (Fig. 2), showed the largest number of recovered OTUs (most of which were shared between them; Supplementary Fig. S1), were relatively close to each other in overall community composition (Fig. 3), and showed relatively minor taxonomic distortion (Fig. 4).

The choice of extraction method for plant samples was more complex. For Arabidopsis, the DNeasy Plant kit had the greatest α diversity and number of recovered genera whereas, for maize, it was the Quick-DNA kit, and soybean showed no clear winner. Similarly, methods that showed little variability in one sample were highly variable in a different sample, giving no clear guide for choosing. A challenge with plant samples is that the extraction method must not only lyse the microorganisms (including endophytes within the host tissue) but must also remove the phenolic compounds and other inhibitors of downstream processing steps (Pollock et al. 2018). These compounds vary from species to species; thus, it may be that some testing is required to identify the best method for each species. Without that data—or for studies that will span a wide range of species—the best advice is, again, to choose a single method and use it consistently but be aware of host-induced biases that can occur. As for which to choose, from this (admittedly small) dataset, it appears that the Quick-DNA may be the most middle-of-the-road choice, avoiding both the exceptionally good result of the DNEasy Plant kit on Arabidopsis (Fig. 2) and the outlier result of PowerSoil on maize (Fig. 4). It also has consistently modest variability within each sample—rarely the best, but also never the worst—again making it a good middle-of-the-road choice.

When dealing with multiple sample types such as bulk soil, rhizosphere, stem, and leaf, the choice of method should be guided by the sample and goals of the experiment. Any conditions that will be compared should be extracted with the same method. This is necessary to avoid condition-specific biases when comparing, for example, bulk soil to rhizosphere or either of those to stem tissue to trace the origin of plant endophytes. If, however, there is no reason to compare across compartments, then each compartment can be extracted with the method that works best for it. This will recover a better snapshot of the community in each compartment but at the cost of not being able to compare between them. It also adds more complexity to the experiment that could result in an error, such as by using the wrong method on some of the samples.

Choice of 16S rRNA gene amplification method (experiment II).

The goal of choosing an amplification protocol is to minimize bias while also minimizing the amount of host DNA in the sample. Host contamination is not generally a problem for soil samples but becomes crucial for samples composed mostly of host tissue (roots, stems, leaves, and so on).

Of the methods we tested, the blocking oligos (specifically, BO_5/7) and discriminating primers did the best job at excluding DNA from organelles, excluding almost all organellar DNA (Fig. 5). BO_3/4 was more modest, at approximately 25% organellar reads. The PNA clamps also showed only modest control (30 to 75% organelle DNA), though it should be noted that we tested only a single concentration of PNAs, and that leaf tissue is especially high in plastid sequence. Higher concentrations might be able to exclude more plastid sequence (at the cost of using more PNA reagent per reaction). In other experiments in our lab, we have also found these concentrations to be sufficient for nonphotosynthetic tissue, such as roots and stems (data not shown).

Comparing amplification of nonplant samples (two soils and a defined community) revealed that only the PNAs maintained a community composition comparable with the universal controls (Fig. 6). The PNA method was also consistently second-best in terms of number of bacterial genera recovered (Supplementary Fig. S3), coming behind the universal primers for soil samples and the discriminating primers for plant samples. Both sets of blocking oligos performed consistently more poorly than the other methods. We also found them much more difficult to work with, with a more complex protocol and a much higher failure rate than the other methods (personal observation). This implies that using blocking oligos effectively in a new species will require significant optimization.

Given these results, we recommend the use of either PNA clamps or plastid-discriminating primers. The PNA clamps are preferred because they minimally distort the underlying community while still increasing target sequences approximately 20-fold; higher concentrations of PNA may be able to improve this further (but at increased reagent cost). Discriminating primers are slightly easier and much less expensive but do result in significant taxonomic distortion. Blocking oligos do not appear to be a good option for researchers seeking out-of-the-box solutions, as we did here. They are very species specific, have a more complex protocol, and appear to distort the taxonomic profile (Fig. 6). The V3-V4 set used here and recommended in the original article also requires longer read lengths than most other methods (paired-end 300 instead of paired-end 250), which also increases costs and limits which Illumina sequencing platforms can be used. The significantly lower cost of blocking oligos and the ability to design them for multiple loci (Agler et al. 2016) means they may still be worthwhile for labs planning to do many experiments in the same species. However, labs seeking off-the-shelf solutions for a few experiments are probably better off with another method.

Implications for plant-associated microbiome work.

Plant-associated microbiome analysis presents unique challenges due to the abundance of host-associated DNA, especially from plant organelles. Probably for this reason, the first release of the Earth Microbiome Project data included almost no plant-associated samples (Thompson et al. 2017). However, plant-associated microbes play crucial roles in global ecosystems and agriculture. Understanding these interactions will be necessary to fully harness them during the 21st century. Although several of the methods compared in this article help probe the plant microbiome, better methods are still needed. Critically, to our knowledge, there are no easy methods to reliably separate plant host DNA from microbes, a necessary step for whole-metagenome profiling of the plant microbiome (especially the plant endosphere). Progress is being made in this direction (Carrión et al. 2019; Nobori et al. 2018), and we encourage further development of these methods to help unveil additional aspects of plant–microbe interaction.

ACKNOWLEDGMENTS

We thank M. Agler for help troubleshooting the blocking oligos protocols.