Genotyping-by-Sequencing Reveals Fine-Scale Differentiation in Populations of Pseudoperonospora humuli

Downy mildew pathogens are members of the Peronosporaceae (Oomycota, Oomycetes, Peronosporales), one of the largest groups of fungi-like organisms that parasitize flowering plants (Palti and Kenneth 1981). Within the genus Pseudoperonospora, P. cubensis and P. humuli are the most economically important pathogens and are closely related genetically as sister taxa (Choi et al. 2005; Mitchell et al. 2011; Runge et al. 2011). P. humuli may infect multiple hosts (Hoerner 1940; Runge and Thines 2012), although the pathogen appears to be best adapted to hop (Humulus lupulus) and is found mainly in hop production regions of the northern hemisphere and in Argentina (Ojiambo et al. 2015; Royle and Kremheller 1981). In contrast, P. cubensis parasitizes numerous plants within the family Cucurbitaceae worldwide and has an uncharacteristically large host range compared with other downy mildew pathogens (Cohen et al. 2015; Ojiambo et al. 2015). Infection of hop by P. cubensis and certain cucurbits by P. humuli may occur with artificial inoculation under experimental conditions (Mitchell et al. 2011; Runge and Thines 2012). However, genetic evidence to date suggests the organisms are distinct species and differentiated by host in natural populations (Kitner et al. 2015; Ojiambo et al. 2015; Runge et al. 2011; Summers et al. 2015; Wallace and Quesada-Ocampo 2017).

P. humuli perennates as mycelium in the perennial root system of hop plants (Coley-Smith 1962, 1965; Royle and Kremheller 1981), contributing to the persistence of the pathogen within infected plants (Coley-Smith 1964) and the polyetic nature of disease outbreaks (Gent et al. 2008; Johnson and Anliker 1985). Hoerner (1940) reported that under controlled conditions, P. humuli could produce sparse sporulation on all plants tested in the genera Cannabis, Celtis, and Urtica. Runge and Thines (2012) reported that under controlled conditions, two isolates of P. humuli could cause disease on Cucumis sativus and two wild cucurbit species, Bryonia dioica and Sicyos angulatus. In contrast, Mitchell et al. (2011) found no evidence that P. humuli causes downy mildew on cucumber or cantaloupe. The frequency and relevance of natural infection of these plants to downy mildew epidemics is unknown but is probably insignificant based on extant genetic data (Mitchell et al. 2011; Runge et al. 2011; Summers et al. 2015; Thomas et al. 2017a).

P. humuli appears to be homothallic (Gent et al. 2017) and oospores are formed abundantly in inoculated plant tissues and under natural conditions (Bressman and Nichols 1933; Chee and Klein 1998; Coley-Smith 1962; Gent et al. 2017; Royle and Kremheller 1981). A few reports indicated successful germination and infection of hop by oospores of P. humuli (Bressman and Nichols 1933; Jones 1934; Magie 1942). However, most attempts to germinate oospores have failed (Gent et al. 2017; Royle and Kremheller 1981) and the role of oosporic inoculum in the development of downy mildew of hop under natural conditions is based largely on circumstantial evidence. Several authors discount the importance of oosporic inoculum of P. humuli as a source of primary inoculum (Coley-Smith 1962; Royle and Kremheller 1981; Skotland and Johnson 1983).

Little is known about the population diversity of P. humuli, in part because of the obligately parasitic lifestyle of the organism, which adds expense and complexity to experiments because of the difficulty of isolate collection, maintenance, and seemingly simple tasks such as DNA extraction. Furthermore, axenic culture is not possible, complicating experimental design and data interpretation because of the potential for contaminating DNA from the host and epiphytic organisms. Chee et al. (2006) assessed genetic diversity in P. humuli using marker systems based on random priming and heterogeneous samples of P. humuli (bulked sporangia from diseased shoots) collected from hop yards in Oregon and Washington. In this approach, the authors found a greater number of genotypes in samples collected from Oregon and suggested that sexual recombination may occur in Oregon but not in Washington. Mitchell et al. (2011) and Mancino (2013) conducted a multigene phylogenetic analysis of a large number of isolates of P. humuli and did not find evidence for population differentiation by region of origin. Wallace and Quesada-Ocampo (2017) identified microsatellite markers from transcriptome sequencing of the sister species, P. cubensis, and demonstrated transferability of the markers in P. humuli. Genotyping using 11 microsatellite markers identified in P. cubensis on 22 isolates of P. humuli originating from mostly the United States, but also two isolates from Europe and Japan, suggested clonality in P. humuli. In contrast, populations of P. cubensis display greater genotypic richness and heterozygosity compared with P. humuli and may be differentiated based on host, region of origin, and mating type (Kitner et al. 2015; Naegele et al. 2016; Ojiambo et al. 2015; Quesada-Ocampo et al. 2012). Two genetic clades are known to exist within P. cubensis (Runge et al. 2011), with distinct host ranges and mating types (Cohen et al. 2013, 2015; Quesada-Ocampo et al. 2012; Summers et al. 2015; Thomas et al. 2017a, b, c). Emergence of virulent pathotypes in Europe was coincident with genetic changes and potential hybridization of lineages (Kitner et al. 2015). In contrast, multiple investigations have found no evidence for pathotypes or races in P. humuli (Hoerner 1940; Royle and Kremheller 1981). In support of this, resistance to downy mildew in hop is quantitative (Henning et al. 2015) and has been remarkably stable in resistant cultivars over multiple decades and continents (Royle and Kremheller 1981).

The ultimate level for study of genetic polymorphism is directly ascertaining DNA for single nucleotide polymorphisms (SNPs) and other variants such as indels because these markers should be reproducible, codominant, and dispersed across the genome (Hartl and Clark 2007). Genotyping-by-sequencing (GBS) offers a cost-effective method for using high-throughput sequencing for genotyping strains (Elshire et al. 2011; Grünwald et al. 2016). This method reduces genome complexity by fragmenting DNA by a restriction enzyme or enzymes, followed by sequencing of each fragment and analysis of these sequences to identify variants. These variants can then be used to ask genetic questions about population patterns such as structure, linkage, or migration. GBS has been adopted for a range of plant pathogens, including Verticillium dahliae and Phytophthora rubi (Grünwald et al. 2017; Milgroom et al. 2014).

In oomycetes, sexual recombination is important for creating genetic diversity within pathogen populations and generation of survival structures that enable persistence during periods of unfavorable conditions (Grünwald et al. 2016). Given the importance of downy mildew on hop and limited understanding of population diversity of P. humuli, there is motivation to explore genetic diversity of the pathogen population to inform hop breeding efforts and also to infer the role of sexual recombination in P. humuli. To this end, the present study sought to quantify the degree of genetic diversity and population differentiation within and among P. humuli populations at multiple spatial scales. We tested the hypotheses that populations of P. humuli were predominately clonal and/or highly inbred but differentiated at the scale of individual hop yards. We expected that GBS would provide unprecedented resolution to find variation within clonal populations in space and time as mutations accumulate and populations diverge.

MATERIALS AND METHODS

Plant materials.

Plants of the downy mildew-susceptible cultivar Pacific Gem were propagated from softwood cuttings and maintained in a greenhouse free of downy mildew. Pacific Gem is highly susceptible to downy mildew and, to our knowledge, is universally susceptible to all isolates of the pathogen. The greenhouse was maintained at 20 to 25°C with a 14-h photoperiod. Plants were grown in 7.6-cm-diameter pots in Sunshine Mix #1 (Sun Gro Horticulture, Hubbard, OR) for 14 to 21 days before inoculation or use in an experiment. Watering was conducted daily and Sunshine Technigro 16-17-17 Plus fertilizer with micronutrients (Sun Gro Horticulture) was applied at each irrigation to promote leaves maximally susceptible to downy mildew.

Isolates and maintenance.

In 2010 and 2011, diseased shoots with symptoms and signs of systemic infection were collected from three hop yards in western Oregon, referred to herein as yards OR-7, Humpert, and Goulet; all were cultivar Nugget (Haunold et al. 1984). The Humpert and Goulet yards were grown by the same producer, whereas OR-7 was grown by a separate producer. All of the yards were planted in the early 1990s from rhizomes, which is typical because the plant is propagated vegetatively to ensure uniformity in a commercial field (Neve 1991). The original planting material for the Humpert and Goulet yards was derived from the same source; the planting material for OR-7 was sourced separately. The centroid of the OR-7 yard was located 8.7 km from the centroid of the Humpert yard and 15.2 km from the Goulet yard; the centroids of the Humpert and Goulet yards were 19.1 km apart.

In each yard, shoots with downy mildew were collected from each of 15 to 19 plants selected arbitrarily. An additional 9 to 13 diseased shoots were collected from each of three heavily diseased plants in each yard. Thus, the sampling schemes had a hierarchical structure, with one level of hierarchy for the 15 to 19 plants selected arbitrarily (i.e., yard scale) and two levels for the intensively sampled plants (i.e., yard scale and plants within yards).

Single-lesion isolates of P. humuli were derived using methods similar to those described by Summers et al. (2015). Diseased shoots were brought to a laboratory, the stems were placed in beakers of water, and each shoot was covered in a plastic bag overnight to induce sporulation. The following day, sporangia were washed from the abaxial leaf surface with sterile water applied from an airbrush sprayer, adjusted to a concentration of approximately 10,000 sporangia per milliliter, reinoculated onto one or several potted plants of Pacific Gem, and incubated in a sealed plastic container as described by Mitchell et al. (2011). The following morning, the plants were removed from the containers and allowed to air dry for 24 h. The plants were returned to the containers and placed into a growth chamber at 13 to 18°C for 6 to 8 days. Sporangia were harvested as described above, diluted to approximately 100 sporangia per milliliter, and inoculated onto three to four new plants. The following week, a plant bearing one lesion was selected and sporangia from the single lesion were harvested and again inoculated onto new plants. This process was repeated two times to reduce the likelihood of a single lesion containing two genotypes of P. humuli and to generate what was deemed a clonal isolate. Thereafter, each isolate was maintained separately on whole plants using the methods of Mitchell et al. (2011). Sporangia were stored at −80°C in Tris-EDTA buffer until DNA was extracted.

Given that P. humuli is an obligate parasite and there was evidence of clonality and/or a high degree of inbreeding as described below, extra steps were taken to include controls and ensure reproducibility of variant loci identified. A subset of the isolates obtained from the three yards noted above, as well as isolates from four other hop yards in Oregon and one yard in Washington, was obtained for the purpose of assessing genotyping variability within clonal progeny derived from a single lesion. For these isolates, diseased shoots were collected arbitrarily from within each yard, induced to sporulate overnight, and reduced to a single-lesion isolate as described previously. These isolates are referred to as the first “generation.” From each of the first-generation isolates, at least two single-lesion isolates were obtained and increased as second-generation isolates. In turn, a single-lesion isolate was obtained from each of the second-generation isolates, resulting in third-generation isolates. Thus, from each of the original isolates, there were two second-generation asexual progeny isolates and two third-generation progeny. To further assess technical errors (noise) associated with genotyping, 21 isolates were sequenced multiple times from the same library preparation. In total, there were 77 isolates included as clones or technical replicates. Similar to the other P. humuli isolates, sporangia were increased, harvested, and stored at −80°C until DNA was extracted.

Further, 15 isolates of P. humuli originating from 12 other hop yards (planted to various cultivars) in New York, Oregon, Vermont, Washington, and Wisconsin were included to estimate genetic differentiation among populations (Supplementary Table S1). We refer to these isolates as “other” populations. DNA also was obtained from two single-lesion isolates of P. cubensis originating from Ohio and North Carolina (isolates MSU-1 and ew045, respectively) representing clade 2 and clade 1, respectively (Runge et al. 2011; Thomas et al. 2017a). These isolates were genotyped and included in the analyses described below. Table 1 provides a summary of the isolates used in the analyses.

TABLE 1. Summary of the number of Pseudoperonospora humuli isolates used in each analysis^a

Yard/source name	Georeferencedb	Hierarchical sampling		Yard/source total	Technical replicates	Total
		Random	Intense
Goulet	Yes	13	30	43	5	48
Humpert	No	15	28	43	5	48
OR-7	Yes	15	31	46	23	69
Other	No			50	38	215

^a Note that the tally of isolates reflects the sum of the isolates retained after quality filtering. Isolates from Goulet, Humpert, and OR-7 hop yards were collected using a hierarchical sampling procedure in which a single isolate was collected from each of 13 to 15 hop plants at random (termed “Random”) and an additional 9 to 13 isolates were collected from each of three heavily diseased plants (termed “Intense”). Isolates labeled as “Other” were obtained from 12 other hop yards to estimate genetic differentiation among the three yards sampled hierarchically versus other populations.

^b The location of one of the intensely sampled hop plants in the Humpert yard was not recorded. Therefore, Mantel tests as described in the text were not possible for this yard.

View as image HTML

DNA extraction.

Sporangia of P. humuli were rinsed from affected leaves and collected by centrifugation to obtain approximately a 200-µl pellet. DNA was extracted from sporangia by the Oregon State University Center for Genome Research and Biocomputing (CGRB). Sporangia were lysed by grinding for 6 min in a Qiagen Tissue Lyser II (Qiagen Sciences, Germantown, MD). Genomic DNA was extracted using a Mag-Bind Plant DNA DS Kit (Omega Bio-Tek, Norcross, GA) and automated on a KingFisher Flex Purification System (Thermo Fisher Scientific, Waltham, MA).

Genotyping and variant calling.

DNA was submitted to the CGRB for GBS library construction (Elshire et al. 2011). In brief, DNA concentration was quantified using quantitative PCR (qPCR) with reagents from the Illumina Library Quantification Kit Universal qPCR Mix (Kapa Biosystems, Wilmington, MA). Samples were digested with ApeK1, and adapters were ligated. Library quality was assessed using the Agilent Bioanalyzer 2100 with the Agilent High Sensitivity DNA Kit (Agilent Technologies, Santa Clara, CA). Barcoded libraries with 48 or 49 isolates were sequenced on either the Illumina 2000 or Illumina 3000 sequencing system (Illumina, San Diego, CA) as 100- or 150-bp single end reads. Five sequencing runs were conducted in total. Sequence data were deposited in the NCBI Sequence Read Archive (accession PRJNA545725).

Reads were mapped to a draft reference genome of P. humuli (isolate 502AA; described in full in Rahman et al. 2019) using BWA aln version 0.7.5a-r405 (Li and Durbin 2009). This reference consisted of 40,506,691 bp arranged on 5,534 contigs. The reference had an N50 of 19,165 and a minimum contig size of 500 bp, and the longest contig was 195,040 bp. Variant calling was performed using TASSEL 3.0.168 (Bradbury et al. 2007) at the CGRB. Default TASSEL settings were used except as follows. The maximum number of good, barcoded reads per lane was increased to 500 million. The maximum number of tags the tags by tax (TBT) can hold while merging was increased to 50 million. The maximum number of sites output per contig (“chromosome” in the terminology used by TASSEL) was increased to 3 million. The minimum SNP call rate for a taxon to be included in the output was set to 0.1. The minimum minor allele frequency for reporting was set to 0.01. Linkage disequilibrium was not used to filter SNPs.

Quality control of variants and samples.

Variants obtained from the TASSEL variant calling pipeline were filtered for quality using vcfR (Knaus and Grünwald 2017) in R software (R Core Team 2018). Samples flagged by the TASSEL pipeline as failed were omitted from downstream analysis. A maximum depth of 100× was imposed to eliminate variants with unusually high coverage potentially attributable to mapping to repetitive regions. To retain variants that were from the base-ploid fraction of the reference, the variants were further filtered on high and low coverage. A 90% confidence interval was created for each library by taking the 0.05 and 0.95 quantiles. Variants above or below this threshold were flagged as missing data. Furthermore, variants were filtered based on a minimum read depth of 4. Isolates with >50% missing data were removed, as were variants with over 10% missing data across the remaining sample. After quality filtering, violin plots were constructed to visualize the depth of sequencing for variant sites to ensure sequencing was relatively unimodal.

Genotyping error is present in high-throughput sequence data as a result of several sources of technical error (Fox et al. 2014; Pfeiffer et al. 2018). Technical error in the data set was characterized in three ways: by repeated sequencing of technical replicates, by comparing variation in genotype calls among sequencing runs (batch effects), and by verification of allelic states in a subset of highly differentiated loci. To assess sequencing error remaining in the quality-filtered data set, histograms were constructed to express the number of pairwise differences in variants between clonal isolates and technical replicates sequenced at least twice. The distribution of variant differences between these samples was visually compared with the distribution of variant differences present in the entire data set. To characterize batch effects, Manhattan plots were constructed to detect patterns of variants with differentiation (based on Hedrick’s G’_ST; Hedrick 2005) using the genetic difference function implemented in vcfR (Knaus and Grünwald 2018). Verification of allelic state in a subset of highly differentiated loci was explored graphically. First, the distribution of G’_ST values among loci was quantified in a histogram. From this histogram, we evaluated in detail a manageable subset of loci that were highly differentiated. We examined the allelic state of loci that had G’_ST values of >0.3 near the extremes of the right tail of the histogram. This threshold was selected arbitrarily simply to subset the data, which resulted in 16 loci. The allelic states of these loci were visualized in matrix plots with reference to hop yard (for the three yards sampled hierarchically) and sequencing batch. As a further quality assurance step to ensure that the genotypes were reproducible, we created a similar matrix of allelic states for the same 16 loci for the technical replicates. This allowed comparison of genotyping error and provided a basis for disentangling the genetic differentiation between hop yards from background noise. Note here that these 16 loci were examined simply as representative variants to ensure reproducibility of the genotype calls, but all variants in the quality-filtered data set were used in subsequent analyses. Finally, to ensure that missing data were not systemically related to particular loci in some isolates, we created and inspected a matrix of allelic states for the first 1,000 variants for all isolates.

Population differentiation.

For the isolates sampled in a hierarchical manner, the degree of genetic differentiation among fields and plants sampled within fields was quantified using analysis of molecular variance (AMOVA) (Excoffier et al. 1992). From the variant call format (VCF) data, a matrix of pairwise distances was created using the function dist.gene() (Paradis et al. 2004) and AMOVA was performed using amova() (Paradis 2010). Two separate AMOVA calculations were performed, one with the isolates sampled arbitrarily across each of the three fields and the second based on the intensive sampling of the three heavily diseased plants in each of the yards. A total of 1,000 permutations were performed to establish significance. Ordination plots based on principal coordinates analysis and a matrix of pairwise distances were constructed to visualize the differentiation of isolates from within each of the three fields, as well as other isolates of P. humuli representing outgroups and the two isolates of P. cubensis. The distance matrix of pairwise distances was created using the function dist.gene() (Paradis et al. 2004), and principal coordinates analysis (classical multidimensional scaling) was performed using cmdscale() and plotted with ordiplot() (Oksanen et al. 2017). Ellipses of 90% confidence intervals were drawn on the ordination plots.

To infer reproductive mode, we calculated linkage disequilibrium to quantify whether variant alleles were randomly associated, as expected in randomly mating populations. The VCF data for the isolates from OR-7, Goulet, and Humpert yards were converted into a format readable by PLINK (Purcell et al. 2007) using VCFtools (Danecek et al. 2011). We then calculated R² for all pairwise combinations of variants using PLINK (Purcell et al. 2007). Output from PLINK was visualized in scatterplots and histograms constructed in R software.

Based on the results of the AMOVA (as described below), there was evidence of genetic differentiation among hop yards and plants within yards. Therefore, the relationship between genetic distance and physical distance of the plants from which the isolates were collected was calculated using Mantel tests (Mantel 1967). Distance matrices were again created with the function dist.gene() (Paradis et al. 2004) for genetic distances, Euclidean geographic distances were created with dist(), and the Mantel test performed using mantel.rtest() (Dray and Dufour 2007). The spatial location of one of the plants in the Humpert yard was not recorded; therefore, the analysis was conducted only for the OR-7 and Goulet yards, as noted in Table 1. Within a yard, significance was determined by conducting 999 permutations of the data as described in Dray and Dufour (2007). The relationship between Euclidean distance among plants and genetic distance was visualized in scatterplots.

RESULTS

Genotyping and variant calling.

Genotyping of 240 isolates resulted in 67,574 raw variants identified before quality filtering in a data matrix that contained 16.0% missing data. Violin plots of variants constructed for each isolate indicated a bimodal distribution typical for GBS data with other oomycete organisms (Knaus and Grünwald 2018), with major peaks near 10 to 30× sequencing depth (DP) and also above 100× depth for the variants sequenced, as reported in the TASSEL VCF file (Supplementary Fig. S1).

Quality control of variants and samples.

After filtering variants, 216 isolates were retained and the final data set consisted of 25,227 variants and a data matrix containing about 2.6% missing data. The distribution of variant read depth after filtering was unimodal, as determined by violin plots. The mean depth of sequencing ranged from 11.5 to 48.8× per sample. We found no evidence that missing data were systematically associated with certain populations (Supplementary Fig. S2). Certain variant loci were missing in most isolates, indicating poorly represented or sequenced variants, and certain isolates had a high degree of missingness, indicating poor sequencing of these isolates. Such variants and isolates were removed from the analysis based on the filtering criteria described above.

Histograms expressing the distribution of variants in the overall data set and subset of technical replicates were approximately normally distributed with a right skew (Fig. 1). The range and median of the distributions were broadly similar. Among the technical replicates, the median number of pairwise differences in variants identified was 41, whereas the median in the full data set was 52. Thus, the majority of genetic variation was attributable to sequencing or other sources of error.

Download as PowerPoint

Although most of the observed genetic variation was attributable to experimental error, we were able to identify variants that were associated with distinct populations. Importantly, these SNPs were highly reproducible in the technical replicates. To identify informative loci that differentiated populations from among the random genetic variation observed in the data set as a result of technical error, Manhattan plots were created to express G’_ST in relation to contig position. As expected, the plots of G’_ST values indicated that the majority of variants had G’_ST values near 0 and few variants were clustered at the same contig position (Fig. 2). However, there were some variants that were highly differentiated. As described previously, we examined 16 loci individually based on the threshold G’_ST we imposed to reduce data complexity in a matrix of allelic states for each isolate from the three hop yards in Oregon and 15 progeny derived from a single clonal isolate considered to be technical replicates (Fig. 3). The progeny derived from the same clonal isolate varied in only two of the 240 isolate × loci combinations (i.e., 15 isolates and 16 loci) and were independent of sequencing batch effects, indicating that the overall genotyping error rate was minimal. In contrast, the variants identified as highly differentiated based on G’_ST were indeed differentiated between hop yards and were not associated solely with sequencing batch. Thus, variant calls were deemed robust and reproducible.

Download as PowerPoint

Download as PowerPoint

Population differentiation.

Representative isolates of both clades of P. cubensis were well differentiated from all isolates of P. humuli based on principal coordinate analysis of pairwise genetic distances (Fig. 4). Among isolates of P. humuli collected from the three hop yards in Oregon, isolates from yard OR-7 were differentiated from isolates originating from the Goulet and Humpert yards. Among the latter two, the isolates from the Humpert yard were, with only one exception, contained within the 90% confidence ellipse of the Goulet yard. Isolates obtained from other yards in the Pacific Northwest or other areas of the United States spanned the 90% confidence ellipses of the three populations from the Oregon hop yards. The patterns in the ordination plots were verified by AMOVA. In the hierarchical sampling of the three hop yards in Oregon with the intensive sampling of heavily diseased plants, populations were differentiated among (P = 0.002) and within (P < 0.0001) hop yards (Table 2). However, only a total of 19.8% of the genetic variance was associated with differences among yards and 20.3% of the variance was among plants within the yard. Similarly, among the isolates collected individually and randomly from diseased plants, there was evidence for genetic differentiation at the yard scale (P < 0.0001) (Table 3).

Download as PowerPoint

TABLE 2. Analysis of molecular variance (AMOVA) results for the genetic variation of Pseudoperonospora humuli isolates collected using a hierarchical sampling approach from heavily diseased plants (referred to as intensively sampled yards) within each of three hop yards^a

AMOVA	Degrees of freedom	Sum of squares	Mean squares	σ2	P value	Component of variance (%)
Among hop yard	2	81,653	40,827	776.6	0.002	19.8
Within hop yard	6	105,443	17,574	899.4	<0.0001	20.3
Error	80	700,770	8,760	8,759.6	–	59.9
Total	88	887,866	10,090	–	–	–

^a Sample sizes used in this analysis included 30, 28, and 31 isolates collected from the hop yards indicated as Goulet, Humpert, and OR-7, respectively. Dashes indicate a statistic was not relevant.

View as image HTML

TABLE 3. Analysis of molecular variance (AMOVA) results for the genetic variation of Pseudoperonospora humuli isolates collected using a hierarchical sampling approach from randomly sampled plants within each of three hop yards^a

AMOVA	Degrees of freedom	Sum of squares	Mean squares	σ2	P value
Hop yard	2	147,912.4	73,956.2	3,894.1	<0.0001
Error	41	699,047.5	17,049.9	17,049.9	–
Total	43	846,959.9	19,696.7	–	–

^a Sample sizes used in this analysis included 13, 15, and 15 isolates collected from the hop yards indicated as Goulet, Humpert, and OR-7, respectively. Dashes indicate a statistic was not relevant.

View as image HTML

Analysis of linkage disequilibrium indicated that the overwhelming majority of variant loci in the populations had a nonrandom association, were strongly linked, and clustered at short genomic distances (Fig. 5). This finding is consistent with the expectation for an inbreeding and/or asexually reproducing population without random mating and assortment of alleles.

Download as PowerPoint

For the intensively sampled plants in the Goulet and OR-7 yards, there was evidence that genetic distance between isolates was significantly associated, albeit weakly, with physical distances of the plants from which the isolates were collected (Fig. 6). Monte Carlo simulations yielded simulated P values of 0.005 for the Goulet yard and 0.001 for the OR-7 yard. However, there was no evidence of an association between genetic distance and physical distance in any of the three yards based on the individual isolates collected randomly from plants (P ≥ 0.372 in all instances).

Download as PowerPoint

DISCUSSION

The sum of the data presented here support the hypothesis of a clonal or highly inbred population structure in P. humuli with small, yet significant, differentiation of populations among hop yards and among plants within yards. Evidence for clonality or inbreeding was demonstrated convincingly in that the number of variants identified was broadly similar between the full data set and the subset of only technical replicates (Fig. 1). Most variant loci had values of G’_ST near 0 and were scattered randomly across contigs, consistent with random error associated with sequencing or other nonbiological sources of variation. Importantly, there was genetic differentiation of populations at multiple spatial scales and we were able to confirm that genotyping was reproducible when isolates were sequenced repeatedly in multiple batches. Furthermore, there was a large magnitude of linkage disequilibrium among variant loci, consistent with repeated inbreeding and/or asexual recombination. Collectively, the overall limited genetic diversity within P. humuli but genetic differentiation in select variants found across the genome is consistent with repeated selfing and/or vegetative reproduction and with gradual differential accumulation of somatic mutations in different populations. This finding is in line with the homothallic mating system in P. humuli (Gent et al. 2017), as inbreeding can reduce genetic diversity and effective population sizes (Charlesworth 2003; Goodwin 1997).

Very little was known about the population structure of P. humuli before this study. One of the few reported studies on the population biology of P. humuli by Chee et al. (2006) reached quite different conclusions and suggested that the population in Oregon was genetically diverse and possibly reproducing sexually. The former research was conducted in the mid-1990s when genotyping methods based on random priming were some of the only tools available. Further, DNA was extracted from a heterogeneous source of inoculum derived from bulked inoculum collected without purification from diseased shoots. Given the lack of reproducibility of random priming methods (Penner et al. 1993) and the difference in how inoculum was obtained, we cannot directly compare the results of the present study to those of Chee et al. (2006).

The finding of limited genetic diversity within P. humuli populations is consistent with the findings of Summers et al. (2015) and Wallace and Quesada-Ocampo (2017). Although the sample sizes in these studies were relatively small, both found that P. humuli was well differentiated from the two clades of P. cubensis but with little evidence of differentiation of populations based on region of origin. Low population diversity of a pathogen in a specific geographic region often is considered circumstantial evidence for a nonnative origin of an organism (Goodwin 1997; Grünwald and Goss 2011; Prospero et al. 2007). The lack of diversity found in P. humuli within and across populations, both in the Pacific Northwest and apparently other production regions too, may indicate that the pathogen was introduced to North America. Hop downy mildew was first described in Asia in 1906 (Miyabe and Takahashi 1906) and was confirmed thereafter in eastern North America in 1909 and Western Europe in 1920 (Royle and Kremheller 1981). The disease was definitively reported in western Washington and Oregon in 1929 and 1930, respectively (Hoerner 1933). The limited genetic diversity of the pathogen observed in the United States and the timing of the disease reports in the early twentieth century point to an introduction of the pathogen. Presumably, P. humuli was imported on infected planting based on the geographic isolation of the hop industry in the U.S. Pacific Northwest. Given that hop is propagated vegetatively, subsequent spread in planting material within and across production regions seems likely and could be the subject of future research.

Although the P. humuli population appears to be clonally reproducing and/or highly inbred and possesses only limited genetic diversity among yards and even across states, GBS allowed us to discover fine-scale genetic differentiation at the yard and plant scales that suggests persistence of the pathogen over time. AMOVA indicated a significant effect of hop yard on the genetic variation of isolates whether collected randomly across diseased plants or intensively from heavily diseased plants. This finding can be explained by the overwintering biology of the pathogen. The primary means of perennation of P. humuli is believed to be as mycelia in the root system and buds of systemically infected plants (Coley-Smith 1962; Skotland 1961; Ware 1926), resulting in chronic infection of the same plant or yard (Coley-Smith 1964; Johnson and Anliker 1985). This systemic infection results in the occurrence of a few, heavily diseased plants in hop yards and pronounced aggregation of shoots with downy mildew among plants (Gent et al. 2011). Thus, founder genotypes of the pathogen introduced in association with infected planting material (Royle and Kremheller 1981) or dispersed from other sources (as is known in P. cubensis; Ojiambo et al. 2015, 2017) can be maintained because of chronic infection and persistence of the pathogen in a subset of heavily diseased plants if these plants are not rogued out (Coley-Smith 1962). Accumulation of somatic mutations over time attributable to systemic infections could result in differentiation of populations at the field and subfield scale, provided that migration of new genotypes between fields is limited (Hartl and Clark 2007).

In other oomycete pathogens, limited migration suppresses gene flow and increases genetic differentiation of populations (Goodwin 1997). Although P. humuli is capable of long-distance dispersal by aerial transport (Gent et al. 2009), most sporangia are expected to be deposited very near the source plant on which sporangia were produced and, paradoxically, also at long distances because of the dynamics of dispersal associated with an accelerating wave (Mundt et al. 2009). Within-field dispersal, colonization, and persistence of chronic infection of plants could produce a stepping-stone model type of gene flow that could result in founder events on individual plants with low within-population (i.e., plant) diversity (Charlesworth 2003). The Mantel tests with the intensively sampled plants support this idea, as the mean genetic distance of isolates within a plant was significantly associated with geographic distances of the plants from which the isolates were collected. We did not detect a significant association between genetic distance and geographic distance with the isolates collected randomly from diseased plants, perhaps reflecting that heterogeneous sources of inoculum can incite downy mildew within a yard but only systemic, chronic infections that result in the most severe disease provide the conditions or time necessary for populations to differentiate.

As noted before, the three hop yards sampled reflected two original sources of planting material: the same source for the Humpert yard and Goulet yard, and a different source for the OR-7 yard. It is striking that nearly all of the isolates from the Humpert yard were contained within the 90% confidence ellipse of the Goulet yard (Fig. 4), yet isolates from these two yards were differentiated from those originating from OR-7. The genetic similarity of the Humpert and Goulet yards is not explained by physical proximity, because these yards were further separated from each other than either was to the OR-7 yard. Given that infected planting material is an important source of primary inoculum of P. humuli, the data may suggest a founder population effect associated with infected rhizomes or potted plants. Based on only three hop yards, this suggestion clearly is speculative and further investigation is warranted. If confirmed, this could have profound implications for distribution of novel strains of the pathogen, such as those possessing fungicide resistance (Gent et al. 2008). Similarly, the lack of genetic variation observed in the 16 isolates obtained from other hop production areas of the United States may suggest dissemination of P. humuli across the country. More study is needed to understand the genetic diversity of the organism outside of the U.S. Pacific Northwest and dispersal routes of the pathogen.

A side note is worth making here on the importance of replication and experimental design when conducting population studies using GBS or other methods that utilize high-throughput sequencing where sequencing error is invariably present. This is especially critical with obligately biotrophic organisms with limited genetic diversity as found here. Detection of genetic variation at the field and subfield levels requires a high-resolution genotyping method, but also rigorous controls and quality assurance steps to ensure confidence in variants that are identified. Technical error is ever present in high-throughput sequencing data and efforts need to be made to quality assure data before analysis (Grünwald et al. 2017). Indeed, most genetic variants identified in the full data set could be attributed to technical error owing to contamination when working with obligate pathogens, sequencing, or other sources (Fig. 1), yet a small number of loci were responsible for the differentiation observed among populations and consistently segregated among populations (Fig. 2). Error rates associated with high-throughput sequencing technology based on sequencing-by-synthesis have been estimated near 0.1 to 0.2% per base (Fox et al. 2014; Pfeiffer et al. 2018). If we assume that the clonal progeny used as technical replicates in this research were genetically identical, the two misclassified variants as shown in Figure 3 suggest an error rate of roughly 0.1% in the quality-filtered data (i.e., 2 genotyping errors out of 16 loci × 15 isolates × 2 alleles per isolate × 4 possible nucleotides). Only because of inclusion of technical replicates throughout this study were we able to ascertain with confidence that the variants observed were consistent and not simply artifacts of batch effects. A challenge for most studies is that standards do not exist for level of missing data, depth of sequence coverage, and other quality assurance metrics. Quality assurance steps are even more imperative for obligately parasitic organisms because of the difficulty of ensuring a reliable reference genome free of nontarget organisms such as epiphytic bacteria associated with the host. We benefited from having a carefully produced and quality-assured reference genome (Rahman et al. 2019), but this is often not available and quality assurance of any reference is essential. Best practices for quality assurance are needed as population genomic studies become more routine and sequencing technologies continue to develop. Until best practices become more standard, approaches to assess error similar to those reported here can provide some confidence that genotype calls are reproducible and biologically based.

ACKNOWLEDGMENTS

Appreciation is extended to Lina Quesada-Ocampo for access to the P. humuli genome sequence data and the P. cubensis isolates used in this research. We also thank individuals that provided other isolates.