Metabarcoding reveals a high diversity of woody host-associated Phytophthora spp. in soils at public gardens and amenity woodlands in Britain

Soil sampling and DNA extraction

Fourteen sites in northern Britain were sampled between October 2014 and September 2015. These included eleven public gardens/arboreta, two public cemeteries and one privately owned upland juniper woodland with no public access. With the exception of one of the public gardens, all sites had previously confirmed Phytophthora outbreaks. The geographical coordinates (not shown here) and approximate altitude (Table S1) were determined for each site from a roughly central location. For ten of the sites the underlying soil type was determined using the Soil Information for Scottish Soils web interface http://sifss.hutton.ac.uk/ (Table S1). For the remaining four sites (sites 2, 4, 5 and 14) the underlying soil type could not be determined due to their location in urban areas (Table S1). At each site, soil samples were collected from around ten trees/shrubs either showing visible symptoms typical of Phytophthora infection such as stem bleeding and foliage browning/bronzing, or which were outwardly healthy but known to be potential Phytophthora hosts. Soil was also sampled from around stumps of trees/shrubs recently felled due to Phytophthora disease. Thus the locations of sampling at each site varied according to the above criteria, with sampling points distributed over as wide an area on each site as possible. For each sampling point the geographical coordinates were recorded as well as the genus/species of tree/shrub present and its condition, either; i) dead, stump ii) dead, standing, iii) live, but showing symptoms typical of Phytophthora infection and iv) live, outwardly healthy. Table S1 shows the number of healthy vs diseased trees/shrubs, and their broad taxonomic grouping, that formed each soil sampling point across the 14 sites.

Each individual soil sample (approximately 300–400 g) was comprised of eight pooled soil cores of 2 cm width × 30 cm depth collected using a soil auger whereby two cores were collected at four points around each individual tree/shrub/stump within one m horizontal distance from the root collar. The eight subsamples were pooled and homogenised by hand in a single grip-seal™ polythene bag and frozen within 12 h of collection. Each soil sample was oven dried at ∼60 °C in aluminium trays for 1–3 days (depending on soil wetness), stirred thoroughly once dry and DNA extracted from three 250 mg subsamples using the PowerSoil® DNA isolation kit (MoBio Laboratories Inc., Carlsbad, CA, USA). A robotic workstation for DNA extraction based on magnetic-particle purification (Kingfisher™ mL Magnetic Particle Processor, Thermo Scientific, Loughborough, UK) was used for the DNA extraction process. Post-DNA extraction clean-up was carried out using either the Jet-Quick™ DNA purification kit (Genomed GmbH, Löhne, Germany) or DNA Clean & Concentrator™ (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions.

Development of positive controls for all Illumina sequencing runs

For the first sequencing run a control genomic DNA mix was set up containing the following 15 Phytophthora species; P. austrocedri, P. cactorum, P. cambivora, P. chlamydospora, P. cinnamomi, P. gonapodyides, P. ilicis, P. kernoviae, P. lateralis, (EU2 lineage), P. obscura, P. plurivora, P. pseudosyringae, P. ramorum, P. syringae and P. boehmeriae (all but the last known to be present in the UK). For the second sequencing run a control genomic DNA mix was set up containing the following ten Phytophthora species; P. cactorum, P. foliorum, P. obscura, P. plurivora, P. rubi, P. siskiyouensis (all known to be present in the UK), and P. boehmeriae, P. castaneae, P. capsici and P. fallax (not known to be present in the UK). Species were included in the control mixes based on being present in Britain and expected to be detected in soil samples, as well as species not present in Britain and therefore not expected to be detected. Species were also chosen to span as wide a range of clades as possible from within these two criteria, with seven clades represented in both mixes. All cultures originated from single hyphal-tip colonies. Genomic DNA was extracted from mycelium harvested from Phytophthora cultures on V8 agar using the Nucleospin Plant Prep II kit (Machery-Nagel, Düren, Germany) according to the manufacturer’s instructions. A SYBR green quantitative real-time PCR assay was carried out on an ABI Prism 7300 Real Time PCR System (Applied Biosystems, Foster City, CA, USA) with the primers 18Ph2F and 5.8S-1R (Scibetta et al., 2012) using DNA of each Phytophthora species diluted 5, 10, 25, 50, 100 and 500×. Cycle threshold (Ct) values at each concentration were determined and DNAs from each species corresponding to Ct values between 22 and 30 were pooled for the DNA control mix; thus species were not expected to be present in equal abundance.

Amplicon PCR

A ∼ 250 bp region of the ribosomal RNA (rRNA) internal transcribed spacer (ITS1) was amplified from each DNA sample using nested PCR with primer pairs 18Ph2F and 5.8S-1R in the first round and ITS6 and 5.8S-1R in the second round according to the protocol of Scibetta et al. (2012), except that proof-reading enzyme KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA) was used for the PCR reaction to minimise errors during PCR. Second round primers were amended with overhang adapters to ensure compatibility with the Illumina index and sequencing adapters. These were: forward overhang; 5′ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG‐ [ITS6] and reverse overhang; 5′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG‐ [5.8S-1R] (Illumina, 2013). For each DNA sample, PCR was carried out in triplicate and all Phytophthora-positive PCR replicates were pooled for downstream processing.

Illumina sequencing library preparation and sequencing

Samples were prepared for sequencing following the protocols described for 16S Metagenomic Sequencing Library Preparation (Illumina, 2013). In brief, this involved clean-up of amplicon PCR using Agencourt® AMpure® XP beads (Agencourt Bioscience, Beverly, MA, USA) followed by index PCR in KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA) to attach dual indices and Illumina sequencing adapters to each sample using the Nextera XT Index Kit (Epicentre, Madison, WI, USA). This step ensured that each sample could be uniquely identified during the sequencing run. A second PCR clean-up (as above) was then carried out and DNA of each sample visualised on an Agilent 2200 TapeStation (Santa Clara, CA, USA) and quantified using a Qubit® fluorimeter (Life Technologies, Paisley, UK). For each sequencing run, 96 samples were pooled for paired-end (2 × 250 bp) sequencing on a single flow cell of an Illumina MiSeq sequencer using the MiSeq v. 2 500 bp Nano kit (Illumina, San Diego, CA, USA) at the James Hutton Institute, Dundee, UK. Following quality control and de-multiplexing, FASTQ files containing reads for each sample were exported for bioinformatics analysis.

Bioinformatics analysis

Sequence data from all 14 sites were processed in an attempt to assign them to clusters of highly similar sequences, referred to hereafter as operational taxonomic units (OTUs), and assigned species identity using the bioinformatics software ‘metapy’ (https://github.com/peterthorpe5/public_scripts/tree/master/metapy) (github commit: 6fd1864). The version of metapy used in this study performed read quality trimming using Trimmomatic (version 0.36) (Bolger, Lohse & Usadel, 2014). Reads were then assembled using PEAR (version 0.9.10) (Zhang et al., 2013). Metapy subsequently ran two sequence analysis tools which performed independently of each other. These were Swarm (version 1.2.19) (Mahé et al., 2014) and Bowtie (version 2.2.5) (Langmead, 2010). Swarm calculates sequence differences between pairs of reads to delineate OTUs using k-mer comparisons and a global pairwise alignment algorithm. Clustering results are refined based on sequence abundance and OTU structure (Mahé et al., 2014). For the analyses reported here the Swarm parameters set to delineate OTUs were ‘t -1 d -1’ and corresponded to approximately one to two mismatches across the sequence. Bowtie is a read-mapping tool and metapy filtered the resulting output for perfect matches. Species identification was based on a version of the custom curated Phytophthora reference database developed by Català, Pérez-Sierra & Abad-Campos (2015) which was updated to include more recently described Phytophthora species plus additional ITS1 sequences for species where intra-species ITS1 variability has been documented. This database version was dated June 2017 (https://github.com/peterthorpe5/public_scripts/blob/master/metapy/data/Riddell_et_al_2019.fasta). Any sequences not assigned species identity from the reference database were compared to the GenBank nt database using BLASTN+ (Altschul et al., 1990) to find the nearest matching sequences. Sequences that did not match to the reference database but which were close to Phytophthora sequences in the GenBank nt database were subjected to phylogenetic analyses to show their position within their clade. Table S3 shows the GenBank accession numbers and isolate identifiers for sequences used in the phylogenetic analyses. Maximum likelihood (ML) trees were estimated from nucleotide sequences using PhyML 3.0 (Guindon et al., 2010), with the general time reversible (GTR) model of sequence change and gamma-distributed rate variation across sites (four categories); 1,000 bootstrap replicates were examined. Midpoint-rooted phylogenies were displayed using FigTree (http://tree.bio.ed.ac.uk/software/figtree/).

Statistical analyses

To improve the balance of the data across sites and reduce the number of factors, host health status was defined as ‘healthy’ or ‘symptoms/stump/dead’ (initial analysis indicated no significant differences among symptomatic, stump or dead samples). Table S1 shows the counts of samples, broken down by the relevant site-level and within-site factors considered within the statistical analysis. The ordination of Phytophthora assemblages was conducted using both abundance data (i.e. number of reads) and absence/presence of Phytophthora species. Phytophthora species were excluded from the ordination analysis when they appeared in only a single sample, and samples were excluded where no Phytophthora species were detected. For abundance data, non-metric multidimensional scaling (NMDS) was conducted on the Bray–Curtis dissimilarity matrix, with centring and Wisconsin square-root transformations. For absence/presence data, a Jaccard distance matrix was calculated using the binary data set. In both cases, ordination was conducted using the metaMDS function in the vegan package version 2.5-2 in R (Oksanen et al., 2018), with three dimensions required to produce acceptable stress levels (stress = 0.153 and 0.126 respectively). For statistical analysis, permANOVA was conducted on matrices for both data sets using the adonis functions in the vegan package in R: for drivers of site-level differences, latitude, longitude, altitude and underlying soil type were used as potential factors and variables, with sites used as strata within the models. A separate analysis was subsequently conducted for the effects of host status and taxonomic group within site, with marginal effects calculated in each case.

Generalised linear mixed-effects models (GLiMMs) were also applied to the data set for all samples (n = 140), but restricted to Phytophthora species present in more than 10 samples, as inclusion of rarer species prevented valid model fit. Data were analysed as absence/presence (i.e. binary) data of each individual Phytophthora species, with binomial error structure and probit link function, with host status/Phytophthora species interaction, host taxonomic group and site-level factors/variables as fixed effects, and sample nested within site as random effects. Analysis of variance (Wald chi-square tests) was applied to the model to assess the significance of the main effects and interactions; non-significant effects were subsequently removed to determine the best-fit model. The best-fit model residuals were assessed using the DHARMa package version 0.2.0 in R (Hartig, 2018). Post hoc analysis (Tukey’s HSD) was used to determine significant differences within significant effects, with outputted marginal means, i.e. averaged over significant factors, representing the proportion of samples with Phytophthora species present.

Baiting

Soils were additionally baited for Phytophthora from one of the botanic garden sites (site 12; Table S1). Soil samples were collected in June 2016 from around five of the trees/shrubs originally sampled in October 2014. These sampling points were chosen based on having yielded a diverse mix of Phytophthora sequences in the metabarcoding analysis. For each baited sample approximately 250 g of soil was obtained by pooling four soil cores collected within one horizontal metre of the tree/shrub stem base. Samples were stored at 4 °C and baiting carried out on the soil at two different states; i) after 1 week storage at 4 °C, and ii) air dried at room temperature (18–22 °C) after the first baiting test followed by 5 weeks storage at 4 °C. Soil was baited inside wounded apples of the variety ‘Granny Smith’, two apples per soil sample, with one apple incubated at room temperature and the other incubated at 16 °C. Soil baiting with plant material was also set up whereby each soil sample was placed in a plastic sandwich box which was filled to twice the depth of the soil with distilled water and left to settle for 24 h. Any floating soil material was then removed and intact leaves of the following species placed to float on the water in each sample box as follows; 2 × Rhododendron ponticum leaves, 2× Quercus robur leaves, 6× Hebe leaves, 3× Hebe cuttings, 1× Juniperus communis cutting and one leaf each of Solanum lycopersicum, Capsicum annuum, Petunia sp., Dianthus sp. and F. sylvatica. Controls were also set up in which each baiting test was replicated in the absence of soil.

Bait plants were examined every 2 days for lesion development. Tissue from lesion margins was plated onto Synthetic Mucor Agar of Elliot, Hendrie & Knights (1966) with the addition of antibiotics as described in Brasier et al. (2005) and incubated at 16 °C. Developing colonies were subcultured onto V8 juice agar (V8A; 2 g CaCO₃, 200 ml V8 juice and 15 g agar in 800 ml distilled water), incubated at 16 °C and, once grown, three 1 cm³ plugs were cut from the edge of each colony and transferred to Petri plates containing fresh pond water pre-filtered through Whatman No. 1 filter paper (GE Healthcare, Little Chalfont, UK). The plates containing pond water were incubated at 16 °C for 3–4 days and checked for the presence of Phytophthora-like sporangia.

DNA was extracted from mycelium harvested from Phytophthora-like cultures on V8 agar using the Nucleospin Plant Prep II kit (Machery-Nagel, Düren, Germany) according to the manufacturer’s instructions. PCR was performed using one µl of either the Phytophthora-specific forward primers Ph2 (Ippolito, Schena & Nigro, 2002) or ITS6 (White et al., 1990), and the universal reverse primer ITS4 (White et al., 1990) at a concentration of 10 mM. Total reaction volume was 25 µl comprising 1.5 µl MgCl₂ (at 0.45 mM), five µl of 5× buffer, 0.5 µl dNTPs (at 0.2 mM), 0.125 µl U Taq DNA polymerase (Bioline, London, UK), 17.5 µl molecular grade water and one µl of template DNA. Amplification was performed in a Veriti 96 well Thermalcycler (Applied Biosystems, Foster City, CA, USA) with initial denaturation at 95 °C for 5 min followed by 35 cycles of 94 °C for 30 s, 55 °C for 30 s, 72 °C for 1 min and a final extension of 72 °C for 7 min. PCR products were purified and sequenced in both directions with the BigDye version 3.1 Ready Reaction Kit on an ABI Prism 3730 capillary sequencer (Applied Biosystems, Foster City, CA, USA). Raw sequences were aligned and edited using Sequencher version 4.8 for Windows, and searched against published ITS sequences in the GenBank nt database using BLASTN+ (Altschul et al., 1990).

Sequencing output and performance of clustering tools in the control reactions

The two sequencing runs, which included samples from the 14 sites and two control reactions, together generated >400,000 good quality sequences that could be considered for analysis. A total of 7,123 assembled sequences were generated from the control reaction comprising 15 Phytophthora species included in the first sequencing run. In the second sequencing run a total of 4,459 assembled sequences were generated from the control reaction comprising ten Phytophthora species. Using the metapy software, the two sequence analysis tools Bowtie and Swarm reported different numbers of true positive and false positive Phytophthora species for each control reaction, and different numbers of species not detected in each reaction (Table 1).

Using the Bowtie tool, species included in the control mixes but not identified were P. boehmeriae, P. gonapodyides, P. ilicis and P. syringae in the first control mix and P. boehmeriae, P. cactorum and P. siskiyouensis in the second control mix. False positive species assigned to OTUs using Bowtie were P. idaei in the first control reaction and P. glovera, P. agathacidica, P. sp. novaeguineae, P. citricola, P. fragariae and P. idaei in the second control reaction (Table 2). The number of reads allocated to an OTU using Bowtie ranged from 258 to one in the first control reaction and 192 to one in the second control reaction (Table 2). For two of the OTUs the false positives were due to shared ITS1 sequences preventing discrimination of very closely related species. These were P. glovera, which has an identical ITS1 sequence to P. capsici, and P. agathacidica and P. sp. novaeguineae which have identical ITS1 sequences to P. castanae. Since Bowtie was used with settings that assign species to OTUs only when there are no mismatches across the alignment, the remaining false positive results were potentially ITS1 sequence variants within the genomes of individual strains included in the control mix. Excluding P. boehmeriae, which is discussed later, the false-negative results and overall low number of reads assigned to each OTU using Bowtie are probably due to differences between the reference sequences in the database and sample DNA. These points are enlarged upon in the discussion.

Using Swarm as the clustering tool, two of the species included in the first control reaction, P. pseudosyringae and P. ilicis, could not be distinguished but rather clustered together into a single OTU with 1,595 reads (Table 2). This was due to high sequence similarity as the ITS1 database sequences for these two species differ by only two single nucleotide deletions in P. ilicis; one is a deleted A in a run of As and one a deleted C. Two other closely related species with highly similar ITS1 sequences, P. pluvialis and P. nemorosa, were also falsely assigned to the P. pseudosyringae and P. ilicis OTU (Table 2). With the exception of P. boehmeriae, OTUs were assigned to all remaining Phytophthora species included in the first control reaction with the number of reads generated per OTU varying from 1,373 to two (Table 2). A false positive, P. idaei, which has a highly similar ITS1 sequence to P. cactorum, was assigned to the P. cactorum OTU, and an OTU containing four reads was also assigned to a false positive species, P. mississippiae (Table 2). These results are likely due to ITS1 sequence variation within the isolates of P. cactorum and P. gonapodyides included in the control mix, as discussed later. In the second control reaction, P. boehmeriae reads were not detected but all remaining species included in the control mix were assigned to OTUs, in some cases together with falsely assigned species of high ITS1 sequence similarity, with read numbers per OTU varying from 1,075 to 33 (Table 2).

Since P. boehmeriae sequences were not detected in either control reaction, PCR and Sanger sequencing was used to obtain the ITS1 sequence from the P. boehmeriae isolate included in the control mixes. The resulting sequence was identical to the P. boehmeriae sequence in the reference database, thus ruling out a sequence mismatch as being the cause for the failure to assign sequences to P. boehmeriae. Assembled sequences which had not clustered to the database using Swarm (approximately 8% of the total number of assembled sequences in each control reaction) were compared to sequences in GenBank using BLASTN+ (Altschul et al., 1990), but none matched any sequences of P. boehmeriae. Subsequent Illumina sequencing and analysis of reads from a pure culture of P. boehmeriae resulted in high numbers of reads being matched to the P. boehmeriae sequence in the reference database using the Swarm tool (D.E.L. Cooke, L. Pritchard, P. Thorpe, E. Randall and B. Clark, 2018, unpublished data) indicating the species was detectable using the pipeline. However, there are 3 and 8 bp differences, respectively, between the P. boehmeriae ITS sequence and the second round PCR primer sequences, ITS6 F and 5.8S-1R (Scibetta et al., 2012). Inefficient primer annealing may therefore have limited the amplification of the P. boehmeriae ITS1 region when included as a component of these species mixes.

The outputs obtained using Swarm as the clustering tool (with a threshold set at 1–2 mismatches across the sequence) allowed for the identification of a greater number of species included in each control reaction compared with the outputs generated using the read-mapping tool Bowtie, which was less sensitive but had fewer false positives (Tables 1 and 2). For example, in some cases where the Swarm output assigned multiple species to an OTU, the Bowtie output assigned a single species to the corresponding OTU (e.g. P. pseudosyringae and P. cactorum/P. idaei in the first control mix and P. plurivora and P. rubi in the second control mix) (Table 2). Therefore, in attempting to capture the full range of diversity in the analysis of sequences from the soil samples, the species data reported in this study are those generated using Swarm as the clustering tool, excluding all singleton sequences. The output generated using the Bowtie tool was taken into consideration when multiple closely related species were assigned to an individual OTU using Swarm. If the Bowtie output from the same sample assigned the reads to just one of these species then only this species was included in the data.

Analysis of sequences from soil samples

Across all 140 soil samples analysed from 14 sites, 112 samples and all sites yielded an amplification product in the nested PCR, with sequences assigned to 35 distinct OTUs in total (Table 3). Trees/shrubs representing 34 plant genera yielded Phytophthora-positive soils from around their bases. Of these, approximately 25% were outwardly healthy, 40% showed foliage symptoms typical of Phytophthora infection and 35% were either dead or were stumps of hosts felled prior to sampling due to previously confirmed Phytophthora infection.

The permANOVA analysis showed that site-level factors and variables (latitude, longitude, altitude and underlying soil type) were not significant drivers of the variance in the distance matrices in either the Phytophthora DNA abundance or absence/presence data (Table S3). Host taxonomic group and health status were not significant drivers within-site, and the only significant factor was site itself, which was highly significant for both the abundance data (F_13,98 = 2.93, p < 0.001) and the presence/absence data (F_13,98 = 2.92, p < 0.001); Fig. 1A plots the site data onto the first two axes of the NMDS for the abundance data, highlighting in colour sites 1, 5, 6, 8, 9 and 10 which have markedly different ordinations compared with the rest of the sites. Figure 1B presents the NMDS plot for the Phytophthora species presence/absence data, highlighting in colour site 8 as markedly different compared with other sites. All sites highlighted in the NMDS plots were found to have within their assemblages either unique species (not found elsewhere in this study) or a high abundance of species consistently found in lesser abundance at other sites. Site 8 in particular yielded DNA of species such as P. psychrophila, P. cryptogea and two unknown species not found in abundance elsewhere (Figs. 1A and 1B; Table 3).

Generalised linear mixed-effects models also showed no improvement in model fit with the inclusion of site-level variables (latitude, longitude, altitude and underlying soil type), host health status (healthy vs stump/dead/symptomatic) or host taxonomic group on the proportion of samples with DNA of Phytophthora species present (applied to those species found in >10 samples in the study) (Table S3). There were significant differences among Phytophthora species in terms of overall DNA presence/absence within a sample (Wald χ² = 129, df = 7, N = 140, p < 0.0001, see Table 4 for post hoc test results).

Across all sites and samples, DNA of P. pseudosyringae/P. nemorosa/P. pluvialis/P. ilicis and P. austrocedri was reported significantly more often in the samples than the six remaining species included in the analysis, except for P. gonapodyides which occurred at a similar frequency to P. austrocedri (Table 4). Given that neither P. pluvialis or P. nemorosa are known to occur in the UK (Hansen et al., 2003; Dick et al., 2014; Brar et al., 2017) and P. ilicis is only known to infect species within the genus Ilex (Buddenhagen & Young, 1957), whereas P. pseudosyringae has been reported previously at 23 separate locations in England, Wales and Scotland (Scanu & Webber, 2016), including one of the sites sampled in this study, we consider that the reads assigned by Swarm to P. pseudosyringae/P. pluvialis/P. nemorosa/P. ilicis most likely represent P. pseudosyringae. This assumption is additionally supported by the assignment of the majority of these sequences to P. pseudosyringae by the read-mapping programme Bowtie, for which metapy assigns reads to species only when the alignment has no mismatches.

Other frequently occurring species identified in the study were P. cambivora, P. cactorum, P. ramorum, P. syringae and P. chlamydospora (Table 4). Of the other species identified in the study, three also represent Phytophthora spp.that are quarantine regulated in the UK (P. ramorum, P. lateralis, P. kernoviae), four are Phytophthora spp.not yet reported in the UK (P. psychrophila, P. pini, P. vulcanica and P. uliginosa/P. europaea/P. flexuosa) and twelve are oomycete sequences with no match to any known species in the custom curated database or in GenBank nt (Table 5). Phylogenetic analyses placed two of these unknown sequences (sp.2 and sp.3) into the recently designated Phytophthora clade 12 (Jung et al., 2017a) (Fig. 2A) and another (sp.10) into Phytophthora clade 1 (Fig. 2B). One of the clade 12 sequences (sp.2) differs from sequences of P. versiformis and P. castanetorum by 1 bp, and from P. quercina and P.sp. ohioensis by 1 bp and a deletion (Table 5) whereas the other clade 12 sequence (sp.3) differs by 1 bp to a sequence (GenBank accession number KP208439) from an unidentified Phytophthora obtained from stream water in the Irati Forest in northern Spain (Català, Pérez-Sierra & Abad-Campos, 2015). The clade 1 sequence (sp.10) differs from any other sequence by at least seven nucleotide differences and two deletions (of 2 and 3 bp). None of the other unknown sequences could be placed within any known Phytophthora clade. Two of these sequences (sp.1 and 6) are highly similar to ‘uncultured Phytophthora’ sequences arising from metabarcoding studies carried out in forests in Spain (Català, Pérez-Sierra & Abad-Campos, 2015; Català et al., 2017). In terms of sequence identity, one sequence (sp.5) is most similar to Phytophthora boehmeriae but differs by 9% (Table 5) and one sequence (sp.7) is most similar to an ‘uncultured oomycete clone’ sequence detected in a metabarcoding study but differs by 11% (Table 5). Three of the sequences (sp.4, 9 and 11) are most similar to the downy mildew species Plasmopara halstedii but differ by 9–11% (Table 4) and two of the sequences (sp.8 and 12) are most similar to known sequences in the genus Pythium, although in all cases the sequences differ by at least 9% (Table 5).

Of the 112 soil samples that yielded Phytophthora DNA, the mean number of species detected per sample was 4.4 (SD = 2.1), with a maximum of eleven and minimum of one species detected per sample. In terms of clade associations, sequences were assigned to Phytophthora species in eight clades; only clades 5 and 9 were not represented (Table 3). Clade 8 was the most represented clade in this study with eight species detected including three well known pathogens of woody hosts (P. ramorum, P. austrocedri, P. lateralis). Of the 23 previously described Phytophthora species detected in this study, 19 produce oospores as the main survival structure and four (P. ramorum, P. chlamydospora, P. cinnamomi and P. lateralis) produce chlamydospores (Table 3).

Baiting

Five of the six soil samples collected from site 12 and stored at 4 °C for 7 days prior to baiting yielded Phytophthora and the following species were isolated; P. austrocedri (UK lineage) from two samples, P. ramorum from one sample, P. cinnamomi from three samples and P. multivora from one sample (Table 6). Following air-drying and 5 weeks storage at 4 °C, three of the six soil samples yielded Phytophthora from baits. These were P. cinnamomi from three samples and P. cryptogea and P. cambivora, each from one sample (Table 6). All six Phytophthora species found by baiting were among the twelve species detected using metabarcoding 20 months earlier (Table 3; site 12).

Methodology

We included two control reactions to test the consistency of PCR amplification of DNA present and found considerable variability in the number of Illumina reads generated per species. However, this is not unexpected given that the Ct values of the DNA dilutions included in the reactions were not consistent across each species in the mix and that P. pseudosyringae and P. ilicis sequences were clustered by Swarm into a single OTU in the first control run due to their ITS1 sequence similarity. Although ITS1 sequence length variation across Phytophthora species may theoretically affect amplification bias, a previous study carried out to determine whether ITS1 sequence length is inversely related to amplification efficiency of ten Phytophthora species showed no relationship between the two (D.E.L. Cooke, 2014, unpublished). DNA from an isolate of P. boehmeriae was included in both control mixes but P. boehmeriae sequences were not detected. Phytophthora boehmeriae is a clade 10 species with an ITS1 sequence that differs markedly (by 15% and 16%, respectively, across 192 bp) from its closest known relatives with ITS sequence data, P. kernoviae and P. morindae. Due to the presence of several base mismatches in the second round primer regions, we suggest that PCR bias, which can cause false-negative results when sequence variation at the universal primer sites prevents efficient annealing (Pawluczyk et al., 2015), was probably the most likely reason for lack of amplification of P. boehmeriae sequences.

In the first control reaction, which included P. gonapodyides in the species mix, four sequence reads were assigned using the Swarm tool to an OTU matching P. mississippiae. Phytophthora mississippiae, which has not previously been recorded in the UK, was recently described from irrigation reservoirs at a plant nursery in Mississippi and is a clade 6 species with a similar ITS sequence to P. gonapodyides (Yang, Copes & Hong, 2013). OTUs were also assigned to P. mississippiae in two soil samples in this study, both of which yielded a much higher number of sequence reads matching P. gonapodyides. Within eukaryotic genomes, rRNA (and hence the ITS1 locus) exists as multiple copies prone to sequence variation with ∼400 copies found in the P. infestans genome (Judelson & Randall, 1998). The sequencing depth afforded by Illumina technology would have allowed reads of ITS1 sequences present in low-copy number to be generated, hence it is likely that the so-called ‘P. mississipiae’ sequences are in fact being amplified from P. gonapodyides and represent a sequence variant present in low-copy number within the genome of the latter species. Similarly, sequence reads matching P. idaei and P. rubi were likely amplified from P. cactorum in the first control mix, and P. fragariae in the second control mix, respectively. It is however a possibility that the misidentifications could have occurred due to errors in PCR as the error rate of Taq polymerase is one in three million. Errors might also occur during Illumina sequencing which has an error rate of one in 1,000 bases (Schirmer et al., 2016).

Another false positive reported in the control reaction using Bowtie was P. citricola whereas an isolate of the closely related species P. plurivora was included in the control mix. P. plurivora was described in a taxonomic reclassification of P. citricola and only two bp across the ITS1 discriminate the two taxa (Jung & Burgess, 2009). Minor sequence variation or database errors may result in misclassification and some sequences of P. plurivora uploaded to GenBank before 2009 are inevitably misidentified as P. citricola.

Five species (P. cactorum, P. gonapodyides, P. ilicis, P. siskiyouensis and P. syringae) included in the control mixes were not identified using the Bowtie tool that requires an exact sequence match to the reference sequence. All five have variable ITS1 sequences reported in GenBank, but only one or two sequences per species were included in the version (June 2017) of the reference database used here. This is an important element for discussion not yet raised in the Phytophthora metabarcoding literature. Although the ITS region has proven valuable in Phytophthora diagnostics and can be used to distinguish the majority of known species, intraspecific and intra-genomic sequence variation is an important consideration for accurate identification. The reference database should include sequences representative of the known variation within species so that the assignment of sample sequences to OTUs, and species, can be based on appropriate similarity thresholds such as the one to two mismatch threshold used by Swarm here. Bowtie was not the preferred tool for the purpose of this study as it requires a perfect read match and does not allow for the fact that ITS1 sequences of species present in environmental samples might not perfectly match the reference database sequences. However, Bowtie was nonetheless useful for elucidating variation in ITS1 sequence reads. Also, where sequences were assigned using Swarm to multiple closely related species, the output from Bowtie could be informative if it had assigned the same sequences to only one species.

Our control results emphasise the need for further validation and refinement of the reference database and the metapy pipeline. Illumina sequencing of the ITS1 region amplified from pure DNA of known species should be used to identify intraspecific ITS1 sequence variants for subsequent inclusion in the reference database. Such information could also be used to test the efficiency of amplification and clustering of different sequence types. Given that the ITS1 region cannot be used to differentiate a number of very closely related Phytophthora species, as reported in this and previous studies (i.e. Català, Pérez-Sierra & Abad-Campos, 2015; Burgess et al., 2017), it would be useful to explore alternative barcoding regions for Phytophthora, such as the cox mitochondrial genes or a single copy nuclear gene that is conserved within species and yet exhibits sufficient interspecific variation to be able to discriminate species. Current advances in long read sequencing technology such as Pacific Biosciences (PacBio) or Oxford Nanopore technology also offer the prospect of increasing the diagnostic potential of new longer gene regions. For example Redondo et al. (2018b) used PacBio sequencing in a metabarcoding study of Phytophthora communities along climatic gradients in Sweden and achieved a high resolution for species identification as a result of being able to sequence both the ITS1 and ITS2 regions. However, PacBio generates fewer sequence reads than Illumina and tends to have higher error rates (Tedersoo, Tooming-Klunderud & Anslan, 2018).

Phytophthora species detected

Taking into account the observed limitations of using ITS1 as the barcode region, we nonetheless consider our approach to have been successful in that the DNA of 23 distinct Phytophthora species was detected in soil samples across the 14 sites sampled. Overall, the diversity of species detected was greater than predicted with some species far more widespread than previously thought. As expected, given the nature and Phytophthora histories of the sampling locations, all Phytophthora species detected, with the exception of P. primulae, have documented associations with woody hosts and/or ornamental shrub species.

The geographical location of each site and underlying soil type had no apparent effect on Phytophthora species’ presence and abundance, which suggests that the Phytophthora communities at each site are largely comprised of introduced species reflecting the range of out planted hosts. There were, however, differences among sites, with site 8 having the most atypical species assemblage. This site was unusual in that it had in the early 2000s been subjected to a very extensive programme of landscaping and replanting, involving the importation of high volumes of containerised plants and soil from continental Europe. The seven most abundant Phytophthora species found in this study were detected at the majority of sites, and in soils associated with a range of tree/shrub genera. All seven species, with the exception of P. austrocedri, have broad host ranges and three of these species, P. pseudosyringae, P. austrocedri and P. ramorum, are recently invasive and having a severe impact on commercial forest and natural woodland landscapes in the UK. One surprising outcome of this study was that none of the Phytophthora species detected in the soil samples were more associated with diseased than healthy hosts. Given that the majority of soils sampled here had been subject to a high turnover of plants, soil and composts over many years, including understorey plantings of annual or perennial herbaceous species, we can assume that the Phytophthora assemblages of any one sample reflect the type and source of material brought in to each site, rather than being directly associated with the host tree/shrub/stump central to each sampling point. These findings support those of Redondo et al. (2018a) who analysed Phytophthora diversity along a gradient of human interference in Sweden and concluded that invasive Phytophthora spp. share a common introduction pathway via out-planting of infected nursery stock in urban environments or other sites of high human interference. Population reservoirs are established at these locations with the more ecologically well adapted species able to spread onwards into less populated areas. Additionally, some of the more abundant species found by Redondo et al. (2018a, 2018b) in Sweden were also found in abundance at our Scottish sites; these included P. cactorum, P. cambivora, P. gonapodyides, P. pseudosyringae, P. ramorum and P. syringae.

Phytophthora pseudosyringae, the most abundant species detected in this study, is a caducous homothallic species capable of aerial dispersal and was first described in 2003 associated with Fagaceae and Alnus hosts in continental Europe (Jung et al., 2003). It is thought to be native to continental Europe, having a scattered but widespread distribution there (Jung et al., 2016). Earlier findings of this species may have been misidentified due to morphological similarities with P. syringae (hence the name ‘pseudosyringae’) (Jung et al., 2003). Following its formal description, P. pseudosyringae was isolated from coastal forests in California and Oregon, associated with the same hosts as P. ramorum (Wickland, Jensen & Rizzo, 2008). In the US the pathogen is genetically almost clonal and is hypothesised to be an introduction from Europe (Linzer et al., 2009). Phytophthora pseudosyringae was first reported in Britain in 2009 causing disease on Nothofagus spp. (Scanu, Jones & Webber, 2012) and Vaccinium myrtillus (Beales et al., 2009), and has subsequently been isolated from F. sylvatica (Scanu & Webber, 2016) and from aerial stem cankers on Larix kaempferi in Britain, often together with P. ramorum (J. Webber and A. Harris, personal communication). More recently, P. pseudosyringae was isolated for the first time from basal cankers on Aesculus hippocastanum in Sweden (Redonodo et al., 2016). Thus the geographical distribution and potential host range of P. pseudosyringae appears to be much broader than initially thought and it clearly represents a threat to numerous tree species in Britain. Although this pathogen was first reported in Britain in 2009, its widespread distribution in the country suggests that it has been here for some time (Scanu & Webber, 2016). Isolates collected in Britain and morphologically identified as P. pseudosyringae had ITS sequences that shared 100% identity with the ITS sequence of the type strain isolated in Germany in 1997 (Scanu & Webber, 2016). Analyses of the genetic diversity of British isolates across a range of loci will help to shed light on the question of whether P. pseudosyringae is native or introduced to this country.

Phytophthora austrocedri, the second most abundant species detected in this study, is causing widespread damage to Juniperus communis in northern Britain (Green et al., 2015) with all British isolates collected so far conforming to a single clonal ‘UK’ lineage (Henricot et al., 2017). The UK lineage differs genetically and morphologically from the clonal ARG lineage which has expanded in southern Argentina to cause extensive mortality of Austrocedrus chilensis (Vélez et al., 2013). To date the known host range of P. austrocedri is limited to species within the Cupressaceae, on which it is frequently intercepted in trade (J. Barbrook, Animal and Plant Health Agency, York, England, personal communication; A. Schlenzig, Science and Advice for Scottish Agriculture, Edinburgh, Scotland, personal communication). In this study, the majority of soil samples yielding P. austrocedri DNA barcode sequences and the two soil samples from which P. austrocedri was baited were collected from around tree/shrub species not documented as hosts for the pathogen. Additionally, ten of the sites represent first-time detections of P. austrocedri. If the DNA originates from live propagules, as suggested by the successful baiting of P. austrocedri from one of the sites where this pathogen has not previously been reported, then this suggests either that P. austrocedri is well adapted to survive in soil in the absence of infected host material, or it can live in association with a much broader range of plant species than is currently thought. Crone et al. (2013a, 2013b) demonstrated the persistence of resting structures of P. cinnamomi in the roots of asymptomatic non-hosts and it is quite possible that other Phytophthora species could do the same.

The third most abundant Phytophthora found in this study was P. gonapodyides, a ubiquitous clade 6 species which flourishes in aquatic habitats and is thought to play a role in breakdown of plant debris (Brasier et al., 2003). This species has been isolated from lesions on diseased trees of various species in Britain (Forest Research Tree Health Diagnostic and Advisory Service) although its role as a pathogen in these cases is unclear. Other abundant Phytophthora spp. detected here include: P. cactorum, P. cambivora and P. syringae, all of which are considered common in Britain, causing disease on a wide range of woody and non-woody hosts (Cooke, 2015); P. chlamydospora, P. cryptogea and P. plurivora, also regarded as cosmopolitan with broad host ranges (http://forestphytophthoras.org/) (Cooke, 2015), and P. ramorum, a quarantine-regulated pathogen causing serious mortality of Larix spp. in Britain and with a host range encompassing many shrubby species (Brasier & Webber, 2010).

Of the remaining Phytophthora species detected in this study, four (P. psychrophila, P. pini, P. vulcanica and P. uliginosa/P. europaea/P. flexuosa) have not previously been reported in the UK. Phytophthora psychrophila, found in highest abundance in this study around a dead C. lawsoniana at site 8, was first described in 2002 following its isolation from soil in European oak forests (Jung et al., 2002). It is considered to be a weak pathogen able to cause necrotic lesions when artificially inoculated on to roots of young Quercus robur (Jung et al., 2002). A subsequent study by Pérez-Sierra et al. (2013) however, found P. psychrophila to be partly responsible for root losses and bark cankers in Quercus forests in eastern Spain, suggesting that this species has the potential to be more aggressive than initially thought.

Phytophthora pini is established in Europe and North America as a pathogen of plants in seven genera including Pinus and Fagus as well as ornamentals and vegetables, and it is regarded as an increasing threat to the horticultural industry (Jung & Burgess, 2009; Kong et al., 2009; Hong et al., 2011). This pathogen is very closely related to P. plurivora and P. citricola and could not be distinguished from the former species based on ITS1 alone by Burgess et al. (2017). The ITS1 sequences of P. pini and P. plurivora included in the reference database used here originate from voucher specimens (GenBank accession numbers: HQ643310 and HQ643312, respectively) and differ by three nucleotides which was sufficient to differentiate the two sequences. Phytophthora plurivora was not detected in any of the soil samples yielding sequence reads of P. pini. Therefore we believe that our finding of the P. pini DNA barcode could indeed represent the first record of P. pini in the UK.

A unique sequence identical to GenBank sequences of P. vulcanica (MF036209-213) was amplified from a soil sample collected from around a healthy Castanea sativa at a single site in this study. Phytophthora vulcanica was described very recently following its isolation from rhizosphere soil of F. sylvatica in Sicily (Jung et al., 2017a). It was weakly pathogenic to young seedlings of F. sylvatica in soil infestation studies, and at this stage is presumed to be endemic to Europe (Jung et al., 2017a).

An OTU containing 2,545 sequence reads assigned to P. uliginosa/P. europaea/P. flexuosa was reported at a single site in a soil sample collected from around a Betula pendula. Phytophthora flexuosa was described recently from Taiwan (Jung et al., 2017b) and is least likely to be the species present here. BLASTN+ analysis of the sequence against GenBank nt showed an equivalent 99% identity with the same single base difference to ITS1 sequences of both P. uliginosa and P. europaea. The question as to whether it is most likely to be P. europaea or P. uliginosa present at the site is an interesting one as both species have been found in rhizosphere soils of Quercus spp. in Europe (Jung et al., 2002), with P. europaea subsequently also isolated from Quercus forests in the USA (Balci et al., 2006). Although P. europaea was found to be only weakly pathogenic on Quercus spp., P. uliginosa is reported to be much more aggressive on these hosts (Jung et al., 2002).

Evidence from DNA barcoding suggests a living organism is present but verification of these new records for Britain should be sought by re-sampling the sites and attempting to obtain living cultures from soil or from host material. Vigilance should also be maintained at these sites for disease symptoms. The same applies to sites where the quarantine-regulated pathogens P. ramorum, P. austrocedri and P. kernoviae have been detected for the first time. Conversely, P. lateralis, a regulated pathogen of Chamaecyparis lawsoniana, was detected in far fewer soil samples than expected given that many were collected from around C. lawsoniana with symptoms of active infection at known P. lateralis outbreak sites. In the Pacific North West of the USA where P. lateralis is invasive and causing high levels of mortality of C. lawsoniana in its native range, the pathogen is rarely recovered from soil even in areas of high infestation (Hansen et al., 2000). It can however be baited from the organic fraction of infested soil, leading the authors to conclude that it persists in soil only on infected root fragments, which may not have been retained in the samples processed in the current study.

Three novel Phytophthora sequences were identified in this study. Unknown sp.2, which is closest to P. castanetorum and P. versiformis in clade 12, was amplified from two soil samples; one collected from around a rhododendron and the other from around a pieris, both located at different sites. The correct species identity of our sp.2 remains unknown at present. Sp.3 also appears to be an unknown clade 12 species found in high read numbers from around a healthy rhododendron and Sp.10, which appears to be a novel clade 1 Phytophthora, was found around a healthy Juniperus communis. Other sequences did not cluster with any known Phytophthora clade and included sequences considered to be of downy mildew or Pythium origin. The amplification of these sequences from some of the soil samples indicates that the primers used here do cross-react with other oomycetes although the extent of cross reaction was low compared with the much higher abundance of unambiguous Phytophthora sequences reported. It should also be noted that two recent phylogenies of the oomycetes have placed the downy mildew species Plasmopara halstedii, which aligns most closely with our unknown sp. 4, sp.9 and sp.11, within Phytophthora close to Phytophthora clade 1 (McCarthy & Fitzpatrick, 2017; Bourret et al., 2018) with the latter authors proposing a new clade 16 within Phytophthora in which Plasmopara halstedii sits (Bourret et al., 2018).

Clearly, one shortfall of using metabarcoding alone for Phytophthora detection is that living cultures are required for absolute verification of the presence and viability of the organism and, in the case of potential new species, for taxonomic studies and host risk analyses to be undertaken. Baiting to obtain living Phytophthora isolates was not initially intended to be part of the current study. However, a subsequent visit to one of the sites 20 months after the initial soil samples were collected afforded the opportunity to try to isolate the Phytophthora species that had been detected based on DNA data. Baiting at this site, which had never previously reported any Phytophthora diseases, yielded six of the 12 Phytophthora species identified at the baited sampling points by metabarcoding. This is not surprising given the time lag between the two sampling periods and the low sensitivity of baiting which has a tendency to produce false-negative results (Vannini et al., 2013). Similarly, Khaliq et al. (2018) isolated only six Phytophthora species by baiting soil from urban parks in Perth, Australia compared with 13 species detected in the same soil samples using metabarcoding. There may also be a sampling bias inherent to baiting which can favour dominant, fast growing Phytophthora species which easily produce zoospores under the baiting conditions and which have a host preference for the baiting material used (Scibetta et al., 2012). Certainly P. cinnamomi dominated the isolates obtained by baiting in this study. Additionally, some Phytophthora species are regarded as unculturable and so the ability to detect them must rely on molecular methods (Català et al., 2017).

Most of the Phytophthora species detected here using metabarcoding produce oospores and not chlamydospores as their main survival structure. Oospores are thick walled and resilient, enabling long-term survival in soil and other substrates without a living host (Ristaino & Gumpertz, 2000). By only sampling soil, our approach might have been biased towards homothallic species that readily produce oospores, whereas water sampling may yield different or a greater variety of species (Català, Pérez-Sierra & Abad-Campos, 2015; Redondo et al., 2018a). Research is currently underway to compare metabarcoding and baiting approaches for analyses of Phytophthora diversity in both soil and water samples from a range of sites across Britain, including adjacent disturbed and less disturbed environments. It is quite likely that future protocols will use a refined combination of methods for screening environmental samples for Phytophthora pathogens.