Methods in Ecology and Evolution

Accepted Article
DR ANDREW DOPHEIDE (Orcid ID : 0000-0003-1554-9832)

Article type : Research Article

Handling editor: Professor Michael Bunce

Handling editor: Professor Michael Bunce

Impacts of DNA extraction and PCR on DNA metabarcoding

estimates of soil biodiversity

Andrew Dopheide1,2,3*, Dong Xie4, Thomas R. Buckley1,3, Alexei J. Drummond4, and Richard D.

Newcomb1,2
1School of Biological Sciences, The University of Auckland, New Zealand

2The New Zealand Institute for Plant & Food Research, Auckland, New Zealand

3Manaaki Whenua - Landcare Research, Auckland, New Zealand

4Center for Computational Evolution, The University of Auckland, New Zealand

*Email: [email protected]. Current address: Manaaki Whenua - Landcare Research, Auckland, New

Zealand.

Running headline: DNA metabarcoding and soil biodiversity estimates

This article has been accepted for publication and undergone full peer review but has not been
through the copyediting, typesetting, pagination and proofreading process, which may lead to
differences between this version and the Version of Record. Please cite this article as doi:
10.1111/2041-210X.13086
This article is protected by copyright. All rights reserved.
 Abstract
Accepted Article
1. Soils are ubiquitous and important components of terrestrial ecosystems, richly populated with

diverse organisms with important ecological roles. DNA metabarcoding is a promising tool for

efficient and taxonomically comprehensive analyses of soil biodiversity, but the outcomes of these

analyses are likely affected by basic methodological factors such as sampling and laboratory

protocols.

2. We investigated the impacts of DNA extraction methodology and multiple PCRs on DNA

metabarcoding biodiversity estimates for multiple taxonomic groups in soil. We applied four DNA

extraction methods with different size constraints in parallel to fixed volumes of soil from a pair of

forest sites, followed by sets of ten individually indexed PCRs and Illumina sequencing of prokaryote

(16S) and eukaryote (18S, fungal 26S, and metazoan COI) barcodes.

3. Methods utilising larger soil volumes resulted in higher biodiversity estimates for arthropods but not

for prokaryotes or micro-eukaryotes, and improved spatial discrimination of metazoan communities

but not prokaryotes. Each DNA extraction method resulted in a biased community composition, and

these biases were consistent across sites and amplicons. Some 36 to 41 % of OTUs were shared

between all four DNA extraction methods and OTUs from many phyla had contrasting abundances

from different methods. From 0.2 to 19 % of OTUs, accounting for 9.4 to 97 % of sequence reads,

were amplified in all ten PCRs from each DNA extraction method and sample. OTUs detected in only a

single PCR accounted for 26 to 41 % of OTUs but only 0.3 to 2.8 % of the sequence reads from each

extraction method and sample. The inferred species richness increased strongly with the number of

PCRs carried out.

4. We conclude that a consistent DNA extraction method should be used for comparisons between

samples, confirm that larger samples of soil (≥ 15 g) should be used for analyses of metazoan

biodiversity than for microbes, and that the number of PCR replicates undertaken strongly affects

species richness estimates.

Key-words: DNA extraction, eukaryotes, metabarcoding, metazoans, PCR replicates, prokaryotes, soil

biodiversity

This article is protected by copyright. All rights reserved.

 Introduction
Accepted Article
DNA metabarcoding techniques have enabled the efficient characterisation and monitoring of prokaryote and

eukaryote communities in terrestrial ecosystems (e.g. Buée et al., 2009; Caporaso et al., 2011; Ji et al., 2013;

Drummond et al., 2015). Soil is a promising substrate for comprehensive DNA metabarcoding analyses of

terrestrial biodiversity, due to the great abundance and diversity of organisms and biological debris that are

present in soils (Giller, 1996; Torsvik & Øvreås, 2002; Buée et al., 2009). Furthermore, soil communities are

key contributors to ecological processes such as decomposition and nutrient cycling, yet knowledge of these

communities lags behind that of above-ground biodiversity (Bardgett & Van Der Putten, 2014; Decaëns,

2010).

DNA metabarcoding analyses involve a series of methodological decisions that are likely to affect measures of

biodiversity (Alberdi et al., 2017; Dickie et al., 2018; Murray et al., 2015). The effects of DNA extraction

methods on DNA yield and detected microbial community composition have been examined for substrates

including oral and rumen samples (Henderson et al., 2013; Vesty et al., 2017), permafrost (Vishnivetskaya et

al., 2014), and soils (Delmont et al., 2011; Kang & Mills, 2006; Wagner et al., 2015). In soils, the detection of

prokaryote communities is affected by site-specific factors such as spatial heterogeneity and the propensity of

cells and molecules to adsorb to soil particles (Lombard et al., 2011), sample size (Kang & Mills, 2006) and

extraction method (Delmont et al., 2011). Eukaryotes, particularly metazoans, are larger and less abundant

but more heterogeneous in morphology and size than prokaryotes, suggesting that different sampling and

DNA extraction strategies are necessary to characterise the biodiversity of these organisms in soil. Taberlet et

al. (2012) proposed a method for recovery of extracellular DNA from large volumes of soil in order to reduce

sample heterogeneity, resulting in lower DNA yield and quality but congruent community composition

compared to a conventional DNA extraction kit (Zinger et al., 2016). Otherwise, the impacts of DNA extraction

method on biodiversity estimates for eukaryote groups in soil have rarely been assessed.

PCR amplification is susceptible to various biases, artifacts, and errors (Acinas et al., 2005; Elbrecht & Leese,

2015; Murray et al., 2015; Polz & Cavanaugh, 1998), the occurrences of which are influenced by factors

including the DNA template diversity (Fonseca et al., 2012), and thermal cycling parameters (Qiu et al., 2001;

Stevens et al., 2013). Stochastic PCR variability may also affect the outcomes of metabarcoding-based

This article is protected by copyright. All rights reserved.

 biodiversity analyses (Leray & Knowlton, 2017; Tedersoo et al., 2010; Zhou et al., 2011), particularly for
Accepted Article
low-abundance template molecules. Replication of PCRs is widely included in metabarcoding analyses to

counter such variability, but the number of replicates that are included is usually arbitrary. Ficetola et al.

(2015) and Piggott (2016) investigated how many replicates were required to confidently detect a target

species, which ranged from two to six (Piggott, 2016), or eight (Ficetola et al., 2015) depending on the

methods used and the detection probability of the target. The effects of PCR replication on multi-taxon

metabarcoding biodiversity estimates have otherwise received little attention. For example, it is unclear how

the number PCRs undertaken affects estimates of soil biodiversity, and whether this varies among different

taxonomic groups.

In this study we applied four different DNA extraction methods with differing size limits in parallel to fixed

volumes of soil collected from a pair of spatially separated forest sites, followed by Illumina sequencing of sets

of ten individually-indexed PCRs for prokaryote 16S, eukaryote 18S, fungal 26S, and metazoan cytochrome c

oxidase subunit I (COI) barcode genes. We investigated whether the DNA extraction method affects the range

of organisms detected or biodiversity estimates for prokaryotes and eukaryote groups, and whether spatially

separated soil communities can be discriminated based on the results of different extraction methods and

taxa. In addition, we examined the variability of PCR amplifications and how the number of PCRs undertaken

affects estimates of prokaryote and eukaryote soil biodiversity.

Methods

Sample collection

Soil samples were collected in the forested Waitakere Ranges Regional Park, located west of Auckland in the

North Island of New Zealand (latitude -36.906, longitude 174.520). The sample site is a 20 x 20 m grid of 16

subplots (each 5 x 5 m), identified by letters A through P. Approximately 2.0 kg of soil was collected from each

of two diagonally opposite corner subplots (D and P, approximately 21 metres apart; Figure S1), by excavating

soil and overlying leaf litter at semi-random locations to a depth of 10 cm. The soil material was sealed in

plastic bags and stored at -80 °C until DNA extraction was undertaken.

This article is protected by copyright. All rights reserved.

 DNA extraction, PCR and sequencing
Accepted Article
Each soil sample was homogenised before removal of four 15 g subsamples, for DNA extractions using four

different methods followed by a fixed number of PCRs (Figure S1). The extraction methods included two

widely used soil DNA extraction kits with contrasting upper size limits, and two methods that can be applied

to samples of arbitrary size that were expected to be more effective for metazoan biodiversity analyses.

In the first DNA extraction method (1.5 g direct), DNA was extracted from ten 1.5 g soil aliquots using the

MoBio Powersoil RNA extraction kit with DNA elution accessory kit (MoBio Laboratories, Carlsbad, CA, USA).

These DNA extracts were further processed through a Nucleospin Soil extraction kit (Macherey-Nagel GMBH

& CO. KG, Dueren, Germany), omitting lysis steps, to remove apparent PCR inhibitors. In the second method

(7.5 g direct), DNA was extracted from two 7.5 g soil aliquots using the MoBio Powermax DNA extraction kit

(MoBio Laboratories, Carlsbad, CA, USA). In the third method (indirect extraction), biological material was

separated from 15 g soil using a modified sugar centrifugation protocol (Freckman & Virginia, 1993;

Velasco-Castrillón et al., 2014), before DNA extraction from the separated material using a Nucleospin Soil

DNA extraction kit (Macherey-Nagel GMBH & CO. KG, Dueren, Germany). In the fourth method (phosphate

buffer extraction), 15 g soil was mixed with saturated phosphate buffer (0.12 M¸ pH 8) for 20 minutes before

DNA extraction from a 100 µl aliquot of buffer liquid using a Nucleospin Soil DNA extraction kit

(Macherey-Nagel GMBH & CO. KG, Dueren, Germany), excluding the lysis steps (Taberlet et al., 2012).

Four barcode regions covering the major taxonomic groups found in soil were amplified by PCR from each

extract: an average 570 b.p. (sd = 7) of the prokaryote 16S V4-V6 region using primers 530F and 1100R

(Acosta-Martínez et al., 2008); 338 b.p. (sd = 37) of the eukaryote 18S V7-V8 region using primers #3-F and

#5-RC (Machida & Knowlton, 2012); 578 b.p. (sd = 37) of the fungal 26S D1-D2 region using primers NL1 and

NL4 (Kurtzman & Robnett, 1998); and a 313 b.p. fragment of the metazoan COI gene, using primers mICOIintF

and jgHCO2198 (Leray et al., 2013) (Table S1). One PCR was carried out for each of the ten 1.5 g direct

extracts, five PCRs were carried out on each of the two 7.5 g direct extracts, and ten PCRs were carried out for

each of the 15 g indirect and 15 g phosphate buffer extracts, giving a total of ten PCRs per gene from each 15 g

soil subsample, for two subplots (80 PCRs in total for each gene). The sets of ten PCRs from the 15 g indirect

This article is protected by copyright. All rights reserved.

 and 15 g phosphate buffer extracts are PCR replicates, as are the sets of five PCRs from 7.5 g direct extracts,
Accepted Article
whereas the PCRs from 1.5 g direct extracts represent replicate subsamples/extractions.

PCRs were carried out in 25 µl reaction volumes containing 0.5 units of KAPA 2G Robust Hotstart polymerase

(Kapa

Biosystems Inc, Boston, Massachusetts), MgCl2 at a concentration of 1.5 mmol l-1, dNTPs at a concentration of

0.2 mmol l-1, forward and reverse primers at a concentration of 0.2 µmol l-1 each, and 20-25 ng of template

DNA.

The temperature cycle used for PCRs was (a) 94 °C for 3 minutes; (b) 25 cycles (16S), 27 cycles (18S and 26S),

or 30

cycles (COI) of 94 °C for 30 seconds, 60 °C (16S), 58 °C (18S), 52 °C (26S), or 48 °C (COI) for 30 seconds, and

72 °C

for 45 seconds; (c) 72 °C for 3 minutes. All PCR mixes were prepared in a dedicated laminar flow hood with a

dedicated set of micro-pipettes. Negative control reactions containing nuclease-free water instead of DNA

template

were carried out in parallel with all sets of PCRs.

PCR amplicons were cleaned up using the Ampure XP magnetic bead system (Beckmann-Coulter, MA, USA)

and pooled in sets of four (one each of 16S, 18S, 26S, and COI amplicons), resulting in a total of 80 amplicon

pools. Dual barcode indices and sequencing adaptors were attached to each amplicon pool using an Illumina

Nextera XT Index kit, followed by a further Ampure XP cleanup step, quantification and normalization of DNA

concentration, and combining of pools for sequencing. The combined pool was sequenced by New Zealand

Genomics Ltd. in a 2 x 300 b.p. Illumina MiSeq run.

Sequence processing and bioinformatics

Paired sequence reads for each amplicon pool were merged using PEAR paired-end read merger version 0.9.3

(Zhang et al., 2014), then separated by primer, trimmed of primers, and relabelled by amplicon pool using

custom python scripts. Few of the paired reads from the prokaryote 16S and fungal 26S amplicons

successfully merged, as the length of these amplicons (including primers) tended to exceed the maximum read

This article is protected by copyright. All rights reserved.

 length of the MiSeq system, and therefore only the (higher-quality) R1 sequences from these amplicons were
Accepted Article
used in analyses.

The merged sequences (18S and COI) were not length trimmed whereas the 16S and 26S R1 sequences were

uniformly trimmed to 250 b.p. to remove the lowest-quality 3ʹ regions, followed by quality filtering, initially at

different stringency levels (removal of sequences with >2.5, > 1.0, or > 0.5 maximum expected errors) to

examine the effects of this factor on alpha diversity estimates (Table S2). For subsequent analyses, any

sequences with maximum expected errors > 1.0 were removed. Any duplicates among the quality-filtered

sequences were removed, followed by clustering of the non-duplicates into OTUs at a 97 % identity threshold

(with singleton sequences either included or excluded), with simultaneous chimera filtering, after which the

quality-filtered sequences were matched to the OTUs at a 97 % identity threshold, using UPARSE and other

functions within USEARCH v8.0.1623 (Edgar, 2013). The taxonomic identity of OTUs was determined using

BLAST to compare sequences with the NCBI Genbank nr nucleotide database. BLAST results were parsed

using the lowest-common-ancestor (LCA) algorithm in MEGAN Metagenome Analyser (Huson et al., 2011),

resulting in a consensus taxonomic assignment for each OTU.

Biodiversity analyses

To compare the biodiversity detected among different samples and DNA extraction methods based on

different taxa, the pool of sequences per PCR was randomly subsampled to the lowest number of sequence

reads per PCR within the results from each amplicon (after exclusion of any PCRs with < 25 % of the mean

sequence abundance per PCR), after which richness and effective diversity metrics (Jost, 2006) were

calculated for taxonomic groups based on the subsampled PCRs from each sample and extraction method

using the R package vegetarian (Charney & Record, 2012). The subsampling and diversity calculations were

repeated ten times for each amplicon.

The overlap between DNA extraction methods was investigated by determining the proportions of OTUs with

minimum abundances of 1 to 10 sequence reads that were detected in the results of multiple extraction

methods. To identify OTUs with differing abundances between samples and DNA extraction methods, OTU

abundances were normalised according to a negative binomial model, as recommended by McMurdie and

Holmes (2014) using the R packages DESeq2 (Love et al., 2014) and phyloseq (McMurdie & Holmes, 2013).

This article is protected by copyright. All rights reserved.

 Log 2-fold differences in normalised OTU abundances between samples and extraction methods were
Accepted Article
identified by negative binomial Wald tests in DESeq2, with correction for multiple inferences using the

Benjamini-Hochberg method for control of the False Discovery Rate (Benjamini & Hochberg, 1995).

Non-metric multidimensional scaling (MDS) ordinations and other multivariate comparisons of community

composition between samples, extraction methods, PCRs and taxonomic groups were carried out using the R

package vegan (Oksanen et al., 2017).

To analyse the effects of multiple PCRs on biodiversity estimates, 25 random combinations of results from one

to ten PCRs were selected for each amplicon and extraction method, followed by calculation of diversity

estimates from each combination of PCRs, both with and without subsampling to an equal number of

sequence reads for every amplicon and extraction method. To investigate whether exclusion of low

abundance OTUs affected these results, these calculations were repeated with OTUs with total abundances < 5

excluded, and with OTUs occurring in < 2 of each set of ten PCRs excluded. In addition, to investigate the effect

of sequencing depth on biodiversity estimates, the total sequence pool from each amplicon and method was

randomly subsampled to ten-percentile intervals from 10 to 100 %, followed by calculation of diversity

estimates. All biodiversity analyses were based on OTU clustering with singleton sequences excluded, unless

noted otherwise. Figures were generated using ggplot2 (Wickham, 2009). All analyses were carried out in

RStudio (RStudio team, 2015) with R version 3.1.1 (R Core Team, 2016).

Results

Overall sequencing results

We obtained a total of 9.8 million R1 and 8.8 million R2 sequence reads, with 1.3 million (26S) to 3.5 million

(18S) sequences per amplicon (Table S2). The 16S dataset was dominated by bacteria, including 16 bacterial

phyla, with Proteobacteria most abundant followed by Actinobacteria, Acidobacteria, and Bacteroidetes

(Table S3, Figure S2). The 18S dataset contained the broadest range of eukaryote taxa, including six metazoan

phyla, six fungal phyla, and 12 protist groups. The 26S dataset included a similar range of fungal and

metazoan phyla but fewer protist groups than the 18S dataset. The COI dataset included OTUs from 11

metazoan phyla, two fungal phyla, and 9 protist groups.

This article is protected by copyright. All rights reserved.

 Biodiversity differences between DNA extraction methods
Accepted Article
We observed differences in biodiversity estimates for eukaryote taxa detected by different extraction

methods, related to organism size (Figures 1, S3-S6). In the 18S dataset, the mean species richness and

effective alpha diversity estimates for single-celled eukaryote groups such as Amoebozoa, Cercozoa, and (to a

lesser extent) Ciliophora were highest in the results of the the smallest extractions (1.5 g direct), and effective

gamma diversity estimates for these groups were highest in the results of ten 1.5 g direct extractions (Figure

1). For example, Cercozoan species richness in 1.5 g direct extract results was 73 and 81 per subplot,

compared to 49-58 for the other methods.

In contrast, mean species richness and effective alpha diversity estimates for 18S Arthropoda OTUs were

highest in the results of 15 g phosphate buffer extractions (e.g. richness of 80 and 94 per subplot compared to

18-64 for the other methods). Similarly, alpha diversity estimates for 18S Nematoda OTUs were highest in the

results of 15 g phosphate buffer and 7.5 g direct extractions. The estimated gamma diversity of 18S

Arthropoda OTUs from 15 g phosphate buffer extractions was similar to estimates from ten 1.5 g direct

extractions or two 7.5 g direct extractions, and lowest from 15 g indirect extractions. Effective beta diversity

estimates were highest in the results of 1.5 g direct extractions for 18S Arthropoda, Nematoda, and

Platyhelminthes OTUs (as well as Basidiomycota), but more consistent among the results of different

extraction methods for micro-eukaryote groups. Similar patterns were observed for Arthropoda biodiversity

estimates from the COI dataset (Figure S6).

Alpha diversity estimates for 16S bacterial phyla differed little among the different extraction methods, except

for Proteobacteria, for which the highest mean effective alpha diversity was detected from 1.5 g direct

extractions (81 and 95 per subplot) and the lowest from phosphate buffer extractions (16 and 27 per subplot)

(Figure S3). Alpha and gamma diversity estimates for 26S Ascomycota OTUs varied more between subplots

than among extraction methods (except for 1.5 direct extractions), whereas alpha diversity of 26S

Basidiomycota OTUs was more consistent among samples and methods (Figure S5).

OTU differences between DNA extraction methods and samples

Based on OTU clustering with singletons excluded, a large minority of OTUs were present in the results of all

four DNA extraction methods (41 % of 16S, 40 % of 18S, 37 % of 26S, and 36 % of COI OTUs) (Figure 2). From

This article is protected by copyright. All rights reserved.

 12 % (18S) to 21 % (COI) of OTUs were present in the results of only one extraction method, while the
Accepted Article
remainder of OTUs occurred in intermediate combinations of DNA extraction methods. Further exclusion of

low abundance OTUs resulted in increasing proportions of OTUs being shared by all methods. For example, 42

% (COI) to 60 % (16S) of OTUs with abundances >= 5 were present in the results of all methods, as were 54 %

(COI) to 78 % (16S) of OTUs with abundances >= 10. In contrast, if singletons were included in OTU clustering,

47 % (COI) to 75 % (16S) of the resulting OTUs were present in the results of only a single DNA extraction

method (Figure S7).

Comparisons of normalised OTU abundances using DESeq2 identified OTUs from a very wide range of

taxonomic groups with significantly differing abundances between methods and samples (Table S4; Figures

S8-S11). Overall, in each dataset there were more OTUs with significantly greater abundance in the results of

1.5 g direct extracts compared to the other extraction methods than vice versa (Table S4). For example, in the

18S results, there were 1027, 1256, and 1095 eukaryote OTUs with greater abundance in 1.5 g direct extract

results than in 7.5 g direct, 15 g indirect, or 15 g phosphate buffer extract results respectively, whereas the

opposite was true of 370, 507 and 307 eukaryote OTUs, respectively. Similar disparities were observed for the

other amplicon datasets.

Variability and spatial discrimination of samples

Non-metric multidimensional scaling (MDS) ordinations indicated that OTU assemblages resulting from PCRs

from the same DNA extraction method and sample were in most cases more similar than those resulting from

different DNA extraction methods and samples (Figures 3, S12-S15). Sets of replicate PCRs from a single DNA

extract were highly consistent. For example, data points representing PCR replicates from each of the 15 g

phosphate buffer and indirect extraction methods were in most cases tightly clustered suggesting similar

composition. Results from the 7.5 g direct extractions from each sample were usually grouped into two

clusters corresponding to sets of five replicate PCRs from each of two extracts, whereas sets of PCRs from

individual 1.5 g direct extractions showed the most heterogeneity.

Strikingly consistent patterns of community structure related to DNA extraction method were evident

between amplicons, taxonomic subsets of each amplicon dataset, and subplots (Figures 3, S12-S15). For

example, ordinations based on all 18S, 26S and COI OTUs were very similar to each other, as were ordinations

This article is protected by copyright. All rights reserved.

 based on the fungal components of these datasets, with consistent patterns across both subplots. 18S and COI
Accepted Article
metazoan OTU ordinations also showed very consistent patterns between the two subplots, as did the 16S

prokaryote ordination, albeit with a different pattern of similarity among extraction methods compared to the

eukaryote datasets (Figure 3). Results from the two different subplots usually formed two non-overlapping

groups regardless of extraction method or taxonomic group, with any overlap generally limited to a subset of

the 1.5 g direct extraction results.

Results from 1.5 g direct extractions had the most overlap of within- and between-subplot pairwise distances,

for all amplicons and taxonomic groups (Figure 4). The greatest separation of within- and between-subplot

pairwise distances was observed for the 15 g indirect or 15 g phosphate buffer extract results, especially for

the 18S and COI metazoan OTUs.

PERMANOVA analyses provided evidence of significant differences between subplots and extraction methods

for all amplicons, along with interactions between subplots and methods (except for the 16S dataset when

Jaccard distance was used; Table 1). Dispersion differences contributed to the contrasts among DNA

extraction methods (mainly due to the greater heterogeneity of 1.5 g direct extract results), but not to

differences between subplots (data not shown).

Variability of PCRs

Considerable variation in OTU abundance, from one or two to thousands of sequences per OTU, was observed

in all amplicon datasets. OTUs that occurred only in the results of a single PCR accounted for 26 to 41 %

(mean 31 %) of the total OTUs detected in each set of ten PCRs (Figure 5). However, these OTUs represented

only from 0.3 to 2.8 % (mean 1.2 %) of the total sequence reads in each dataset. OTUs that occurred in all ten

PCRs, in contrast, accounted for 0.2 to 19 % (mean 11 %) of the total OTUs, but 9.4 to 97 % (mean 77 %) of

the total sequence reads from each set of ten PCRs. These consistently-detected OTUs accounted for the

smallest proportions of sequences among sets of ten 1.5 direct extractions and PCRs, particularly for the 26S

results (25 % and 44 % per subplot compared to 83 to 92 % for the other methods) and COI results (9.4 %

and 28 % per subplot compared to 59 to 89 % for the other methods). The proportions of OTUs that occurred

in at least half of each set of PCRs ranged between 16 % and 35 % for the 1.5 direct extractions to between 32

This article is protected by copyright. All rights reserved.

 % and 42 % for the 15 g indirect extractions. In most cases, a higher proportion of OTUs occurred in all ten
Accepted Article
PCRs than in intermediate numbers of PCRs.

Within each set of ten PCRs from a single DNA extraction method, the ten most abundant OTUs were usually

the same and ranked in the same abundance order in each PCR, even when substantial variation in the total

abundance of sequence reads was observed between PCRs (Figures S17-S20). These OTUs were most

consistent among replicate PCRs from 15 g indirect and phosphate buffer extractions, and least consistent

among sets of PCRs from 1.5 g direct extractions.

Effects of multiple PCRs on biodiversity estimates

Species richness estimates increased as increasing numbers of PCRs were pooled together (Figures 6, S21).

The species richness resulting from ten pooled PCRs was typically some two- to four-fold higher than the

richness resulting from individual PCRs. Species accumulation curves flattened as the number of PCRs

increased, suggesting that extrapolated richness might plateau in the range of 10 to 20 pooled PCRs. Effective

alpha diversity increased at lower rates from a lower base, levelling off at intermediate numbers of pooled

PCRs, with the largest increases (some two- to three-fold) observed in pooled sets of PCRs from 1.5 g direct

DNA extractions. When an equal number of sequence reads were drawn from increasing numbers of PCRs,

there were only slight increases in richness for some datasets as the number of PCRs increased. Excluding

OTUs with total abundances < 5 resulted in slightly reduced maximum biodiversity estimates, but otherwise

very similar curves as those resulting from all OTUs (Figure S22), whereas excluding OTUs that occurred in

fewer than two of each set of ten PCRs resulted in curves plateauing in the range of 7-10 PCRs (Figure S23).

Biodiversity estimates for incremental subsamples of the total sequence pool from each set of ten PCRs

resulted in very similar trends to the estimates based on combinations of different numbers of PCRs (Figure

S24).

This article is protected by copyright. All rights reserved.

 Discussion
Accepted Article
Our results show that sample size in combination with DNA extraction method and the number of PCRs have

important effects on biodiversity estimates for different taxonomic groups within soil. The alpha diversity

estimates observed for different methods and taxonomic groups confirm that extraction of DNA from larger

amounts of soil tends to result in higher estimates of metazoan biodiversity (Taberlet et al., 2012). Although

the highest metazoan gamma diversity estimates tended to result from sets of ten 1.5 g direct extractions, the

higher estimates of arthropod alpha diversity in individual PCRs from 7.5 g direct and 15 g phosphate buffer

extractions suggest that ten separate extractions using either of the latter methods would result in more

comprehensive samples of metazoan biodiversity than the former. Larger samples also resulted in the most

discrimination of within- and between-subplot distances for eukaryotes (especially metazoans), whereas

overlapping distances based on 1.5 g soil DNA extracts suggested that this is not ideal for comparing

metazoan communities among habitats. Of course, this may be due to each set of 1.5 g direct extracts and

PCRs representing ten small subsamples of soil, whereas the results from the larger extraction methods

represent replicate PCRs from one (or two) large samples, resulting in greater heterogeneity among the

former than the latter. The 15 g indirect extractions tended to produce lower estimates of metazoan

biodiversity than the 15 g phosphate buffer extractions, despite the former method being targeted at

microfauna (Velasco-Castrillón et al., 2014). This suggests that the sieving and centrifugation steps involved in

this method may have resulted in the exclusion of metazoans that are represented as body fragments or free

DNA molecules, rather than whole organisms.

In contrast to the metazoan biodiversity estimates, our results suggest that DNA extraction size is a less

important consideration for analyses of prokaryote or micro-eukaryote biodiversity. Indeed, the smallest DNA

extraction sizes resulted in the highest biodiversity estimates for prokaryotes and some micro-eukaryote

groups, whereas for other micro-eukaryote phyla there was little difference among the methods. This is

consistent with previous studies, which have shown that soil volumes < 0.5 g can provide representative

samples of bacterial communities (Kang & Mills, 2006; Ranjard et al., 2003), whereas volumes ≥ 1 g are

necessary for representative sampling of fungal communities (Ranjard et al., 2003).

This article is protected by copyright. All rights reserved.

 These results are consistent with the suggestion by Taberlet et al. (2012) that larger samples (i.e. larger DNA
Accepted Article
extraction sizes) are necessary to obtain representative biodiversity estimates for metazoans in soil than for

micro-organisms. This can be explained by differences in the sizes, abundances and distributions of micro-

and macro-organisms, with metazoans being larger, less abundant and more sparsely distributed than

prokaryotes or micro-eukaryotes. Of course, larger samples than 15 g can be processed using the indirect and

phosphate buffer methods, and the optimal sample size for analysis of metazoans in soil may well lie beyond

the range considered in this study.

Although we observed no obvious differences in the taxonomic range of phyla or classes detected by each

DNA extraction method, the strikingly consistent patterns of community structure between the results of

different DNA extraction methods applied in different subplots suggest that each method results in the

detection of a biased assemblage of OTUs. These results show that a consistent DNA extraction methodology

should be used if biodiversity is to be compared among samples and ecosystems, otherwise extraction method

biases may cause misleading comparisons. Application of all methods together would likely result in detection

of the widest range of organisms, as is the case for analyses of soil microbes (Delmont et al., 2011).

Comparisons of OTU abundances between the results of different DNA extraction methods suggested that

many more OTUs were significantly more abundant in the results of 1.5 g direct extractions compared to each

other method than vice versa, for all amplicons. OTU abundances are influenced by PCR amplification

efficiency, which is in turn affected by factors including inhibitory molecules (Miller et al., 1999) and template

complexity (Fonseca et al., 2012), which are likely to vary among different extraction methods. In this case,

most of the 1.5 g direct DNA extracts initially failed to amplify, but were subjected to an additional inhibitor

removal step that resulted in greatly improved PCR outcomes. The 1.5 g direct DNA extraction method may

also have resulted in less diverse pools of DNA template molecules than the other methods, due to the sample

volume. Together, these factors may have caused improved amplification of template molecules from the 1.5 g

direct extractions, resulting in higher OTU abundances compared to the other methods.

The barcode identification of individual amplicon pools enabled examination of amplification consistency

among sets of PCRs. The overall multivariate similarity of PCRs from the same DNA extract is reassuring,

suggesting consistent amplification of molecules in each case. The consistency of the most abundant OTUs

within each set of PCRs is particularly striking, showing that certain templates within each extract are always

This article is protected by copyright. All rights reserved.

 amplified to high abundance, even when the overall efficiency of PCR amplification differs (see Porazinska et
Accepted Article
al., 2009; Amend et al., 2010; Egge et al., 2013; Elbrecht & Leese, 2015). However, there was increasing

variability of less abundant OTUs between PCRs, as demonstrated by an inverse relationship between the

total number of OTUs, and their abundance and occurrence among PCRs. This can be partly attributed to the

complex composition of soil DNA extracts, with molecules derived from numerous organisms with varying

sizes, abundances, and gene copy numbers. These complex templates are amplified with variable efficiency by

PCR, due to factors including primer/template mismatches, variation of template AT/GC content, polymerase

characteristics and PCR cycling conditions (Acinas et al., 2005; Elbrecht & Leese, 2015; Polz & Cavanaugh,

1998). Low abundance OTUs may also result from PCR or sequencing errors (Leray & Knowlton, 2017), or

over-splitting of OTUs (Clare et al., 2016).

How many PCRs are necessary for estimating soil biodiversity? Previously, Ficetola et al. (2015) suggested

that at least eight PCR replicates are required to correctly infer the presence of specific organisms with a low

detection probability, while Piggott (2016) concluded that two to six PCR replicates were required to

confidently detect an endangered freshwater fish. In this case, the species accumulation curves based on

combining multiple PCRs showed some evidence of flattening, suggesting that between ten and twenty PCRs

may have resulted in maximal OTU detection. Of course, PCR and sequencing artefacts may have affected

these estimates (Leray & Knowlton, 2017). A common strategy to limit these effects is the exclusion of low

abundance sequences or OTUs (Elbrecht & Leese, 2015; Vesty et al., 2017), yet our estimates were derived

from OTU clustering with singleton sequences excluded, and further removal of OTUs with abundances < 5

resulted in very similar conclusions. A more conservative filtering approach enabled by individually indexed

PCRs is to exclude OTUs based on their occurrences in multiple PCRs (Alberdi et al., 2017). In this case,

species accummulation curves based on OTUs that occurred in at least two of each set of PCRs plateaued in

the range of seven to ten, with richness estimates reduced by some 30 %.

Our results show that the optimal number of PCRs will depend on the amenability of the intended targets to

PCR amplification along with assumptions about OTU filtering, with detection of the most readily amplified

organisms assured from a single PCR, but detection of less consistently-amplified targets requiring five to ten

(or more) PCRs. The observed richness estimates were strongly affected by the number of PCRs that were

included. This clearly shows that balancing the number of PCRs that are carried out for different samples is an

This article is protected by copyright. All rights reserved.

 important consideration if species richness comparisons are intended, in addition to the use of consistent
Accepted Article
DNA extraction methods. These results also imply that previous analyses based on triplicate PCRs (and

smaller DNA extraction sizes, e.g. Drummond et al., 2015) are likely to have underestimated the extent of soil

biodiversity.

Our conclusions are based on soil material from a particular forest site, and the extent to which they apply to

soils in different ecosystems and locations is unclear. It may be appropriate to optimise DNA extraction and

PCR methods on a case-by-case basis in the form of pilot studies. Together with other recently published

examinations of sampling, laboratory, and analysis methods (Alberdi et al., 2017; Dickie et al., 2018; Murray et

al., 2015), our results provide a guide for future DNA metabarcoding studies of soil biodiversity.

Acknowledgements

Funding for this research was provided by the New Zealand Tertiary Education Commission in the form of

Centre of Research Excellence funding for the Allan Wilson Centre for Molecular Ecology and Evolution, and

by the University of Auckland in the form of a Doctoral Scholarship received by AD.

Sample collection in the Waitakere Ranges Regional Park was carried out under a research permit provided

by Auckland Council (WS593). We thank Georgianne Griffiths for assistance with sample collection, New

Zealand Genomics Limited (NZGL) for carrying out the Illumina sequencing, and New Zealand eScience

Infrastructure (NeSI) for providing high-performance computing resources and support for this research.

Author Contributions

AD, AJD, TRB and RDN designed the research; AD carried out the sample collection, DNA extractions, PCR

amplifications and sequencing library preparation; AD and DX carried out the bioinformatic processing, OTU

identifications and biodiversity analyses; AD led the writing of the manuscript; All authors contributed to the

drafts and gave final approval for publication.

This article is protected by copyright. All rights reserved.

 Data accessibility
Accepted Article
Illumina sequence data is available from the NCBI Sequence Read Archive under accession SRP148718.

https://www.ncbi.nlm.nih.gov/sra/SRP148718

Analysis code, metadata, and processed data used in this research are available from Github. DOI:

10.5281/zenodo.1403909

References

Acinas, S. G., Sarma-Rupavtarm, R., Klepac-Ceraj, V., & Polz, M. F. (2005). PCR-induced sequence artifacts and bias: Insights

from comparison of two 16S rRNA clone libraries constructed from the same sample. Applied and Environmental

Microbiology, 71(12), 8966–8969. http://doi.org/10.1128/AEM.71.12.8966

Acosta-Martínez, V., Dowd, S., Sun, Y., & Allen, V. (2008). Tag-encoded pyrosequencing analysis of bacterial diversity in a

single soil type as affected by management and land use. Soil Biology and Biochemistry, 40(11), 2762–2770.

http://doi.org/10.1016/j.soilbio.2008.07.022

Alberdi, A., Aizpurua, O., Gilbert, M. T. P., & Bohmann, K. (2017). Scrutinizing key steps for reliable metabarcoding of

environmental samples. Methods in Ecology and Evolution, 2017(June), 1–14.

http://doi.org/10.1111/2041-210X.12849

Amend, A. S., Seifert, K. A., & Bruns, T. D. (2010). Quantifying microbial communities with 454 pyrosequencing: Does read

abundance count? Molecular Ecology, 19(24), 5555–5565. http://doi.org/10.1111/j.1365-294X.2010.04898.x

Bardgett, R. D., & Van Der Putten, W. H. (2014). Belowground biodiversity and ecosystem functioning. Nature 515,

505(7528), 505–511. http://doi.org/10.1038/nature13855

Benjamini, Y., & Hochberg, Y. (1995). Controlling the False Discovery Rate: a Practical and Powerful Approach to Multiple

Testing. Journal of the Royal Statistical Society: Series B (Methodological), 57(1), 289–300.

Buée, M., Reich, M., Murat, C., Morin, E., Nilsson, R. H., Uroz, S., & Martin, F. (2009). 454 Pyrosequencing analyses of forest

soils reveal an unexpectedly high fungal diversity. New Phytologist, 184(2), 449–456.

http://doi.org/10.1111/j.1469-8137.2009.03003.x

Caporaso, J. G., Lauber, C. L., Walters, W. A., Berg-Lyons, D., Lozupone, C. A., Turnbaugh, P. J., … Knight, R. (2011). Global

patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proceedings of the National Academy

of Sciences of the United States of America. http://doi.org/10.1073/pnas.1000080107

Charney, N., & Record, S. (2012). vegetarian: Jost diversity measures for community data. Retrieved from

This article is protected by copyright. All rights reserved.

 https://cran.r-project.org/package=vegetarian

Clare, E. L., Chain, F. J. J., Littlefair, J. E., Cristescu, M. E., & Deiner, K. (2016). The effects of parameter choice on defining
Accepted Article
molecular operational taxonomic units and resulting ecological analyses of metabarcoding data. Genome, 59(11),

981–990. http://doi.org/10.1139/gen-2015-0184

Decaëns, T. (2010). Macroecological patterns in soil communities. Global Ecology and Biogeography, 19(3), 287–302.

http://doi.org/10.1111/j.1466-8238.2009.00517.x

Decaëns, T., Jiménez, J. J., Gioia, C., Measey, G. J., & Lavelle, P. (2006). The values of soil animals for conservation biology.

European Journal of Soil Biology, 42(SUPPL. 1), S23–S38. http://doi.org/10.1016/j.ejsobi.2006.07.001

Delmont, T. O., Robe, P., Cecillon, S., Clark, I. M., Constancias, F., Simonet, P., … Vogel, T. M. (2011). Accessing the soil

metagenome for studies of microbial diversity. Applied and Environmental Microbiology, 77(4), 1315–24.

http://doi.org/10.1128/AEM.01526-10

Dickie, I. A., Boyer, S., Buckley, H. L., Duncan, R. P., Gardner, P. P., Hogg, I. D., … Weaver, L. (2018). Towards robust and

repeatable sampling methods in eDNA-based studies. Molecular Ecology Resources, (May).

http://doi.org/10.1111/1755-0998.12907

Drummond, A. J., Newcomb, R. D., Buckley, T. R., Xie, D., Dopheide, A., Potter, B. C. M., … Nelson, N. (2015). Evaluating a

multigene environmental DNA approach for biodiversity assessment. GigaScience, 4(1), 46.

http://doi.org/10.1186/s13742-015-0086-1

Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nature Methods, 10(10),

996–8. http://doi.org/10.1038/nmeth.2604

Egge, E., Bittner, L., Andersen, T., Audic, S., de Vargas, C., & Edvardsen, B. (2013). 454 pyrosequencing to describe microbial

eukaryotic community composition, diversity and relative abundance: a test for marine haptophytes. PLoS ONE,

8(9), e74371. http://doi.org/10.1371/journal.pone.0074371

Elbrecht, V., & Leese, F. (2015). Can DNA-based ecosystem assessments quantify species abundance? Testing primer bias

and biomass-sequence relationships with an innovative metabarcoding protocol. PLoS ONE, 10(7), 1–16.

http://doi.org/10.1371/journal.pone.0130324

Ficetola, G. F., Pansu, J., Bonin, A., Coissac, E., Giguet-Covex, C., De Barba, M., … Taberlet, P. (2015). Replication levels, false

presences and the estimation of the presence/absence from eDNA metabarcoding data. Molecular Ecology

Resources, 15(3), 543–556. http://doi.org/10.1111/1755-0998.12338

Fonseca, V., Nichols, B., & Lallias, D. (2012). Sample richness and genetic diversity as drivers of chimera formation in nSSU

metagenetic analyses. Nucleic Acids Research, 40(9), e66. http://doi.org/10.1093/nar/gks002

Freckman, D. W., & Virginia, R. A. (1993). Extraction of nematodes from Dry Valley Antarctic soils. Polar Biology, 13(7),

This article is protected by copyright. All rights reserved.

 483–487. http://doi.org/10.1007/BF00233139

Giller, P. S. (1996). The diversity of soil communities, the “poor man’s tropical rainforest.” Biodiversity and Conservation,
Accepted Article
5(2), 135–168. http://doi.org/10.1007/BF00055827

Henderson, G., Cox, F., Kittelmann, S., Miri, V. H., Zethof, M., Noel, S. J., … Janssen, P. H. (2013). Effect of DNA extraction

methods and sampling techniques on the apparent structure of cow and sheep rumen microbial communities. PloS

One, 8(9), e74787. http://doi.org/10.1371/journal.pone.0074787

Huson, D. H., Mitra, S., Ruscheweyh, H. J., Weber, N., & Schuster, S. C. (2011). Integrative analysis of environmental

sequences using MEGAN4. Genome Research, 21(9), 1552–1560. http://doi.org/10.1101/gr.120618.111

Ji, Y., Ashton, L., Pedley, S. M., Edwards, D. P., Tang, Y., Nakamura, A., … Yu, D. W. (2013). Reliable, verifiable and efficient

monitoring of biodiversity via metabarcoding. Ecology Letters, 16(10), 1245–57. http://doi.org/10.1111/ele.12162

Jost, L. (2006). Entropy and diversity. OIKOS, 113(2), 363–375.

Kang, S., & Mills, A. L. (2006). The effect of sample size in studies of soil microbial community structure. Journal of

Microbiological Methods, 66(2), 242–50. http://doi.org/10.1016/j.mimet.2005.11.013

Kurtzman, C. P., & Robnett, C. J. (1998). Identification and phylogeny of ascomycetous yeasts from analysis of nuclear large

subunit (26S) ribosomal DNA partial sequences. Antonie van Leeuwenhoek, 73(4), 331–71. Retrieved from

http://www.ncbi.nlm.nih.gov/pubmed/9850420

Leray, M., & Knowlton, N. (2017). Random sampling causes the low reproducibility of rare eukaryotic {OTUs} in {Illumina}

{COI} metabarcoding. PeerJ, 5, e3006. http://doi.org/10.7717/peerj.3006

Leray, M., Yang, J. Y., Meyer, C. P., Mills, S. C., Agudelo, N., Ranwez, V., … Machida, R. J. (2013). A new versatile primer set

targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for

characterizing coral reef fish gut contents. Frontiers in Zoology, 10(1), 34. http://doi.org/10.1186/1742-9994-10-34

Lombard, N., Prestat, E., van Elsas, J. D., & Simonet, P. (2011). Soil-specific limitations for access and analysis of soil

microbial communities by metagenomics. FEMS Microbiology Ecology, 78(1), 31–49.

http://doi.org/10.1111/j.1574-6941.2011.01140.x

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-Seq data with

DESeq2. Genome Biology, 15, 550. http://doi.org/10.1186/s13059-014-0550-8

Machida, R. J., & Knowlton, N. (2012). PCR Primers for metazoan nuclear 18S and 28S ribosomal DNA sequences. PLoS

ONE, 7(9), e46180. http://doi.org/10.1371/journal.pone.0046180

McMurdie, P. J., & Holmes, S. (2013). phyloseq: an R package for reproducible interactive analysis and graphics of

microbiome census data. PLoS ONE, 8(4), e61217. http://doi.org/10.1371/journal.pone.0061217

McMurdie, P. J., & Holmes, S. (2014). Waste not, want not: why rarefying microbiome data is inadmissible. PLoS

This article is protected by copyright. All rights reserved.

 Computational Biology, 10(4), e1003531. http://doi.org/10.1371/journal.pcbi.1003531

Miller, D. N., Bryant, J. E., Madsen, E. L., & Ghiorse, W. C. (1999). Evaluation and optimization of DNA extraction and
Accepted Article
purification procedures for soil and sediment samples. Applied and Environmental Microbiology, 65(11), 4715–24.

Retrieved from

http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=91634&tool=pmcentrez&rendertype=abstract

Murray, D. C., Coghlan, M. L., & Bunce, M. (2015). From benchtop to desktop: important considerations when designing

amplicon sequencing workflows. PLoS ONE, 10(4), e0124671. http://doi.org/10.5061/dryad.2qf0t

Oksanen, J., Blanchet, F. Guillaume Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., … Wagner, H. (2017).

vegan: Community Ecology Package. Retrieved from http://cran.r-project.org/package=vegan

Piggott, M. P. (2016). Evaluating the effects of laboratory protocols on eDNA detection probability for an endangered

freshwater fish. Ecology and Evolution, 6(9), 2739–2750. http://doi.org/10.1002/ece3.2083

Polz, M. F., & Cavanaugh, C. M. (1998). Bias in template-to-product ratios in multitemplate PCR. Applied and Environmental

Microbiology, 64(10), 3724–30. Retrieved from

http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=106531&tool=pmcentrez&rendertype=abstract

Porazinska, D. L., Giblin-Davis, R. M., Faller, L., Farmerie, W., Kanzaki, N., Morris, K., … Thomas, W. K. (2009). Evaluating

high-throughput sequencing as a method for metagenomic analysis of nematode diversity. Molecular Ecology

Resources, 9(6), 1439–1450. http://doi.org/10.1111/j.1755-0998.2009.02611.x

Qiu, X., Wu, L., Huang, H., Donel, P. E. M. C., Palumbo, A. V, Tiedje, J. M., & Zhou, J. (2001). Evaluation of PCR-generated

chimeras, mutations, and heteroduplexes with 16S rRNA gene-based cloning. Applied and Environmental

Microbiology, 67(2), 880–887. http://doi.org/10.1128/AEM.67.2.880

R Core Team. (2016). R: A language and environment for statistical computing. Vienna: R Foundation for Statistical

Computing. Retrieved from https://www.r-project.org/

Ranjard, L., Lejon, D. P. H., Mougel, C., Schehrer, L., Merdinoglu, D., & Chaussod, R. (2003). Sampling strategy in molecular

microbial ecology: Influence of soil sample size on DNA fingerprinting analysis of fungal and bacterial communities.

Environmental Microbiology, 5(11), 1111–1120. http://doi.org/10.1046/j.1462-2920.2003.00521.x

RStudio team. (2015). RStudio: Integrated development environment for R. Boston, MA: RStudio, Inc. Retrieved from

http://www.rstudio.org

Stevens, J. L., Jackson, R. L., & Olson, J. B. (2013). Slowing PCR ramp speed reduces chimera formation from environmental

samples. Journal of Microbiological Methods, 93(3), 203–205. http://doi.org/10.1016/j.mimet.2013.03.013

Taberlet, P., Prud’Homme, S. M., Campione, E., Roy, J., Miquel, C., Shehzad, W., … Coissac, E. (2012). Soil sampling and

isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies. Molecular

This article is protected by copyright. All rights reserved.

 Ecology, 21(8), 1816–20. http://doi.org/10.1111/j.1365-294X.2011.05317.x

Tedersoo, L., Nilsson, R. H., Abarenkov, K., Jairus, T., Sadam, A., Saar, I., … Kõljalg, U. (2010). 454 Pyrosequencing and
Accepted Article
Sanger sequencing of tropical mycorrhizal fungi provide similar results but reveal substantial methodological

biases. New Phytologist, 188(1), 291–301. http://doi.org/10.1111/j.1469-8137.2010.03373.x

Torsvik, V., & Øvreås, L. (2002). Microbial diversity and function in soil: from genes to ecosystems. Current Opinion in

Microbiology, 5(3), 240–5. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/12057676

Velasco-Castrillón, A., Schultz, M. B., Colombo, F., Gibson, J. a E., Davies, K. A., Austin, A. D., & Stevens, M. I. (2014).

Distribution and diversity of soil microfauna from East Antarctica: assessing the link between biotic and abiotic

factors. PLoS ONE, 9(1), e87529. http://doi.org/10.1371/journal.pone.0087529

Vesty, A., Biswas, K., Taylor, M. W., Gear, K., & Douglas, R. G. (2017). Evaluating the impact of DNA extraction method on

the representation of human oral bacterial and fungal communities. PLoS ONE, 12(1), 1–13.

http://doi.org/10.1371/journal.pone.0169877

Vishnivetskaya, T. A., Layton, A. C., Lau, M. C. Y., Chauhan, A., Cheng, K. R., Meyers, A. J., … Onstott, T. C. (2014). Commercial

DNA extraction kits impact observed microbial community composition in permafrost samples. FEMS Microbiology

Ecology, 87(1), 217–230. http://doi.org/10.1111/1574-6941.12219

Wagner, A. O., Praeg, N., Reitschuler, C., & Illmer, P. (2015). Effect of DNA extraction procedure, repeated extraction and

ethidium monoazide (EMA)/propidium monoazide (PMA) treatment on overall DNA yield and impact on microbial

fingerprints for bacteria, fungi and archaea in a reference soil. Applied Soil Ecology, 93, 56–64.

http://doi.org/10.1016/j.apsoil.2015.04.005

Wickham, H. (2009). ggplot2: elegant graphics for data analysis. New York, NY: Springer New York.

Zhang, J., Kobert, K., Flouri, T., & Stamatakis, A. (2014). PEAR: a fast and accurate Illumina Paired-End reAd mergeR.

Bioinformatics, 30(5), 614–20. http://doi.org/10.1093/bioinformatics/btt593

Zhou, J., Wu, L., Deng, Y., Zhi, X., Jiang, Y., Tu, Q., … Yang, Y. (2011). Reproducibility and quantitation of amplicon

sequencing-based detection. The ISME Journal, 5(8), 1303–1313. http://doi.org/10.1038/ismej.2011.11

Zinger, L., Chave, J., Coissac, E., Iribar, A., Louisanna, E., Manzi, S., … Taberlet, P. (2016). Extracellular DNA extraction is a

fast, cheap and reliable alternative for multi-taxa surveys based on soil DNA. Soil Biology and Biochemistry, 96(May),

16–19. http://doi.org/10.1016/j.soilbio.2016.01.008

This article is protected by copyright. All rights reserved.

 Tables
Accepted Article
Table 1: PERMANOVA analysis statistics for communities detected by metabarcoding analysis of four

amplicons and four DNA extraction methods applied to 15 g soil samples from two subplots separated by

approximately 21 m. Values are PERMANOVA F statistics with p-values in parentheses.

Distance metric Factor 16S 18S 26S COI

subplots 161.115 (≤0.001) 64.771 (≤0.001) 55.958 (≤0.001) 74.386 (≤0.001)
Effective beta methods 150.964 (≤0.001) 41.795 (≤0.001) 29.908 (≤0.001) 18.315 (≤0.001)
subplots:methods -6.966 (1) 9.614 (≤0.001) 4.197 (≤0.001) 3.368 (0.003)
subplots 9.1297 (≤0.001) 9.58 (≤0.001) 13.3305 (≤0.001) 15.5486 (≤0.001)
Jaccard methods 4.8829 (≤0.001) 6.1941 (≤0.001) 5.4355 (≤0.001) 4.6277 (≤0.001)
subplots:methods 2.3076 (≤0.001) 2.8773 (≤0.001) 3.0168 (≤0.001) 2.8965 (≤0.001)

Figure legends

Figure 1: Alpha, beta and gamma 0 (richness) and 1 (effective diversity, calculated as the

exponential of Shannon entropy; Jost, 2006) estimates for 18S phyla detected in the results of four

DNA extraction methods and sets of ten PCRs applied to 15 g soil subsamples from two subplots.

Each data point represents the estimated mean diversity across ten PCRs, based on data randomly

subsampled to equal read depths for each sample, method, and PCR. The subsampling and diversity

estimates were repeated ten times. Forward (/) and back (\) slashes represent separate estimates

for samples from subplots D and P, respectively. Only phyla with mean estimated species richness ≥

8 are shown. Chytridiomyc. and Basidiomyc. denote Chytridiomycota and Basidiomycota,

respectively.

This article is protected by copyright. All rights reserved.

 Figure 2: The proportions and taxonomic composition of OTUs per amplicon with total abundances
Accepted Article
≥ 2, 3, 5, or 10 that were detected in the results of different DNA extraction methods and sets of ten

PCRs applied to 15 g soil subsamples from two subplots. Labels on the y axis indicate the

combinations of DNA extraction methods in which OTUs were detected. A: 1.5 g direct extractions;

B: 7.5 g direct extractions; C: 15 g indirect extractions; D: 15 g phosphate buffer extractions; All: all

four extraction methods.

Figure 3: Multivariate community structure (effective beta distance) of (a) 16S prokaryote, (b) 18S

protist, (c) 18S fungal, (d) 18S metazoan, (e) 26S fungal, and (f) COI metazoan OTU assemblages

detected by four DNA extraction methods and sets of ten PCRs applied to 15 g soil subsamples from

two subplots separated by approximately 21 m.

Figure 4: Pairwise multivariate sample similarity (effective beta distance) of OTU assemblages and

major taxonomic groups detected by four DNA extraction methods and sets of ten PCRs applied to

15 g soil subsamples from within the same subplot (red) or from different subplots separated by

approximately 21 m (blue).

Figure 5: The abundance and occurrence between sets of ten PCRs for 16S, 18S, 26S and COI

amplicons from four DNA extraction methods applied to 15 g soil subsamples from two subplots.

Figure 6: Alpha diversity estimates resulting from pooling combinations of one to ten PCRs (left) or

from subsampling an equal number of sequence reads from pooled combinations of one to ten PCRs

(right) for 16S, 18S, 26S and COI genes from four DNA extraction methods applied to 15 g soil

subsamples from subplot D. Gold and blue points represent estimated species richness and effective

alpha diversity (calculated as the exponential of Shannon entropy; Jost, 2006), respectively.

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.