RESEARCHFree Access icon

Population Structure of Phytophthora plurivora on Rhododendron in Oregon Nurseries

    Affiliations
    Authors and Affiliations
    • Nicholas C. Carleson1
    • Valerie J. Fieland1
    • Carolyn F. Scagel2
    • Jerry E. Weiland2
    • Niklaus J. Grünwald2
    1. 1Department of Botany & Plant Pathology, Oregon State University, Corvallis, OR
    2. 2Horticultural Crops Research Unit, United States Department of Agriculture, Agricultural Research Service, Corvallis, OR

    Published Online:https://doi.org/10.1094/PDIS-12-18-2187-RE

    Abstract

    Phytophthora plurivora is a recently described plant pathogen, formerly recognized as P. citricola. Recent sampling of Pacific Northwest nurseries frequently encountered this pathogen, and it has been shown to be among the most damaging Phytophthora pathogens on ornamentals. We characterized the population structure of P. plurivora in a survey of four Oregon nurseries across three different counties with focus on Rhododendron hosts. Isolates were identified to the species level by Sanger sequencing and/or a PCR-RFLP assay of the internal transcribed spacer (ITS) region. We used genotyping-by-sequencing to determine genetic diversity. Variants were called de novo, resulting in 284 high-quality variants for 61 isolates after stringent filtering. Based on Fst and AMOVA, populations were moderately differentiated among nurseries. Overall, population structure suggested presence of one dominant clonal lineage in all nurseries, as well as isolates of cryptic diversity mostly found in one nursery. Within the clonal lineage, there was a broad range of sensitivity to mefenoxam and phosphorous acid. Sensitivity of the two fungicides was correlated. P. plurivora was previously assumed to spread clonally, and the low genotypic diversity observed within and among isolates corroborated this hypothesis. The broad range of fungicide sensitivity within the P. plurivora population found in PNW nurseries has implications for managing disease caused by this important nursery pathogen. These findings provide the first perspective into P. plurivora population structure and phenotypic plasticity in Pacific Northwest nurseries.

    Phytophthora plurivora is a recently described species that forms part of what used to be described as the Phytophthora citricola complex of species (Jung and Burgess 2009). The P. citricola complex was recently split into eight species (P. plurivora, P. pachypleura, P. pini, P. multivora, P. capensis, P. caryae, P. acerina, and P. citricola sensu stricto), and potentially includes other taxa (Bezuidenhout et al. 2010; Brazee et al. 2017; Ginetti et al. 2014; Henricot et al. 2014; Jung and Burgess 2009; Scott et al. 2009). All species formerly part of this species complex are now placed in clade 2 of the Phytophthora genus (Hong et al. 2011).

    P. plurivora is a dominant pathogen in the Pacific Northwest nursery industry located predominantly in the states of Washington and Oregon in the United States and the province of British Columbia in Canada (Parke et al. 2014). A recent intensive sampling over 4 years in four nurseries identified P. plurivora as being the most common species of Phytophthora in communities sampled in Oregon nurseries, making up around 33% compared with P. cinnamomi (26%), P. syringae (19%), and P. citrophthora (11%) among others. This study focused on sampling plant roots and foliage from healthy and symptomatic plants, potting media, used containers, water, greenhouse soil, and container yard substrates. A more recent sampling that focused solely on foliage of healthy and symptomatic rhododendron showed it to be the second most dominant Phytophthora pathogen after P. syringae (Knaus et al. 2015).

    P. plurivora is not necessarily present in other nursery environments. For example, Patel et al. (2015) did not find it in Florida nurseries from symptomatic host samples. They found predominantly P. cinnamomi and P. parasitica, which thrive in warm, tropical climates, unlike P. plurivora. Another recent study focusing on sampling Phytophthora from irrigation reservoirs in Mississippi and Alabama nurseries did not encounter P. plurivora (Copes et al. 2015). Schlenzig et al. (2015) sampled P. plurivora from Rhododendron, Pieris, and Syringa, but it was not the dominant species found in Scottish nurseries. In addition, Prospero et al. (2013) sampled Swiss nurseries and found P. plurivora to be dominant second to P. ramorum.

    P. plurivora has a wide host range and is known to infect woody plants in nurseries, as well as forest stands. Woody ornamentals commonly affected include Rhododendron, Pieris, and Syringa (Parke et al. 2014; Prospero et al. 2013; Weiland et al. 2018). Forest stands and landscape specimens of beech, maple, oak, ash, and spruce have been increasingly impacted since the discovery of the pathogen (Jung and Burgess 2009; Mrázková et al. 2013; Parke et al. 2014; Rytkönen et al. 2013; Schlenzig et al. 2015; Weiland et al. 2010). The decline of Fagus sylvatica (European beech) has been associated with P. plurivora and climate extremes (Jung 2009). Weiland et al. (2010) showed isolates from F. sylvatica were equally aggressive on lilac. Isolates of P. plurivora from the Pacific Northwest have been shown to cause root rot on Rhododendron (Weiland et al. 2018). Free-living P. plurivora has also been retrieved from soil and streams (Huai et al. 2013).

    A homothallic mating system was previously reported for P. plurivora (Jung and Burgess 2009). Homothallic species are able to produce sexual spores through self-fertilization and the presence of only one mating type (Ni et al. 2011). Homothallism is more common in Phytophthora than heterothallism, although the mating system is a mosaic across the genus and is suspected to have evolved several times during species radiation (Kroon et al. 2004). It has been hypothesized that self-fertilization reduces genetic diversity, driving heterozygosity down over time (Charlesworth 2003). Indeed, a self-fertilizing diploid species would theoretically lose 50% of its heterozygosity in each generation (Goodwin 1997).

    Previous genotyping efforts on P. plurivora have shown that rates of heterozygosity are similar to other homothallic Phytophthora species including P. alni (subsp.) uniformis and P. sojae (Schoebel et al. 2014). However, recent genome sequencing of the homothallic P. cactorum showed that heterozygosity was high enough such that genome assembly proved challenging (Yang et al. 2018). A genome for P. plurivora was not available at the time the work of this study was conducted, and previous studies used microsatellite data to characterize population structure (Vetukuri et al. 2018). Genome-wide sequencing might more accurately capture the genotypic diversity of this prevalent plant pathogen.

    Little is known about the phylogeographic and evolutionary history of P. plurivora. A recent study investigated the global population structure of P. plurivora (Schoebel et al. 2014). This study included 14 strains from streams, ponds, and nurseries in Oregon. It found that the pathogen is most likely native to Europe and was spread around the world through nursery exports of diseased plant material. Additionally, coalescent analyses showed the pathogen population size is expanding along with its geographic range. However, population structure has yet to be understood; genetic clustering did not correlate to geographic groups of isolates. Furthermore, the null hypothesis of random mating, as tested using the index of association r¯d (Agapow and Burt 2001), could only be rejected in six geographic populations. The West Coast region of the United States, including the Oregon populations, had the highest r¯d value in the study, indicating an almost entirely clonally reproducing population (Schoebel et al. 2014).

    Genotyping-by-sequencing (GBS) is a state-of-the-art technique for genotyping based on restriction enzyme digests and short read sequencing (Ellegren 2014). This technique, and other methods for sequencing a reduced representation of the genome, yields a high number of single nucleotide polymorphisms (SNPs) at a moderate cost (McCormack et al. 2013). SNPs discovered using reduced representation methods have been applied in plant pathology and microbial ecology for describing population structure, evolutionary history, and linkage or association analyses (Grünwald et al. 2016). Recent studies used GBS to characterize populations of the broad-host-range pathogen Verticillium dahliae (Milgroom et al. 2014), the potato late blight pathogen Phytophthora infestans (Hansen et al. 2016), and the raspberry pathogen Phytophthora rubi on the U.S. West Coast (Tabima et al. 2018).

    Common fungicides for Phytophthora disease management on trees or other woody hosts include mefenoxam and various phosphonates such as fosetyl-Al (Keinath 2007; Matheron and Porchas 2000; Weiland et al. 2009). A topical application of mefenoxam was shown to inhibit canker expansion in almond trees inoculated with P. citricola s. l. when treated postinfection with low variability in host sensitivity (Browne and Viveros 2005). This study, and others on P. citricola conducted before it was split into various Phytophthora species, likely include what we now know as P. plurivora as well as other species. The only study to date focusing on the recently described P. plurivora found phosphorous acid to be an effective fungicide at a concentration of 34 µg/ml reducing in vitro radial mycelial growth by 50% in the one isolate tested (Dalio et al. 2014). In Weiland et al. (2009), one isolate previously identified as P. citricola, but now identified as P. plurivora (isolate 37-04a), was sensitive to mefenoxam at 3.0 µg/ml and phosphorous acid at 47 µg/ml.

    We isolated 69 strains from Rhododendron hosts in four Oregon nurseries to assess fungicide resistance and genetic diversity of P. plurivora populations. Fungicide testing of mefenoxam and phosphorous acid were carried out on a random population sample. We used GBS to gain insights into genetic diversity, and to determine mode of reproduction and differentiation of populations across nurseries. A secondary objective was to optimize variant filtering procedures to minimize the effects of genotyping error to produce a highly accurate representation of population structure. This work has implications for transportation and sale of nursery stocks, chemical management of disease, and for designing GBS experiments for Phytophthora pathogens.

    Materials and Methods

    Population sampling and species identification.

    Sampling of P. plurivora is complicated by the fact that several Phytophthora species can be isolated from the same host (Knaus et al. 2015; Parke et al. 2014). Thus, extensive sampling of nontarget species was required to obtain a representative population level sample of P. plurivora (Grünwald et al. 2017). Isolates were collected over several months across several years (2007–2015) until we obtained a large enough sample size for analysis. We sampled from four Oregon nurseries (14–21 strains per nursery; nurseries are labeled A, B, C, and H) from different plants using a W transect pattern every 5–10 m in different blocks of host plants focusing on rhododendron with symptoms of root rot (chlorosis, stunting, and wilting). The four nurseries sampled were less than 150 km apart. Although we tried to sample hierarchically to obtain subpopulations in several blocks of plants per nursery, we did not succeed. Thus, our isolates are combined by nursery and only diversity across nurseries can be assessed (Table 1). P. plurivora was isolated by plating necrotic tissues from the root, the root collar region, or leaves according to previous methods (Knaus et al. 2015; Weiland et al. 2018). At nursery H, some isolates were obtained from the soil (Parke et al. 2014).

    Table 1. Populations and samples of Phytophthora plurivora sampled in Oregon nurseries from Rhododendron leaves, roots, or potting soil. N is the number of individuals sampled from each nursery, and NF denotes how many of those were retained after stringent filtering of the GBS data. The following diversity statistics were calculated based on this size. MLG is the number of multilocus genotypes observed, allowing one allelic difference between isolates. eMLG is the number of expected MLGs at a sample size of 4 based on rarefaction. The standard error of this statistic is given in the next column, SE. The next columns are the following diversity statistics considering the MLG diversity in each nursery: Shannon-Wiener Index (H), Stoddart and Taylor’s Index (G), Simpson’s index (lambda), Simpson’s Index after correcting for differing sample sizes (lambda*), evenness (E5), and Nei’s unbiased gene diversity (Hexp). The standardized index of association r¯d was calculated on clone-corrected data and the adjacent P values were calculated with a null distribution created from 999 permutations of the data.

    P. plurivora isolates were identified to species using a PCR-RFLP assay as described previously (Knaus et al. 2015). Briefly, the restriction assay uses polymerase chain reaction amplicons of the internal transcribed spacer (ITS) region digested with the enzyme AciI and resolved on a 2% agarose gel. P. plurivora could be unambiguously identified to species based on reference strains for which we sequenced both the ITS and cox spacer regions (Grünwald et al. 2011). Isolates with a collector other than N. Grünwald were identified by ITS sequence identity (Supplementary Table S1).

    DNA extraction and genotyping by sequencing (GBS).

    All samples (isolates and technical replicates) were placed on a 96-well plate and lyophilized. Whole DNA was extracted from 69 unique isolates and 27 technical replicates using a KingFisher Flex magnetic particle-processing robot (Thermo Scientific, Waltham, MA) at the Center for Genome Research and Biocomputing (CGRB) at Oregon State University using the lyophilized tissue. DNA was extracted using the Omega MagBind plant DNA kit (M1130, Omega Bio-tek, Norcross, GA). The technical replicates were used in downstream analyses to adjust stringency of variant filtering and to minimize genotyping error. We employed a two-enzyme combination GBS protocol that uses a rare-cutting enzyme and a common-cutting enzyme to greatly reduce genome complexity (Poland et al. 2012). Genomic DNA samples were digested with a combination of the PstI (CTGCA/G) and MspI (C/CGG) restriction enzymes as well as an SbfI (CCTGCA/GG) and MspI combination (New England Biolabs, Ipswich, MA). The details for restriction, ligation, and library preparation are deposited in two excel spreadsheets per enzyme combination at OSF (Carleson and Grünwald 2018, https://github.com/grunwaldlab/PhPlurivora_OR_Nurseries_GBS). After barcoding with Illumina adapters, the equimolar, pooled library was sequenced on the Illumina HiSeq 3000 (San Diego, CA) with 150-bp paired-end reads.

    Read processing and SNP calling.

    Each of the experiments were demultiplexed using Stacks version 1.47 with default parameters (process_shortreads.pl) (Catchen et al. 2011, 2013). All samples were genotyped and filtered including all technical replicates. For the purposes of this study, a technical replicate is either DNA extracted from the same culture and placed in two wells, or DNA from one isolated individual that was cultured on different plates. If two replicates of the same isolate were retained after filtering, the replicate with the most reads retained after demultiplexing was selected. We expect that isolates from the same plant are likely to result in sampling the same genotype; thus, three such isolates were excluded from analysis to eliminate bias. In all allelic calculations, missing values (NAs) were removed.

    SNP calling was performed de novo using Stacks version 1.47 (denovo_map.pl) (Catchen et al. 2011, 2013). A reference genome for P. plurivora was not available at the time work was completed (Vetukuri et al. 2018). Default parameters were used for SNP calling, except for the following parameters: minimum coverage (m = 5), maximum mismatches between a single individual (M = 5), and maximum mismatches during catalog building (n = 6). These parameters were selected from recommendations in Paris et al. (2017). All isolates were assumed to be diploid.

    Variants were stringently filtered in R version 3.4.1 using the vcfR package (R Core Team 2017; Knaus and Grünwald 2017). The 31 technical replicates sequenced were matched up to their respective samples to determine that quality control methods were producing accurate genotype information. Allele calls were first censored based on quality score, censoring those with very high or very low read depth. Any allele censored during quality control is equivalent to a missing genotype of “./.” observed directly in a VCF file. Samples and variants were then removed based on percentage of data missing. Any variants that had a minor allele count of 1 (minor allele frequency = 0.011) were then removed. Next, variants with any missing data were either omitted or not to test the effects of missing data on population structure. After this last step of variant omission, all monomorphic variants were removed. Data were transformed into genlight and genind objects usding vcfR (vcfR2genlight, vcfR2genind). All scripts and functions for read processing and SNP filtering are deposited at OSF (Carleson and Grünwald 2018, https://github.com/grunwaldlab/PhPlurivora_OR_Nurseries_GBS). These scripts are annotated to provide more in-depth information for how quality-control was performed.

    Fungicide resistance.

    A select group of P. plurivora isolates (nurseries A, B, and C only) were tested for fungicide resistance to mefenoxam (n = 21) and phosphorous acid (n = 19). Methods were modified from those described in Weiland et al. (2009). Each isolate × fungicide pairing was tested in triplicate in two separate trials to provide experimental repetition. Plugs of inocula (5-mm-diameter) were cut from the edges of cultures and transferred to plates containing either mefenoxam (Subdue MAXX, 22% active ingredient, Syngenta Crop Protection, Greensboro, NC; concentrations tested: 0, 0.1, 10, 100, and 1,000 µg a.i./ml) or phosphorous acid (mono- and di-potassium salts of phosphorous acid, Alude, 45.8% a.i., Cleary Chemical Corp., Dayton, NJ; concentrations tested: 0, 100, 400, 700, 1,000, and 5,000 µg a.i./ml). Plates were incubated at 22 ± 1°C for 5 days. The radial growth was then calculated as the average of two measurements along two perpendicular axes, and 5 mm was subtracted from the average to account for the initial plug size.

    Statistical and population genetic analysis.

    Population genetic analysis was conducted in R using various packages (R Core Team 2017). The R package poppr (Kamvar et al. 2014) was used to convert from genind to genclone objects, assign strata and populations, inspect basic statistical properties of the data, calculate index of association (r¯d), and create minimum spanning networks as calculated with Prevosti’s distance. This distance is the fraction of allelic differences between two isolates (Prevosti et al. 1975). It was also used to estimate diversity (expected number of multilocus genotypes eMLG, Shannon-Wiener index H, Stoddart and Taylor’s index G, and Simpson’s index lambda, evenness E5, and gene diversity Hexp) (Grünwald et al. 2003). eMLG estimates the expected number of MLGs in a population based on rarefaction down to the number of isolates in the least-sampled population. H and G both measure richness and adjust for evenness, increasing as the number of genotypes increase with G tending to be larger than H, and lambda estimates the probability (from 0 to 1) that any two genotypes randomly sampled from a population are different. In a population where all genotypes are present in equal abundance to each other, E5 will be closer to one, while a population with one dominant genotype and several different genotypes in very low abundances will be closer to zero. Lastly, Hexp, or expected heterozygosity, expresses the overall probability (from 0 to 1) that any pair of alleles, randomly sampled from the population, will be different.

    The R package adegenet (Jombart 2008; Jombart and Ahmed 2011) was used to conduct a Discriminant Analysis of Principal Components (DAPC), optimize the a score for optimal numbers of principal components, and conduct principal component analysis (PCA). The R package ade4 (Dray and Dufour 2007) was used to visualize the DAPC in a scatter plot. To calculate Fst and other fixation indices between isolates, the R package strataG (Archer et al. 2017) was used. The R package ggplot2 (Wickham 2009, 2016) was used to visualize r¯d simulations, PCA, and fixation indices. The R package igraph (Csardi and Nepusz 2006) was used to visualize the minimum spanning network using the GEM force-directed layout algorithm (Frick et al. 1994). All scripts and functions for data analysis can be found on OSF (Carleson and Grünwald 2018).

    Fungicide sensitivity to mefenoxam or phosphorous acid was calculated as EC50 values from plots of the culture diameters at each fungicide concentration (Saville et al. 2014). The R package drc (Ritz et al. 2015) was used for generating dose-response curves and determining EC50 values using a log-logistic model for each isolate × fungicide combination. A three-parameter model was used if an isolate showed no growth at the highest concentration of a given fungicide (with the assumption that the lower limits of the data are 0) and a four-parameter model was used if an isolate exhibited growth at the highest concentration. Mean EC50 were calculated across three replicates in two trials, and frequency data were plotted to show variability in mean fungicide sensitivity. A Spearman’s rank correlation coefficient (R) was calculated with the R base package to assess the relationship between sensitivity to mefenoxam and sensitivity to phosphorous acid for isolates with EC50 values for both fungicides (n = 19).

    Results

    Genotyping-by-Sequencing data.

    We sequenced 96 samples, including 69 unique isolates and 27 technical replicates, with PstI-MspI and SbfI-MspI restricted libraries on the Illumina HiSeq 3000 which resulted in 155,083 raw SNPs and 81.66% missing data. Since the SbfI-MspI sequencing resulted in a lower mean read count, all analyses were performed using sequence data from the PstI-MspI restricted library. No reads were removed for low quality during demultiplexing prior to SNP calling. After SNP calling, data were stringently filtered based on both hard cut-offs and percentile thresholds for abnormally low or high sequence coverage at variant loci to reduce bias that can inflate diversity, or cause isolates to cluster together solely based on data quality. Both samples and variants were filtered with increasing stringency until all remaining technical replicates for each isolate showed 0 pairwise genetic distance. Using sequence coverage, missingness, and minor allele frequency as quality control metrics, we filtered our dataset down to 284 high confidence SNPs in 87 samples. Of the 87 samples, 61 were unique isolates, corresponding to 23 multilocus genotypes (Table 1, Supplementary Table S2). The mean read depth of every sample was above 30×, and the overall mean was 57.2× coverage (Supplementary Fig. S1).

    Variation in genetic diversity of P. plurivora in nursery populations.

    Most nurseries had low genotypic diversity where roughly one-fourth to one-half of the individuals belonged to the same clone (Table 1). Nursery H had the highest genotypic diversity, while nursery A had the largest sample size and the most MLGs. The Shannon-Wiener index, Stoddart and Taylor’s index, and Simpson’s index were all highest in nursery H, followed by nurseries A, C, and then B. Evenness followed the same pattern as genetic diversity. Furthermore, nurseries B and C appeared to have similarly low diversity statistics compared with nurseries A and H. This apparently paired behavior was most pronounced for the Shannon-Wiener index. Stoddart and Taylor’s index showed the greatest differences between nurseries, while Simpson’s index, with or without accounting for different sample sizes directly, showed a more even distribution of diversity among the four nurseries. Nursery H was the most even, followed by nurseries A, C, and B.

    To test the hypothesis that P. plurivora is reproducing clonally, we estimated linkage disequilibrium by calculating the index of association r¯d. It was significantly different from 0 in nurseries B, H, and C (Table 1). The r¯d for nursery A was 0.023. This value is close to 0, indicating that we cannot reject the null hypothesis of linkage equilibrium (P > 0.2). In nurseries B, H, and C, as well as overall, we can reject the null hypothesis of linkage equilibrium genome-wide (P < 0.05), suggesting high rates of selfing and clonal reproduction.

    Population structure among P. plurivora nursery populations.

    The structure of the P. plurivora populations was dominated by the presence of one clonal lineage. A total of 52 isolates were clustered nearly, or on top of each other in a principal component analysis (PCA) while nine isolates did not fall into this lineage (Fig. 1). Those few isolates that appeared distinct from the clade containing the clonal lineage were also highly genetically distinct from each other. Despite the high genetic distance, these isolates are genetically related and were considered to cluster together into a clade of cryptic diversity (Supplementary Fig. S2). Even though there were several isolates that were genetically distinct from the clonal lineage, there were some isolates that were exact clones, indistinguishable from technical replicates (H-1 and H-2, H-8 and H-15). These were grouped into the same MLG at a threshold of a Prevosti’s distance of less than or equal to 0.01 using the farthest neighbor algorithm.

    Fig. 1.

    Fig. 1. Principal component analysis of 61 Phytophthora plurivora isolates in four Oregon nurseries. Principal components were calculated using 284 SNPs. Sample labels that are distant from their corresponding point are connected by lines. The multicolored spot on the top left (labeled ‘Clone, n = 52’) represents most of the samples falling into one clonal lineage. Only samples that were not a part of this proposed dominant clonal lineage are labeled. Each axis is labeled with the percent of total variance explained by that dimension.

    Download as PowerPoint

    Samples from all four nurseries had members belonging to the dominating clonal lineage. The isolates of marked diversity were mostly from nursery H, from which half the collected isolates were diverse and the other half belonged to the clonal lineage. However, there was one isolate from nursery C that was highly distant from the rest of the isolates collected from the nursery (C-8), sharing between 20 and 22 common alleles out of the 284 SNPs. In contrast, the rest of nursery C shared between 282 and 284 common alleles. C-8 appeared to be most closely related to H-3, sharing 281 common alleles. But this pairwise difference was not small enough to place them into the same MLG according to the threshold based on Prevosti’s distance specified above.

    Despite the high similarity between members within the dominant clonal lineage, MLGs observed within the clone clustered together by nursery (Fig. 2). This was particularly the case for nurseries A, B, and C. Isolates from nurseries H and C were shared in two MLGs, and one of those MLGs also contained isolates from nursery B. All technical replicates clustered into the same MLG.

    Fig. 2.

    Fig. 2. Minimum spanning network of 23 Phytophthora plurivora multilocus genotypes from four Oregon nurseries with 284 SNPs. The dominant clonal lineage is shown connected with large edges with very low genetic distance. Each multilocus genotype is displayed as a circle, with number of isolates sharing that MLG used as the scale. Colors of the circles show the nursery an isolate was collected from, and the ratio of colors within each circle is proportional to the ratio of isolates from each nursery. Edge width and grayscale are calculated using Prevosti’s distance. Edges representing a distance of 0.35 or greater were removed.

    Download as PowerPoint

    Our analysis of molecular variance (AMOVA) revealed that 29% of the genetic variance was observed between nurseries (P = 0.004) and 65% within nurseries (P = 0.099) (Table 2). Fst was highest between nursery A and all other nurseries (ranging from 0.181 to 0.584), followed by nursery H (0.129–0.584) (Fig. 3). Overall, differentiation ranged from small (H versus C: 0.129) to very high (H versus A: 0.584). However, this Fst calculation was limited in power due to the high level of homozygosity. Across the 284 SNPs 26 clone-corrected individuals, there are a total of 7,384 genotyped loci. Only 36 of the loci were heterozygous. The other 7,358 loci were homozygous for either the major or minor allele.

    Table 2. Analysis of molecular variance (AMOVA) for clone-corrected nursery populations of Phytophthora plurivora. Significance of variance was tested from 999 permutations of the data.

    Fig. 3.

    Fig. 3. Pairwise Fst values between clone-corrected Phytophthora plurivora populations from four Oregon nurseries (A, B, C, H). Values were averaged from 284 SNPs in the 23 multilocus genotypes, stratified into 26 MLGs after accounting for MLGs shared between nurseries.

    Download as PowerPoint

    In previous, unpublished work, we observed a segregating site within the P. plurivora strains for the ITS region with either a T or Y (ambiguous C/T) at position 396 in our alignment. This ITS polymorphism did not correlate with either population structure or phylogenetic clustering. The collection source (leaf, root, stem, or soil) also did not appear to be an artifact driving the split in population structure between clone and cryptic diversity. The isolates collected from soil were from nursery H, but isolates from soil were observed in both the dominant clonal lineage and the clade of cryptic diversity.

    Within the cluster of MLGs that are highly genetically related, there was a broad range of phenotypic response to fungicide ranging from sensitive (radial growth EC50 < 1 µg/ml) to more resistant (radial growth EC50 > 100 µg/ml) for both mefenoxam and phosphorous acid (Fig. 4, Supplementary Table S3). We observed a significant, positive relationship between sensitivity of the two fungicides (P < 0.05), indicating that an isolate with high sensitivity to one fungicide was likely to also be highly sensitive to the other fungicide tested (Fig. 5). One outlier isolate was excluded from this analysis because its mean EC50 of phosphorous acid (716 µg/ml) was more than fourfold larger than that of the next closest isolate (161 µg/ml). All isolates with EC50 values >100 µg/ml of mefenoxam were isolated from nursery A, whereas isolates with EC50 values >100 µg/ml of phosphorous acid were found at each of the three nurseries tested.

    Fig. 4.

    Fig. 4. Distribution of effective fungicide concentration that inhibited radial growth of Phytophthora plurivora by 50% (EC50). A, mefenoxam (n = 21). B, phosphorous acid (n = 19). Mean EC50 were calculated across three replicates in two trials and plotted on a log10 scale. Isolates were tested from nurseries A (red), B (yellow), and C (green).

    Download as PowerPoint
    Fig. 5.

    Fig. 5. Spearman rank correlation of effective fungicide concentration that inhibited mean radial growth of Phytophthora plurivora (n = 19) by 50% (EC50) for the two fungicides mefenoxam and phosphorous acid. Isolates are shaped according to the nursery from which they were collected. Mean EC50 were calculated across three replicates in two trials and plotted with log10-adjusted axes. Two isolates with radial growth measurements for only one fungicide were excluded from the analysis.

    Download as PowerPoint

    Discussion

    We studied populations of the common nursery pathogen P. plurivora in the Pacific Northwest. We sampled populations in four nurseries in the state of Oregon and characterized genetic diversity using genotyping-by-sequencing. Overall, populations were mostly clonal, of low diversity, and with most variance associated with differences within nurseries. Two nurseries were more diverse than the other two nurseries. However, the most salient result is the observation of one clade of a dominant clonal lineage as well as a clade of cryptic diversity mostly observed in one nursery.

    Cryptic diversity did not appear to be entirely a consequence of nursery sampling. Members of the diverse clade were mainly from nursery H, but nursery C had one isolate clustering with an MLG from nursery H in the cryptic diversity clade. Furthermore, nursery H was split in half, with seven isolates of cryptic diversity and seven isolates belonging to the dominant clonal lineage. Despite the presence of diverse isolates in nurseries H and C, linkage disequilibrium as calculated with the index of association was high, suggesting that these populations are highly selfing or clonal (Grünwald and Hoheisel 2006; Halkett et al. 2005; Kamvar et al. 2014). Thus, we propose that the observed cryptic diversity is an artifact of sampling bias. Members of the cryptically diverse clade likely belong to their own clonal lineage, or lineages, that were not sampled, or not sampled deeply enough to observe membership.

    We tested how different individuals within the clonal lineage respond to the fungicides mefenoxam and phosphorous acid. While genetic similarity is high within the clonal lineage, we observed a broad range of phenotypic response to both fungicides. This indicates that resistance is likely based on a complex genetic structure with many loci contributing to the quantitative trait, rather than a simple system of an individual being either highly sensitive or highly resistant. The substantial differences between isolates indicate resistance is not likely to have been inherited from a parent population. If the ancestral source population for the outbreaks in Oregon nurseries was resistant to mefenoxam or phosphorous acid, the trait would likely be more common within the clonal lineage. Given low differentiation within the clonal lineage, it is more likely that fungicide resistance evolved recently. Fungicide sensitivity has been found to decrease in Phytophthora spp. after continued fungicide use (Gisi and Cohen 1996; Grünwald et al. 2006; Parra and Ristaino 2001). A range of resistance in P. cinnamomi to a fungicide was shown to have developed within a single clonal lineage (Dobrowolski et al. 2008). Deeper sampling and resistance screening of both nursery and wild populations in the region is needed to determine if chemical practices influenced the evolution of resistance in multiple nurseries. Genetic structure of resistance may be resolved by deeper sequencing, revealing how resistance evolves and if it evolved multiple times in different P. plurivora populations.

    A previous phylogeographic analysis of global P. plurivora populations fully supports our results (Schoebel et al. 2014). Using STRUCTURE, Schoebel et al. inferred the most likely number of groups to be K = 4. In that study, individuals belonging to two of the four groups found globally were observed in Oregon, and for the isolates that were genotyped in both studies, there were no cases in which samples belonged to the same group in Schoebel et al. (2014) but to different clades in the present study. In other words, Schoebel et al. (2014) inferred similar population structure in Oregon nurseries using k-means clustering as the present study did with a larger sample size focusing on Oregon. The present study provides greater resolution of these two groups, showing diversity among clones. It also confirms the presence of a single polymorphic site in the ITS gene, though this site did not correlate with any demographic or population structure.

    Our results confirm previous studies that P. plurivora is a homothallic oomycete (Jung and Burgess 2009). The standardized index of association rejected the null hypothesis of linkage equilibrium, e.g., random mating, in three out of four nurseries. Nursery A is the exception and is characterized by several MLGs that are a low distance apart. These MLGs were collapsed in other nurseries, but enough difference remained in nursery A to maintain diversity. Many of these individuals that are slightly differentiated from the dominant genotype were sequenced with technical replicates where genetic distance equals 0 among those replicates. This adds support to the hypothesis that the observed diversity is real rather than an artifact of genotyping error. This slight genetic differentiation could be leaving a signature that causes r¯d to not reject the null hypothesis of random mating, despite the organism truly not reproducing with random mating. Failure to reject the null hypothesis of random mating does not equate to confirmation of sexual recombination.

    Additional support for clonality and selfing is provided through the level of heterozygosity. Heterozygosity was very low; less than 5% of the observed, clone-corrected polymorphic loci were heterozygous. This extreme homozygosity has been reported in P. plurivora before and serves as further evidence of severe inbreeding, since theoretically, an organism that selfs will lose half of its heterozygosity each generation (Goodwin 1997; Schoebel et al. 2014). With only a few heterozygous positions in each individual, and low sample sizes after clone-correction, the high Fst values of the present study are unsurprising. High Fst values correspond with a reduction in heterozygosity compared with the expectation for populations under Hardy-Weinberg equilibrium, and high rates of selfing may be responsible for the dramatically low heterozygosity observed.

    A component of error in studies using SNPs is genotyping error. In population genetic studies, some genotyping error is to be expected (Grünwald et al. 2016). This error can be introduced at any of the many steps in the GBS protocol, including DNA isolation, use of restriction enzymes, sequencing on the Illumina machine, and variant calling (Narum et al. 2013; Nielsen et al. 2011; Pompanon et al. 2005). To address genotyping error, we implemented several common methods of filtering the data in the R package vcfR and published custom functions online (Carleson and Grünwald 2018; Knaus and Grünwald 2017). Read depth was the initial parameter we filtered on; we had lower confidence in genotypes called at loci with lower coverage, but loci with coverage that is too high may indicate a repetitive region prone to misalignment which may lead to errors in SNP and genotype calling (Andrews et al. 2016; Nielsen et al. 2011). High missingness is another indicator of poor genotyping accuracy we used to filter our data (Neale and Purcell 2008). In our case, the final dataset had no missing data. Setting a minimum minor allele count of two across all isolates reduced bias of rare alleles in the dataset further. Another consideration is that rare alleles observed in only a single isolate might well be sequencing error. Given this process, we undoubtedly eliminated some SNPs that were not subject to genotyping error (false negatives), but enough markers were retained to test the hypotheses of interest to us.

    To validate our variant filtering approach, we used several technical replicates consisting of multiple samples of DNA isolated from the same individual and genotyped using GBS. The present study was able to leverage the power of nearly 30 technical replicates, verifying the overall topology as well as the diversity within a clone. By filtering increasingly stringently until all technical replicates had a distance of 0, we were able to minimize false positive SNP calls while maximizing the number of true positives. This approach utilized techniques to censor and filter data to produce enough SNP calls to be informative while eliminating any genetic distance within each technical replicate thus providing high confidence into our MLG grouping and the identification of a dominant clonal lineage.

    The use of de novo SNP calling presents a limit to the power of this study. Without genomic coordinates, we do not know the position of SNPs in the genome. While we discovered about 300 SNPs with high confidence, we can only use these for population genetic inferences, but not to infer differential selection on housekeeping or effector genes. Further, it does not allow us to design PCR primers amplifying regions that characterize the dominant clone. We are also not able to determine if certain SNPs fall into the same linkage block, which would affect inferences made about the mating system.

    This work provides novel insights into the population structure and phenotypic diversity of P. plurivora on the U.S. West Coast. Using SNPs inferred de novo, we found presence of one dominant clonal lineage with limited diversity and a less frequent clade of cryptic diversity limited mostly to one nursery. More work is needed to establish the importance of the cryptic clade and the evolution of fungicide resistance in P. plurivora from U.S. nurseries.

    Acknowledgments

    We thank Karan Fairchild, Meg Larsen, and Caroline Press for maintenance of the cultures and general technical support. We thank Bryan Beck for maintaining cultures and testing fungicide sensitivity. We thank Javier Tabima and Brian Knaus for advice in handling genotyping-by-sequencing data. We thank the Center for Genome Research and Biocomputing (CGRB) at Oregon State University for access to their high-performance computing facility, sequencing and bioinformatic support.

    The author(s) declare no conflict of interest.

    Literature Cited

    The author(s) declare no conflict of interest.

    Funding: This project was supported in part by funds from USDA-ARS CRIS Project 2072-22000-041-00-D, 2072-22000-043-00-D, and the USDA-ARS Floriculture Nursery Initiative.