Genotyping by sequencing suggests overwintering of Peronospora destructor in southwestern Québec, Canada.

Several Peronospora species are carried by wind over short and long distances, from warmer climates where they survive on living plants to cooler climates. In eastern Canada, this annual flow of sporangia was thought to be the main source of Peronospora destructor responsible for onion downy mildew. However, the results of a recent study showed that the increasing frequency of onion downy mildew epidemics in eastern Canada is associated with warmer autumns, milder winters, and previous year disease severity, suggesting overwintering of the inoculum in an area where the pathogen is not known to be endogenous. In this study, genotyping by sequencing was used to investigate the population structure of P. destructor at the landscape scale. The study focused on a particular region of southwestern Québec-Les Jardins de Napierville-to determine if the populations were clonal and regionally differentiated. The data were characterized by a high level of linkage disequilibrium, characteristic of clonal organisms. Consequently, the null hypothesis of random mating was rejected when tested on predefined or nonpredefined populations, indicating that linkage disequilibrium was not a function of population structure and suggesting a mixed reproduction mode. Discriminant analysis of principal components performed with predefined population assignment allowed grouping P. destructor isolates by geographical regions, while analysis of molecular variance confirmed that this genetic differentiation was significant at the regional level. Without using a priori population assignment, isolates were clustered into four genetic clusters. These results represent a baseline estimate of the genetic diversity and population structure of P. destructor.

INTRODUCTION

Population genetics is essential to understand how ecological and evolutionary processes influence plant pathogens' temporal and spatial genetic diversity (Grünwald et al., 2017; Milgroom, 2017). This genetic diversity in plant pathogen populations results from migration, random genetic drift, mutation, recombination, or natural selection (Milgroom, 2017). It is generally characterized by the presence of polymorphisms that can be identified through genetic markers, the most commonly used being microsatellites and single‐nucleotide polymorphisms (SNPs). Some markers, such as microsatellites and simple‐sequence repeats, are well suited to study model organisms (i.e., a well‐studied organism for which many tools and resources exist) because they require prior genetic knowledge. Others, such as random amplified polymorphic DNA and amplified fragment length polymorphisms, allow the amplification of DNA fragments of unknown sequence (Milgroom, 2017; Russell et al., 2017). Nevertheless, in both cases, the approach does not allow for deep genomic coverage, potentially limiting the study of organisms with no or limited prior genetic knowledge. Recent developments in high‐throughput DNA sequencing tools such as genotyping by sequencing (GBS) (Elshire et al., 2011) allow for simultaneous exploration of genome‐wide genetic variation over hundreds of loci in a large number of individuals or populations. Hence, these tools enable a finer understanding of the relationships between genetic, temporal, and geographic variations for nonmodel organisms (Narum et al., 2013). Moreover, the availability of R codes (The Comprehensive R Archive Network) and resources (Jombart, 2008; Jombart et al., 2010; Kamvar et al., 2014; Knaus & Grünwald, 2017; Paradis et al., 2017) facilitates the use of these approaches for nonmodel organisms with a predominantly asexual mode of reproduction.

Many crop pathogens are members of the Peronosporaceae family (Crandall et al., 2017). In this family, particular attention has been given to members of the genus Phytophthora, especially Phytophthora infestans (Danies et al., 2015; Goss et al., 2014; Hansen et al., 2016; Knapova & Gisi, 2002; Montarry et al., 2010; Montes et al., 2016; Sjöholm et al., 2013; Widmark et al., 2007), Phytophthora sojae (Cai et al., 2019; Gally et al., 2007; Stewart et al., 2013, 2015; Wu et al., 2017), and Phytophthora ramorum (Gagnon et al., 2016; Goss et al., 2009; Ivors et al., 2006; Mascheretti et al., 2008; Prospero et al., 2004, 2007; Vercauteren et al., 2010). Microsatellites and GBS have been used to provide more accurate estimates of evolutionary processes by identifying patterns in genetic structures, describing spatial and temporal patterns in pathogen populations and clonal dynamics, inferring the role of sexual recombination, and identifying the origin of local populations. For example, GBS data revealed smaller genetic clusters within some of the dominant P. infestans clonal lineages. In addition to being sourced from neighbouring regions through long‐distance dispersal, these results suggest overwintering within regions of the United States where it was not expected (Hansen et al., 2016).

Despite their importance, with more than 700 species, obligate biotrophic downy mildews have received less attention in population genetics research than other groups of pathogens. The adaptations that led to their biotrophic lifestyle have compromised the efficiency of metabolic pathways necessary for saprophytic growth and thus altered their ability to grow on culture media (Thines & Choi, 2015). Hence, the inability of maintaining these organisms on artificial media limited their study. Typically, downy mildews produce numerous wind‐dispersed sporangia (asexual spores), spreading the inoculum over short or long distances during the growing season. It is generally accepted that downy mildew pathogen populations are mostly clonal (Bhattarai et al., 2020; Blanco‐Meneses et al., 2018; Gent et al., 2019; Lyon et al., 2016; Summers, Gulliford, et al., 2015). However, they are also capable of genetic recombination through sexual reproduction (homothallism and heterothallism), producing oospores, which are long‐term survival structures that can, like Peronospora destructor, survive up to 25 years in soils (McKay, 1957). Hence, given their abundant production of sporangia, their short infection–sporulation cycles, and their high host plant specificity, downy mildew pathogens are particularly destructive and of concern for several crops of economic importance, notably Plasmopara viticola in grapevine (Carisse et al., 2020), Bremia lactucae in lettuce (Dhar et al., 2019; Fall, Van der Heyden, Beaulieu, et al., 2015), Pseudoperonospora cubensis in cucurbits (Granke & Hausbeck, 2011; Granke et al., 2014; Summers, Adair, et al., 2015), Peronospora effusa in spinach (Klosterman et al., 2014), and P. destructor in onion (Fujiwara et al., 2021; Van der Heyden, Bilodeau, et al., 2020; Van der Heyden, Dutilleul, et al., 2020).

Variations between isolates of a given downy mildew species with respect to pathogenicity or aggressiveness on specific varieties of their host are commonly observed. For some species, genetic diversity has led to the description of races or pathotypes (Crandall et al., 2017). This is the case for P. effusa, for which at least 14 races related to aggressiveness on given spinach cultivars have been identified (Feng et al., 2014). P. cubensis is another example where population markers have been associated with mating types or clades, which are associated with host specificity (Cohen & Rubin, 2012; Wallace et al., 2020). In other cases, such as Pseudoperonospora humuli, population markers suggested genetic differentiation in populations at a very small scale: the hop yard (Gent et al., 2019). In addition to helping refine the taxonomy of these pathogens, the knowledge of population genetics is useful to improve biosurveillance (Rahman et al., 2020) and, ultimately, cultural practices and control of downy mildew diseases. However, as for many of the downy mildews, this type of genetic information is not yet available for P. destructor, the causal agent of onion downy mildew (ODM).

ODM has been reported worldwide in almost all the regions where onions are grown (Thines & Choi, 2015). It is particularly damaging in cool and humid climates, making it one of the main threats to onion production in northern Europe and North America. In Canada, onion is an important vegetable with a farm gate value of >C$95 million (Mailvaganam, 2017). The production is largely concentrated in the provinces of Québec and Ontario. In Québec, onion is produced in basins of muck soils (chernozem) located about 50 km southwest of Montreal, the most important ones being those of Sherrington (SHERR), Napierville (NAP), and Sainte‐Clotilde (CLO). In these areas, initial inoculum sources are unknown, and P. destructor is generally considered a periodically introduced onion pathogen. However, during the last decade and after a 10‐year ODM‐free period, ODM epidemics have been more frequent, started earlier, and were more severe than in the previous two decades (Van der Heyden, Dutilleul, et al., 2020). In a recent study, the polyetic component of ODM epidemiology was shown to be mostly influenced by the precipitation regime, the regional disease incidence towards the end of the previous year, warmer temperatures in the previous autumn, and warmer winters (Van der Heyden, Dutilleul, et al., 2020). These results support the hypothesis that P. destructor may be overwintering locally, although the authors could not, from the data available, determine whether sexual reproduction was involved or not. Potential sources of inoculum are systemically infected plants (onion sets, volunteer plants, or waste piles), airborne sporangia, and oospores present in soil or plant debris (Hildebrand & Sutton, 1980; McKay, 1957; Palti, 1989). Real‐time quantitative PCR (qPCR) analysis of soil samples suggested the presence of quantifiable DNA in almost 11% of the soil samples analysed, especially samples taken from fields with a short crop rotation (Fujiwara et al., 2021; Van der Heyden, Bilodeau, et al., 2020). In addition, these results suggested local overwintering but did not provide information about population structures.

Given the importance of ODM and the lack of knowledge on P. destructor population diversity, there was a need to explore the genetic variation of P. destructor populations and to provide evidence of local overwintering. In this study, GBS was used to investigate the population structure of P. destructor at the landscape scale. Thus, the objective of this research was to estimate the level of genetic diversity and differentiation between and within populations. The study focused on a particular region of southwestern Québec (Canada)—Les Jardins de Napierville—where isolates were collected between 2016 and 2019 and analysed together with Californian isolates as an outgroup to determine if the populations were clonal and regionally differentiated. Here, we hypothesize that although P. destructor populations are predominantly clonal, these will be differentiated at the scale of the production basins if P. destructor survives the winter locally. A fine resolution such as that obtained by GBS will allow the identification of mutations that accumulate within clonal populations as they diverge in space. The results presented in this study represent a baseline estimate of the genetic diversity and population structure of P. destructor and will improve our understanding of its epidemiology.

RESULTS

Genotyping and variant calling

The Ion Proton protocol used at the Institute of Integrative Biology and Systems (IBIS) provided an average of 32.5 million reads per sequencing chip, for a total of 52.3 million and 77.9 million reads for the first and second run, respectively (Figure S1). The difference in the number of reads between the two runs was largely attributable to the first chip, which produced only 15 million reads, as opposed to the other three, which produced just under 40 million reads each. After demultiplexing, four isolates with low or no reads were dropped, and 50.9 million and 76.1 million reads were retained for the first and second run. Ultimately, the average number of reads per isolate was 0.5 million (Figure S2a). After alignment to the reference genome (Natesan et al., 2020), an average of 0.24 million reads per isolate was kept, and the mean sequencing depth per isolate was 6.75× (Figure S2b,c).

Quality control and filtering

Prior to quality filtering, genotyping of the remaining 112 isolates yielded a total of 45,792 raw variants in a data set containing 21% missing data. Violin plots showing the distribution of sequencing depth among the variants for each isolate suggest a unimodal distribution, with data generally clustered around a median of 10× sequencing depth (Figure S3). After the filtering steps, an additional nine isolates that did not meet the quality criteria were dropped, resulting in a data set containing 103 isolates (Table 1), 5,335 variants, only 2.6% missing data, and a mean sequencing depth per isolate of 17.3× (standard deviation 8.99). The data showed important levels of heterozygosity indicated by a clear peak close to 50% allelic frequency (Figure S4). Linkage disequilibrium analysis suggested that a large proportion of the SNPs had a nonrandom association and were strongly linked (Figure 1a). This linkage did not appear to decay with physical distance (Figure 1b). After removing all SNPs in perfect linkage disequilibrium, 1340 variants were retained. At this step, there was no evidence that missing data were associated with a given population, that variant loci were missing across isolates, or that given isolates had too much missing data (Figure S5).

TABLE 1

Description of the isolates of Peronospora destructor kept after quality filtering

Sample	Year	Grower	Field	CROP	Location	Latitude	Longitude
Qc1	2016	VANW	32	Dry_bulb	SHERR	45.11681	−73.31416
Qc2	2016	BARA	CAMS1	Dry_bulb	CLO	45.15027	−73.40735
Qc3	2016	LECL	30	Bunching	SHERR	45.12344	−73.38038
Qc4	2016	BARR	B5A	Bunching	SHERR	45.09575	−73.35803
Qc5	2017	VANW	20	Dry_bulb	SHERR	45.11662	−73.30671
Qc7	2016	VANW	52	Dry_bulb	SHERR	45.13151	−73.34419
Qc8	2016	BARR	AB1	Bunching	CLO	45.17563	−73.40075
Qc9	2017	VANW	51	Dry_bulb	SHERR	45.13151	−73.34611
Qc11	2016	SOLE	20	Dry_bulb	SHERR	45.11090	−73.38294
Qc12	2017	VANW	31	Dry_bulb	SHERR	45.11420	−73.31577
Qc13	2017	GUEJ	G5	Dry_bulb	SHERR	45.11228	−73.32520
Qc15	2017	VANW	71	Dry_bulb	NAP	45.11505	−73.25197
Qc16	2016	BARR	TMB5	Bunching	CLO	45.15516	−73.40779
Qc17	2018	VANW	52	Dry_bulb	SHERR	45.13114	−73.34353
Qc19	2017	ZUMS	3	Bunching	SHERR	45.13662	−73.37155
Qc20	2017	ZUMS	D2	Bunching	SHERR	45.12199	−73.36379
Qc22	2016	VANW	51	Dry_bulb	SHERR	45.13044	−73.34417
Qc23	2017	GUEJ	B1	Dry_bulb	SHERR	45.12267	−73.37911
Qc24	2017	ZUMS	D7	Bunching	SHERR	45.12047	−73.36288
Qc25	2017	GUEJ	B4	Dry_bulb	SHERR	45.12085	−73.37868
Qc28	2017	GUEJ	G1	Dry_bulb	SHERR	45.11314	−73.32696
Qc29	2017	GUEJ	L5	Dry_bulb	SHERR	45.12881	−73.38487
Qc30	2016	VANW	60	Dry_bulb	NAP	45.11515	−73.25109
Qc31	2016	DELF	C48	Dry_bulb	NAP	45.16394	−73.12679
Qc32	2017	DELF	B36	Dry_bulb	NAP	45.17559	−73.12484
Qc33	2019	DELF	B30	Shallot	NAP	45.17559	−73.12484
Qc34	2019	DELF	B30	Shallot	NAP	45.17583	−73.12418
Qc35	2019	GUEJ	G1	Dry_bulb	SHERR	45.11207	−73.32627
Qc36	2019	GUEJ	G1	Dry_bulb	SHERR	45.11257	−73.32737
Qc37	2019	JSTC	3	Dry_bulb	CLO	45.15072	−73.41980
Qc38	2019	JSTC	3	Dry_bulb	CLO	45.15022	−73.42092
Qc39	2019	LAFI	5	Dry_bulb	CLO	45.14529	−73.44826
Qc40	2019	LAFI	5	Dry_bulb	CLO	45.14492	−73.44770
Qc42	2019	VANW	8	Dry_bulb	SHERR	45.10591	−73.30955
Qc43	2019	VERN	221	Dry_bulb	NAP	45.13999	−73.27191
Qc44	2019	VERN	221	Dry_bulb	NAP	45.14076	−73.26942
Qc49	2019	JSTC	J4	Dry_bulb	CLO	45.13986	−73.43418
Qc50	2019	JSTC	J4	Dry_bulb	CLO	45.13987	−73.43606
Qc52	2019	JSTC	J4	Dry_bulb	CLO	45.13913	−73.43502
Qc53	2019	ZUMS	ND2	Bunching	SHERR	45.12224	−73.36353
Qc45	2016	JSTC	J4	Dry_bulb	CLO	45.13913	−73.43502
Qc46	2019	REMR	R31	Bunching	CLO	45.15657	−73.38100
Qc47	2019	ZUMS	ND2	Bunching	SHERR	45.12224	−73.36353
Qc48	2019	REMR	R31	Bunching	CLO	45.15680	−73.38147
Qc51	2019	JSTC	J4	Dry_bulb	CLO	45.14079	−73.43529
Qc55	2019	BARR	TMB2	Bunching	CLO	45.15443	−73.40884
Qc57	2019	BARR	ORR2	Bunching	SHERR	45.08856	−73.34918
Qc60	2019	JSTC	LABBE1	Dry_bulb	CLO	45.13882	−73.43336
Qc61	2019	JSTC	LABBE1	Dry_bulb	CLO	45.14196	−73.42873
Qc62	2019	VANW	5	Dry_bulb	SHERR	45.11398	−73.30016
Qc63	2019	VANH	J3	Dry_bulb	NAP	45.12609	−73.23789
Qc64	2019	VANH	J3	Dry_bulb	NAP	45.12582	−73.23847
Qc65	2019	VANW	8	Dry_bulb	SHERR	45.10601	−73.30834
Qc66	2019	VANW	8	Dry_bulb	SHERR	45.10547	−73.30847
Qc67	2019	VANW	8	Dry_bulb	SHERR	45.10531	−73.30927
Qc68	2019	VANW	8	Dry_bulb	SHERR	45.10635	−73.30901
Qc69	2019	REMR	M1	Bunching	CLO	45.17042	−73.39599
Qc70	2019	REMR	M1	Bunching	CLO	45.17078	−73.39497
Qc71	2019	REMR	M1	Bunching	CLO	45.16991	−73.39576
Qc72	2019	JSTC	3	Dry_bulb	CLO	45.14960	−73.42096
Qc73	2019	JSTC	3	Dry_bulb	CLO	45.14941	−73.42337
Qc74	2019	JSTC	3	Dry_bulb	CLO	45.14895	−73.42228
Qc75	2019	GUEA	PB4	Dry_bulb	CLO	45.14633	−73.37494
Qc76	2019	BARR	GF28	Bunching	CLO	45.13927	−73.40280
Qc77	2019	LECL	9_11	Bunching	CLO	45.14633	−73.37494
Qc78	2019	LECL	9_11	Bunching	CLO	45.14633	−73.37494
Qc79	2019	LAFI	6_7	Dry_bulb	CLO	45.14591	−73.45108
Qc80	2019	LAFI	6_7	Dry_bulb	CLO	45.14562	−73.45117
Qc81	2019	LAFI	6_7	Dry_bulb	CLO	45.14495	−73.45029
Qc82	2019	LAFI	6_7	Dry_bulb	CLO	45.14501	−73.44979
Qc83	2019	GUEJ	L1	Dry_bulb	SHERR	45.13449	−73.37605
Qc84	2019	GUEJ	B1	Dry_bulb	SHERR	45.12173	−73.37968
Qc85	2019	GUEJ	B4	Dry_bulb	SHERR	45.12044	−73.37927
Qc86	2019	GUEJ	B4	Dry_bulb	SHERR	45.12143	−73.37844
Qc87	2019	GUEJ	B4	Dry_bulb	SHERR	45.12138	−73.37919
Qc88	2019	DELF	B30	Shallot	NAP	45.17623	−73.12471
Qc89	2019	DELF	B30	Shallot	NAP	45.17608	−73.12562
Qc90	2019	DELF	B30	Shallot	NAP	45.17515	−73.12522
Qc91	2019	DELF	B30	Shallot	NAP	45.17437	−73.12486
Qc92	2019	DELF	B30	Shallot	NAP	45.17432	−73.12570
Qc93	2019	DELF	B30	Shallot	NAP	45.17353	−73.12524
Qc94	2019	DELF	B30	Shallot	NAP	45.17352	−73.12627
Qc95	2019	DELF	B30	Shallot	NAP	45.17284	−73.12630
Qc96	2019	DELF	B30	Shallot	NAP	45.17406	−73.12657
Qc97	2019	DELF	B30	Shallot	NAP	45.17437	−73.12486
Qc98	2019	LAFI	5	Dry_bulb	CLO	45.14512	−73.44828
Qc99	2019	LAFI	5	Dry_bulb	CLO	45.14567	−73.44899
Qc100	2019	LAFI	5	Dry_bulb	CLO	45.14612	−73.44887
Qc101	2019	LAFI	5	Dry_bulb	CLO	45.14529	−73.44873
Qc102	2019	LAFI	5	Dry_bulb	CLO	45.14558	−73.44837
Qc103	2019	LAFI	5	Dry_bulb	CLO	45.14580	−73.44934
Qc104	2019	LAFI	5	Dry_bulb	CLO	45.14589	−73.44903
Qc105	2019	GUEJ	G1	Dry_bulb	SHERR	45.11244	−73.32758
Qc106	2019	GUEJ	G1	Dry_bulb	SHERR	45.11188	−73.32755
Qc107	2019	GUEJ	G1	Dry_bulb	SHERR	45.11192	−73.32688
Qc108	2019	GUEJ	G1	Dry_bulb	SHERR	45.11132	−73.32650
Qc109	2019	GUEJ	G1	Dry_bulb	SHERR	45.11154	−73.32603
Qc110	2019	GUEJ	G1	Dry_bulb	SHERR	45.11099	−73.32563
Qc111	2019	GUEJ	G1	Dry_bulb	SHERR	45.11111	−73.32485
Qc112	2019	GUEJ	G1	Dry_bulb	SHERR	45.11161	−73.32531
CA1	2019	CA	2019	CA	CA	–	–
CA2	2019	CA	2019	CA	CA	–	–
CA3	2019	CA	2019	CA	CA	–	–

SHERR, Sherrington; CLO, Sainte‐Clotilde; NAP; Napierville; CA, California.

FIGURE 1

Histogram showing the amplitude of linkage disequilibrium (r 2) between the different loci in Peronospora destructor isolates (a) and the distribution of linkage disequilibrium as a function of genetic distance between associated loci (b)

Manhattan plots were used to represent the distribution of the G′ST statistic as a function of genomic position for the complete data set and the targeted subset. These plots showed that most variants had G′ST values near zero for both the complete data set and the targeted subset (Figure 2). However, among all the variants obtained, some were also found to be highly differentiated. In the end, 158 and 168 informative loci were kept for use in subsequent analyses for the complete and the targeted data set.

FIGURE 2

Manhattan plots showing the distribution of the G′ST statistics as a function of genomic position constructed with (a) the complete data set and (b) the targeted subset. These results show that most of the variants had G′ST values near 0 (no or low differentiation), while a reasonable subset of well‐differentiated variants was also present in the data

Genetic differentiation of P. destructor populations

The analysis of molecular variance (AMOVA) performed with the regional data set revealed that the populations were significantly differentiated at the regional level (p = 0.002), with 17.93% of the genetic variance observed between regions (Table 2). Similarly, for the targeted data set, a significant (p = 0.005) share of the genetic variance (23.98%) was observed between fields (Table 3). The result of the AMOVA revealed a significant but limited effect of year on population structure (p = 0.041), with 8.47% of the variation observed among years (Table S1). The AMOVA performed with the crop as classification factor to define populations did not show significant differentiation among crops (p = 0.581) (Table S2).

TABLE 2

Results of the analysis of molecular variance (AMOVA) conducted with Peronospora destructor isolates from the arbitrary sampling carried out in Québec (regional data set)

AMOVA	df	Sum of squares	Mean squares	σ^2	Variance component (%)	p value
Between regions	2	150.374	75.187	3.066	17.93	0.002
Within regions	66	926.149	14.032	14.032	82.07
Total	68	1076.524	15.831	17.099

TABLE 3

Results of the analysis of molecular variance (AMOVA) conducted with Peronospora destructor isolates collected from the three heavily infested fields in southwestern Québec (targeted data set)

AMOVA	df	Sum of squares	Mean squares	σ^2	Variance component (%)	p value
Between fields	2	187.672	93.836	6.988	23.98	0.005
Within fields	28	620.410	22.158	22.158	76.02
Total	30	808.082	26.936	29.145

The pairwise AMOVAs for the regional data set suggested a significant level of genetic differentiation among local populations. The largest genetic differentiation was found between the NAP and SHERR populations = 0.284, p = 0.006), whereas the genetic differentiation between the SHERR and CLO populations, though smaller, was also significant = 0.153, p = 0.004). In contrast, genetic differentiation between the CLO and NAP populations was not significant = 0.006, p = 0.292) (Table 4). Similar pairwise AMOVA results were obtained for the targeted and Québec data sets, with the same interpretation in terms of significant or nonsignificant genetic differentiation between local populations (Table 4). When the regional data set was split into 2019 and 2016–2018, the pairwise AMOVA results showed that the NAP and SHERR populations were significantly differentiated regardless of the grouping, while the NAP and CLO populations were not differentiated. The SHERR and CLO populations were significantly differentiated for the targeted group in 2019 but not in 2016–2018 (Table 4).

TABLE 4

Pairwise index of differentiation for the regional, targeted, and Québec data sets

Population pair	Regional		2016–2018		2019		Targeted		Québec b
	Φ_PT a	p value	Φ_PT	p value	Φ_PT	p value	Φ_PT	p value	Φ_PT	p value
NAP–SHERR	0.284	0.006	0.486	0.008	0.408	0.009	0.176	0.039	0.150	0.011
NAP–CLO	0.006	0.292	0.127	0.171	0.084	0.196	0.110	0.120	0.026	0.275
SHERR–CLO	0.153	0.004	0.025	0.440	0.499	0.001	0.345	0.006	0.314	0.0001

is the pairwise index of population differentiation.

Includes all samples from the targeted and regional data sets.

Genetic structure of P. destructor populations

In the discriminant analysis of principal components (DAPC) with a priori population assignment, 19 principal components (PCs) were retained for the complete data set (highest mean of success [HMS] = 0.660, root mean square error [RMSE] = 0.345). This result indicated that isolates obtained from California are well differentiated from all the isolates collected from the three muck soil basins in southwestern Québec (Figure 3a). For the Québec data set, 25 PCs were retained in the DAPC (HMS = 0.635, RMSE = 0.402) (Figure 3c). For this data set, 76.9% of the isolates were correctly assigned to the populations defined a priori. More specifically, 85.4%, 67.6%, and 72.2% of correct classification were obtained for the SHERR, CLO, and NAP populations, respectively. The posterior membership probabilities of <90% suggested admixed genotypes in the populations (Figure 3d). The results of the unweighted paired group method with arithmetic mean (UPGMA) dendrogram based on Provesti's genetic distance are consistent with the results obtained with both AMOVA and DAPC; it reveals two main clusters, CA and Québec, with the latter being further subdivided into two other clusters, CLO‐SHERR and NAP (Figure 3b).

FIGURE 3

Discriminant analysis of principal components (DAPC) conducted with a priori population assignment. In (a), the analysis was conducted on the complete data sets with 19 principal components and two discriminant functions, chosen following a cross‐validation procedure. In (b), the dendrogram of the relationships among Peronospora destructor populations is based on Provesti's genetic distance. In (c), the analysis was conducted on the Québec data set, with 25 principal components and two discriminant functions, chosen following a cross‐validation procedure. In (d), a bar graph shows the probability that each isolate belongs to a designated population

The genetic structure of P. destructor populations was also tested using nonparametric DAPC with no a priori population assignment. The comparison of Bayesian information criterion (BIC) values for 1 to 10 clusters suggests a rapid decrease in the BIC value from k = 1 to k = 4, suggesting four genetic clusters containing 24, 13, 42, and 24 isolates (Figure 4a). However, the plot of the discriminant function shows three clear clusters with some other clusters embedded within them (Figure 4b). The first five eigenvalues accounted for 93.4% of the total variation. The UPGMA dendrogram built using Provesti's distance reveals two main clades of P. destructor isolates that did not necessarily correlate with population assignment (Figure S6).

FIGURE 4

Discriminant analysis of principal components (DAPC) conducted with no a priori population assignment. (a) Boxplot of the Bayesian information criterion (BIC) obtained for different numbers of clusters from 999 permutations, indicating the presence of four genetic clusters. (b) Scatter plot of the discriminant functions representing the main four genetic clusters of Peronospora destructor identified by DAPC

The heterozygosity value was significantly (p < 0.05) greater for the NAP population (H obs = 0.358) than for the SHERR (H obs = 0.315) and CLO (H obs = 0.314) populations. The standardized association index leads to an unambiguous rejection of the null hypothesis of sexual recombination in the four populations, including California (p < 0.001) (Table 5). When tested using, as populations, the four genetic clusters detected by DAPC with no a priori population assignment, the null hypothesis of random mating was also rejected. The calculated values ranged from 0.017 to 0.160 (p < 0.01), except for cluster 2, for which  = −0.027 (not significant). These results provide further evidence to reject the null hypothesis of random mating.

TABLE 5

Summary statistics for genetic diversity and linkage measurements for Peronospora destructor populations

Population	Heterozygosity a	Median pairwise genetic distance	r¯_d (p value)
SHERR	0.315 (0.280)	0.129	0.314 (0.01)
CLO	0.314 (0.292)	0.146	0.252 (0.01)
NAP	0.358 (0.343)	0.183	0.306 (0.01)
CA	0.359 (0.259)	0.129	0.739 (0.01)

Values in parentheses are expected heterozygosity.

Provesti's pairwise genetic distance.

For the Québec data set (r = 0.177) as well as the targeted data set (r = 0.239), the null hypothesis of absence of association between genetic and geographic distances was rejected by the Mantel test (Figure S7). The positive sign of the relationship (i.e., the higher one type of distance, the higher the other) indicated genetic isolation by distance. As expected, the Mantel correlogram showed a decrease of the Mantel statistic with increasing geographic distance. At spatial distances of <5.25 km, the Mantel statistic was significantly (p < 0.05) positive (Figure 5). At larger spatial distance classes, isolates had higher genetic distances than expected at random (p ≥ 0.05) (Figure 5).

FIGURE 5

Mantel correlogram plotting the correlation coefficient (r) between genetic distance and geographic distance in seven distance classes. The red and green squares indicate a significant correlation while the blue squares indicate that the correlation is not significant

DISCUSSION

To our knowledge, this is the first population genetic study of the obligate biotrophic plant pathogen P. destructor. Considering that the sporadic nature of the disease seems to evolve towards a more endemic style and that the disease is difficult to manage, knowing the sources of initial inoculum is crucial. This study aimed at providing evidence for local overwintering of P. destructor populations. Isolates obtained from single ODM lesions collected from the main Québec onion production areas were analysed by GBS, following whole‐genome amplification. Whereas dominant asexual reproduction leads to limited recombination and clonal population structures, local overwintering introduces genetic differentiation and a gradual accumulation of mutations within populations belonging to different environments. Our results support the hypothesis of a clonal population structure with limited but significant differentiation among production zones.

Clonal populations showed high linkage disequilibrium, meaning that the nonrandom association of alleles from different loci is higher than expected (Slatkin, 2008). As expected, the proportion of alleles in perfect linkage was high and occurred over large regions of the P. destructor genome, supporting the hypothesis of clonality (Milgroom, 2017). High linkage disequilibrium was also reported for P. humuli, known to have both clonal and homothallic reproduction modes (Gent et al., 2017, 2019). Such high linkage disequilibrium may introduce redundancy, especially for large genomic data from populations with small effective sample size, and this may limit the resolution in identifying informative loci (Bhattarai et al., 2020; Calus & Vandenplas, 2018). In our study, most of the variant loci had G′ST values close to 0 and were randomly distributed throughout the genome even after heavy data pruning, which is characteristic of random sequencing errors and abiotic variation sources (Gent et al., 2019). However, after thorough variant filtering, we still found significant genetic differentiation among P. destructor populations in southwestern Québec.

To date, there is very little information about population genetics for P. destructor. Nonetheless, our results are consistent with previous studies in Peronosporales reporting limited population genetic diversity. In particular, it has been shown that P. humuli populations are weakly differentiated and that eastern American populations tend to be genetically clustered. In contrast, western American populations are less differentiated (Summers, Gulliford, et al., 2015). Wallace and Quesada‐Ocampo (2017) suggested that P. humuli populations are less differentiated than P. cubensis populations, and their results were confirmed by Gent et al. (2019). Low genetic differentiation was reported for Australian populations of P. viticola, while populations from North America, identified as the centre of origin, were shown to have higher genetic diversity (Taylor et al., 2019).

Low genetic diversity is generally considered to be an indicator of recently introduced or re‐emerging plant pathogens in a given region (Grünwald & Goss, 2011; Grünwald et al., 2016). In his monograph, Yarwood (1943) reported the important impact of ODM on Californian onion seed production, mentioning that epidemics were irregular in frequency and variable in severity. P. destructor was sporadically reported in North America (i.e., New York State, Province of Ontario) in the 1970s and 1980s and gave rise to the most important outbreaks since the 1950s in 1977–1979 (Hildebrand & Sutton, 1982; Smith et al., 1985). Since then, P. destructor has been encountered only sporadically in these regions. In 1999, P. destructor was reported in sweet onion for the first time in Georgia (Langston & Sumner, 2000). This report was followed by two severe epidemics in that region in 2007 and 2012 (Parkunan et al., 2013). In southwestern Québec, the scenario was different: Important epidemics were reported between 1989 and 1993, followed by 11 consecutive years with almost no observed symptoms (Van der Heyden, Dutilleul, et al., 2020). Since the mid‐2000s, however, there has been a re‐emergence of P. destructor, with outbreaks reported every year and occurring earlier each year, suggesting that the pathogen has become endemic, unlike other production areas in North America (Van der Heyden, Dutilleul, et al., 2020). Thus, the limited genetic diversity observed in our study supports a recent re‐emergence of P. destructor in southwestern Québec.

Despite a predominant clonal mode of reproduction, the analysis of GBS data revealed significant genetic differentiation between the geographical regions studied, including southwestern Québec and California. The results obtained with the DAPC applied using a priori defined populations showed that the California isolates are clustered together and are more genetically distant than the isolates collected from southwestern Québec. These findings were expected because there is no or little plant material or seeds exchanged between these two regions. Our results also showed that the NAP and SHERR isolates clustered together, while the CLO isolates were less differentiated. These results are supported by the UPGMA dendrogram built from Nei's genetic distances, which first classified populations in two groups (Québec and California) and then the first group in two subgroups (CLO/SHERR and NAP).

Results of the AMOVA indicated a significant effect of the region and field on genetic variation when sampling was regional (Table 2) or targeted (Table 3), which implies that the pathogen persists over time (Gent et al., 2019). Moreover, the pairwise AMOVAs showed genetic differentiation between the same populations, whether with the regional or targeted data set. Similar results were found for P. humuli, which exhibited low but significant genetic differentiation at a finer scale, with almost 20% of the genetic variance associated with hop yards all located in western Oregon (Gent et al., 2019). Concerning other factors that might contribute to shaping of the population structure, we observed no effect of Allium species (A. cepa, A. fistulosum, or A. aggregatum), unlike other Peronosporales such as P. cubensis that showed host‐related genetic differentiation (Wallace et al., 2020). The AMOVA results suggested a limited but significant effect of year on genetic differentiation, which might indicate occasional influx of new genotypes from adjacent producing regions or the presence of sexual reproduction within the populations.

The re‐emergence of ODM in southwestern Québec has been associated with warmer weather conditions in autumn and winter, suggesting the overwintering of P. destructor inoculum (Van der Heyden, Dutilleul, et al., 2020). However, the form in which P. destructor survives from one season to the next remains poorly understood. The carryover of P. destructor as mycelium in dormant bulbs has been proposed as one mode of overseasoning (Hildebrand & Sutton, 1980; Yarwood, 1943). Systemically infected onions used for planting or left in the field after harvest could be the source of seasonal epidemics. Warmer weather conditions could also favour volunteer plants' survival and consequently help maintain the inoculum in the region, which is supported by the limited genetic diversity and clonal population structure that we observed.

A solely asexual mode of reproduction can hardly explain regional ODM epidemics in southwestern Québec (Van der Heyden, Dutilleul, et al., 2020). In our study, although genetic diversity between muck soil basins was significant, a greater genetic diversity was observed within a specific basin. Such genetic structure in P. destructor populations provides evidence supporting the role of local sources of inoculum to initiate ODM epidemics. High levels of linkage disequilibrium are generally considered an indicator of clonality. However, in addition to a lack of recombination, linkage disequilibrium can also be caused by selection, population admixture, or random genetic drift (Milgroom, 2017). In our study, because the standardized index of association was significantly different from zero, we rejected the null hypothesis of random mating for the four populations. When tested against the genetic clusters identified with no a priori defined population, the null hypothesis of random mating was also rejected, indicating that linkage disequilibrium was not a function of population structure but arose from a mixed mode of reproduction (Koenick et al., 2018). In addition, the observed differences between DAPC results obtained with and without a previously defined population also point towards population admixture.

The role of random mating and oospore production in P. destructor is insufficiently documented, especially in North America. It is generally accepted that the production of oospores is possible, and it is known that they do not occur regularly, can survive for a surprisingly long period (up to 25 years), and need a very long maturation period (up to 4 years) before they can germinate (McKay, 1937, 1957; Yarwood, 1943). McKay (1957) stated that oospores were not produced every year, but they could be found in enormous numbers when they were. Thus, rare events of sexual recombination leading to the production of oospores could constitute an important inoculum reservoir that could persist locally for several years. Hence, it is likely that P. destructor has a predominantly asexual mode of reproduction, with occasional sexual reproduction events.

It is not known whether P. destructor is heterothallic or homothallic, but it is known that the amount of heterozygosity is low in homothallic organisms, especially when populations are established for several generations (Goodwin, 1997), and in populations with a strictly asexual mode of reproduction (Goodwin, 1997; Van Poucke et al., 2021). However, the results obtained in this study showed that P. destructor displayed significant levels of heterozygosity, indicating that it would be heterothallic. These results are consistent with the levels of heterozygosity observed in other heterothallic oomycete species. For example, the observed levels of heterozygosity were 0.37 on average for P. effusa (Lyon et al., 2016) and 0.21 for Peronospora belbahrii (Thines et al., 2020), whereas it varied between 0.01 and 0.14 for different heterothallic Phytophthora species and between 0.00 and 0.06 for homothallic species (Goodwin, 1997).

Airborne Phytophthora spp. and Peronospora spp. sporangia are known to travel very long distances, but it is accepted that the majority of the produced inoculum is deposited locally over short distances (Aylor, 1986; Aylor et al., 2001; Aylor & Taylor, 1983; Corredor‐Moreno & Saunders, 2019; Fall, Van der Heyden, Brodeur, et al., 2015; Gent et al., 2009; Granke & Hausbeck, 2011). Our Mantel correlogram results support the hypothesis that most of the inoculum is dispersed over short spatial distances. The exponential‐like decrease showing higher Mantel statistic values in the first spatial distance classes suggests that there are spatial patches of genetic similarities, perhaps smaller than the grain of our study (i.e., within the muck soil basins) (Diniz‐Filho et al., 2013). Limited migration can significantly affect gene flow, leading to a local decrease in genetic diversity. The limited genetic diversity observed in P. destructor combined with significant diversity within populations and the fact that these patches of genetic similarity were observed could result from a recent founder event (Goodwin, 1997). These results are consistent with the presence of a symptom‐free period of nearly 10 years observed in southwestern Québec between 1993 and 2004, followed by a gradual increase in the magnitude of ODM epidemics (Van der Heyden, Dutilleul, et al., 2020).

In conclusion, the results obtained in this work can serve as a baseline estimate to describe the genetic diversity and population structure of P. destructor and thus contribute to our understanding of Peronosporales ecology. Because of limited numbers of isolates, we could not pinpoint the origin of the initial population but could provide plausible population structures and reproduction modes. Keeping in mind that population genetic analysis provides only indirect evidence of sexual or asexual reproduction modes, we conclude that the genetic structure of P. destructor populations in southwestern Québec is consistent with the accumulation of mutations leading to slow population divergence in time and space, while occasional sexual reproduction events also occur. Other factors, such as fungicide resistance, may contribute to shaping of the population structure, and these should be included in further studies (Vogel et al., 2020). Using a finer sampling scale in future research would help confirm the genetic similarity of neighbouring spatial patches and searching and finding oospores would be valuable to understand the role of overwintering inoculum. In addition, a comparative genomics study of isolates belonging to the four populations that we studied would be necessary to deepen our knowledge of the organism. Finally, our results should be considered for the management of P. destructor in onion production, especially for defining priority areas for intervention when the first symptoms are reported, crop rotation planning, and crop residue management.

EXPERIMENTAL PROCEDURES

Population sampling

The Jardins de Napierville county has about 10,000 ha of cultivated muck soil, of which around 2000 ha are dedicated to the cultivation of Allium each year (mainly dry bulb onions, French shallots, and green onions). The muck soils are distributed in several basins, the three most important being SHERR, NAP, and CLO. Hence, the populations are defined based on the geographic location of the muck soil basin in which they were sampled.

P. destructor isolates were collected from diseased onion leaves bearing fresh sporangia from 2016 to 2019. Sampling was performed randomly across the three muck soil basins in 41 fields belonging to 15 farms when agronomists and crop specialists reported symptoms. In 2019, a targeted sampling was conducted in three heavily infected fields, from which 10 to 12 isolates were collected on the same sampling date. In addition to the isolates collected in Québec as described, three isolates were obtained from Hartford County, California. For all isolates, sporangia were taken from single lesions by carefully touching them with a sterile culture swab (Fisher Scientific) and stored in 95% ethanol at −20℃.

DNA preparation and genotyping by sequencing

As the number of sporangia was limited in many isolates, a multiple displacement amplification‐based whole‐genome amplification (WGA) method was used to increase the genomic DNA for all isolates (Casso et al., 2019; Han et al., 2012). To do this, the Repli‐g single cell kit was used according to the manufacturer's instructions (Qiagen). After the WGA procedure, the presence of P. destructor DNA was tested using qPCR as described in Van der Heyden, Bilodeau, et al. (2020) and subsequently it was quantified using a Qubit fluorometer with the dsDNA Qubit HS test kit (ThermoFisher). The DNA concentration was adjusted to 10 ng/µl before being sent to the Genomic Analysis Platform of the IBIS for the GBS sequencing workflow (Elshire et al., 2011). Briefly, DNA was digested with ApeK1 for genome complexity reduction prior to adapter ligation. All libraries were sequenced on Ion Torrent P1V3 chips, targeting 60 to 80 million single‐end reads per chip. Each chip contained 58 isolates, and two rounds of sequencing (two separate chips) were conducted for each library for a total of four chips.

Data processing, variant calling, and quality control

Demultiplexing was performed using process_radtag from Stacks v. 2.53 (Catchen et al., 2011, 2013) and barcodes and sequencing adapters were trimmed using BBDuk v. 38.86 from the BBTools suite (Bushnell, 2014). Reads were mapped to a draft genome of P. destructor (GenBank accession no. WBRY00000000; Natesan et al., 2020) using Bowtie2 v. 2.4.1 with the “very‐sensitive‐local” option (Langmead & Salzberg, 2012). Identification of SNPs in the metapopulation for each locus was performed using gstack (Stacks v. 2.53). The program provided a first catalogue containing the consensus sequence for each locus and a second one containing genotyping data. These two files were then read by the population program (Stacks v. 2.53) to generate a Variant Calling Format (VCF) file. To be processed, a locus had to be present in at least three populations and at least in 50% of all individuals.

SNPs were further filtered for quality using the R package vcfR (Knaus & Grünwald, 2017). First, violin plots were created to inspect the distribution of sequencing depth for variants from each isolate. Isolates with average sequencing depth lower than 4 were excluded from downstream analysis, and a 95% confidence interval of the depth distribution was created to remove variants with unusually high sequencing depth. Variants with more than 15% and isolates with more than 55% missing data were also removed from the data set. A Manhattan plot based on Hedrick's G′ST statistics was built using the genetic distance function available in vcfR to visualize the distribution of the variant's differentiation among the filtered data sets. For each isolate, heterozygosity was explored by plotting a histogram showing the frequency distribution of allele balance for the first and second most abundant allele, as described by Knaus et al. (2020). The VCF data were converted into a format (.bed, .fam, and .bim files) readable by the software PLINK v. 1.9 (Purcell et al., 2007) and linkage disequilibrium was calculated as r 2 between SNPs for each contig and 10‐Mb window. The VCF data were further filtered to remove SNPs in perfect linkage disequilibrium (linkage disequilibrium pruning) using an in‐house python script. Finally, uninformative loci were removed using the informloci() function implemented in poppr with a cut‐off value of 2% and a minor allele frequency of 0.01 (Kamvar et al., 2014, 2015).

Genetic differentiation and population structure of P. destructor populations

The data were organized into four data sets. First, the “complete data set” included all the isolates retained after the filtration steps. The second one, the “Québec data set,” was obtained by removing the Californian isolates from the complete data set. The third data set, called the “regional data set,” corresponded to isolates obtained by random sampling in the three muck soil basins when agronomists and crop specialists reported symptoms. Finally, the fourth data set, referred to as the “targeted data set,” corresponded to isolates obtained in 2019 by sampling in the three severely infested fields. It must be noted that the isolates in the regional and targeted data sets do not overlap.

To determine the genetic structure among P. destructor populations, the AMOVA method was applied using poppr.amova() from the R package ade4 (Dray & Dufour, 2007). AMOVAs were performed separately with the regional data set and with the targeted data set, with significance assessed by performing 999 permutations for each AMOVA. In addition, for both data sets, pairwise AMOVA was performed to calculate the pairwise index of population differentiation . These analyses were conducted with the clone‐corrected and nonclone‐corrected data, with similar results for both; therefore, only the noncorrected data were retained (data not shown).

DAPC (Jombart et al., 2010) was performed on the complete data set and the Québec data set. Because there is a risk of missing useful information or overfitting the data if too few or too many PCs are retained, the number of PCs was chosen using the function xvalDAPC() from the R package adegenet (Jombart et al., 2010). Cross‐validation was performed with a training set corresponding to 90% of the data and a validation data set made up of the remaining 10%. The procedure was repeated 1000 times at each step of the PC retention. The retained number of PCs should provide the highest HMS and the lowest RMSE, with greater importance given to the lowest RMSE. A dendrogram based on the UPGMA was built in poppr by grouping populations using Nei's genetic distance.

To test for the underlying population structure, k‐means clustering was performed with find.clusters() from the R package adegenet (Jombart, 2008) to group isolates in clusters ranging from 1 to 10 clusters. The value of k that returned the most parsimonious fit based on the BIC corresponded to the so‐called theoretical number of genetic clusters. DAPC was then used with no a priori grouping of isolates into genetic clusters. A UPGMA dendrogram was built by grouping individual isolates based on Provesti's genetic distance in poppr.

The standardized index of association was computed to estimate the degree of linkage within each population defined a priori or not using the R package poppr (Kamvar et al., 2014, 2015). To test the null hypothesis of random mating, the statistical assessment of estimates was performed with 999 permutations; the rejection of this hypothesis is generally considered to be associated with an asexual mode of reproduction (Agapow & Burt, 2001).

The Mantel test was used with the regional and targeted data sets to analyse the relationship between genetic divergence and geographic distance (Mantel, 1967). For each of the two Mantel tests, two matrices were created: a matrix of genetic distances calculated with dist.diss() in poppr (Kamvar et al., 2014) and a matrix of geographic distances obtained with dist() in ade4 (Dray & Dufour, 2007). The Mantel tests themselves were performed with the function mantel.randtest() available in ade4, with p values computed from 9999 permutations (Dray & Dufour, 2007). A Mantel correlogram was produced for the Québec data set to test for spatial autocorrelation between pairs of isolates at seven distance classes ranging from 1500 to 15,000 m. The analysis was conducted as described by Legendre and Legendre (2012) using the function mantel.correlog() available in the R package vegan (Oksanen et al., 2016). For each Mantel correlogram ordinate, the Pearson correlation coefficient was computed and tested for significance with 9999 permutations.

FIGURE S1 Violin plots showing the distribution of the total number of reads per sample for the two Ion Proton runs (two chips per run) performed in this study

Click here for additional data file.

FIGURE S2 Distribution of (a) the number of raw reads per sample, (b) the number of reads kept after alignment to the reference genome, and (c) the average sequencing depth per sample after alignment to the reference genome

Click here for additional data file.

FIGURE S3 Violin plots showing the distribution of sequencing depth among the loci for each isolate after the filtering steps

Click here for additional data file.

FIGURE S4 Example of the distribution of allele balance at heterozygous positions for selected Peronospora destructor isolates

Click here for additional data file.

FIGURE S5 Heat map showing the mean sequencing depth for each sample (columns) and for each of the 1,340 retained loci (lines) after the linkage disequilibrium filtering step

Click here for additional data file.

FIGURE S6 UPGMA dendrogram for the Québec isolates based on Provesti's genetic distance. Isolate colours correspond to the genetic clusters identified by discriminant analysis of principal components (DAPC)

Click here for additional data file.

FIGURE S7 Relationship between genetic and geographic distances (Mantel test) between Peronospora destructor isolates for (a) the Québec subset and (b) the targeted subset

Click here for additional data file.

TABLE S1 Results of the analysis of molecular variance (AMOVA) conducted with Peronospora destructor isolates from the arbitrary sampling carried out in Québec (regional data set)

Click here for additional data file.

TABLE S2 Results of the analysis of molecular variance (AMOVA) conducted with Peronospora destructor isolates from the arbitrary sampling carried out in Québec (regional data set)

Click here for additional data file.