Population Genetic Analyses of Botrytis cinerea Isolates From Michigan Vineyards Using a High-Throughput Marker System Approach

As sequencing costs continue to decrease, new tools are being developed for assessing pathogen diversity and population structure. Traditional marker types, such as microsatellites, are often more cost effective than single-nucleotide polymorphism (SNP) panels when working with small numbers of individuals, but may not allow for fine scale evaluation of low or moderate structure in populations. Botrytis cinerea is a necrotrophic plant pathogen with high genetic variability that can infect more than 200 plant species worldwide. A panel of 52 amplicons were sequenced for 82 isolates collected from four Michigan vineyards representing 2 years of collection and varying fungicide resistance. A panel of nine microsatellite markers previously described was also tested across 74 isolates from the same population. A microsatellite and SNP marker analysis of B. cinerea populations was performed to assess the genetic diversity and population structure of Michigan vineyards, and the results from both marker types were compared. Both methods were able to detect population structure associated with resistance to the individual fungicides thiabendazole and boscalid, and multiple fungicide resistance (MFR). Microsatellites were also able to differentiate population structure associated with another fungicide, fluopyram, while SNPs were able to additionally differentiate structure based on year. For both methods, AMOVA results were similar, with microsatellite results explaining a smaller portion of the variation compared with the SNP results. The SNP-based markers presented here were able to successfully differentiate population structure similar to microsatellite results. These SNP markers represent new tools to discriminate B. cinerea isolates within closely related populations using multiple targeted sequences.

As sequencing costs continue to decrease, new tools are being developed for assessing pathogen diversity and population structure. Traditional marker types, such as microsatellites, are often more cost effective than single-nucleotide polymorphism (SNP) panels when working with small numbers of individuals, but may not allow for fine scale evaluation of low or moderate structure in populations. Botrytis cinerea is a necrotrophic plant pathogen with high genetic variability that can infect more than 200 plant species worldwide. A panel of 52 amplicons were sequenced for 82 isolates collected from four Michigan vineyards representing 2 years of collection and varying fungicide resistance. A panel of nine microsatellite markers previously described was also tested across 74 isolates from the same population. A microsatellite and SNP marker analysis of B. cinerea populations was performed to assess the genetic diversity and population structure of Michigan vineyards, and the results from both marker types were compared. Both methods were able to detect population structure associated with resistance to the individual fungicides thiabendazole and boscalid, and multiple fungicide resistance (MFR). Microsatellites were also able to differentiate population structure associated with another fungicide, fluopyram, while SNPs were able to additionally differentiate structure based on year. For both methods, AMOVA results were similar, with microsatellite results explaining a smaller portion of the variation compared with the SNP results. The SNP-based markers presented here were able to successfully differentiate population structure similar to microsatellite results. These SNP markers represent new tools to discriminate B. cinerea isolates within closely related populations using multiple targeted sequences.
The genus Botrytis is highly genetically diverse with more than 30 species that differ in morphology, ecology, biology, and host range (Walker, 2016). The genus Botrytis was generally considered as a single complex species until the late 1990s when it was subdivided into two clades, one clade contains Botrytis spp. that infect mostly monocots and some dicots, while the second clade contains Botrytis spp. that infect a wide host range of eudicots; B. cinerea falls under this second clade (Staats et al., 2005). Population structure and genetic variations have been studied in B. cinerea populations, and new morphologically identical or similar species were identified. Recently a number of cryptic species causing gray mold that lives in sympatry with the B. cinerea complex have been identified on a variety of hosts (Li et al., 2012;Saito et al., 2016;Dowling et al., 2017;Rupp et al., 2017b;Harper et al., 2019). These new cryptic species are more likely considered as host or region specific. Formerly known as B. cinerea Group I, B. pseudocinerea isolates are morphologically identical and live in sympatry with B. cinerea (Fournier et al., 2005;Walker et al., 2011). Genetic polymorphisms in transposable element presence and a group of genes including Bc-hch, erg27, sdh, and cyp51 between B. cinerea and B. pseudocinerea provided evidence for the differentiation of B. pseudocinerea as a new species (Fournier et al., 2003;Walker et al., 2011;Plesken et al., 2015). Within B. cinerea sensu stricto, a large variability in genetic and phenotypic diversity, and host specialization have also been observed (Corwin et al., 2016;Mercier et al., 2019;Soltis et al., 2019;Meng et al., 2020).
Generally, the population structure in B. cinerea was detected to vary between different hosts (Fournier and Giraud, 2008;Walker et al., 2015) or year (Walker et al., 2015;Delong et al., 2020), while less or no variation was detected at the region level (Fournier and Giraud, 2008;Karchani-Balma et al., 2008;Esterio et al., 2011;Wessels et al., 2013;Walker et al., 2015). On the other hand, Muñoz and Moret (2010) found that genetic diversity was high, and population structure varied when comparing isolates of B. cinerea at a continent level. The molecular marker method that is used to investigate genetic differences can also contribute to the observed genetic differences (Walker et al., 2011).
Various molecular markers have been used to investigate the genetic variability and population structure in B. cinerea. Fournier et al. (2003) developed a PCR-RFLP method that differentiated B. pseudocinerea from the B. cinerea complex based on the polymorphism detected in the Bc-hch gene. Other genes that are Botrytis spp. specific were used to differentiate different populations on table grape and blueberry using Sanger sequencing (Staats et al., 2005(Staats et al., , 2007Saito et al., 2016). Recently, DeLong et al. developed a set of microsatellite markers spanning the genome to characterize Botrytis populations (2020). Commonly, the presence/absence of TEs is used for the determination of Botrytis sp. population diversity (Kretschmer and Hahn, 2008;Fekete et al., 2012;Wessels et al., 2016;Hu et al., 2018). Previous studies that used TEs or microsatellites as markers for population differentiation showed that different hosts can be dominated by different populations such as grape and pomegranate, and were dominated by transposa isolates that have both TEs (Váczy et al., 2008;Johnston et al., 2014;Delong et al., 2020;Testempasis et al., 2020). The fungicide resistance profile is different between different strains of Botrytis (Martínez et al., 2008;Leroux et al., 2010;Johnston et al., 2014;Delong et al., 2020). However, studies of population variability in relation to fungicide resistance profile showed limited to no association with the population structure (Wessels et al., 2016;Campia et al., 2017;Hu et al., 2018;Delong et al., 2020).
Several studies in other systems have compared the use of single-nucleotide polymorphism (SNP) or sequencing-based markers and microsatellite markers to describe population structure and diversity (Fischer et al., 2017;Lemopoulos et al., 2018). In one such study, SNP diversity estimates and microsatellite heterozygosity in Arabidopsis were not significantly correlated, but genetic differentiation among populations was correlated (Fischer et al., 2017). Similarly in trout, 16 microsatellites performed similarly to >4,000 SNPs at measuring genetic differentiation, but SNPs were more accurate at estimating individual level heterozygosity (Lemopoulos et al., 2018).
Frontiers in Microbiology | www.frontiersin.org However, similar studies have found that studies, where a low number of SNPs >300 were used, had lower power than microsatellites (Vali et al., 2008). Yet most of these studies involved animals or plants with larger plant genomes than Botrytis (Ciani et al., 2013;Singh et al., 2013). A study on Plasmodium vivax, a malarial parasite, demonstrated that 146 high-quality SNPs using an amplicon sequencing approach were more informative than microsatellite markers (Fola et al., 2020).
Because different Botrytis sp. can exhibit differences in fungicide resistance profile, it is critical to understand the pathogen population structure in different environments. This will allow development of better disease management schemes.
There is no information about B. cinerea population structure and genetic diversity in Michigan grapevine. Therefore, this is the first study to investigate the genetic diversity of the B. cinerea population in Michigan grapes in combination with fungicide sensitivity phenotypic characteristics. Our objectives were to study the population structure of MI isolates of B. cinerea related to year, location, fungicide resistance, and to compare the use of microsatellite and amplicon-based sequencing SNP strategies to quantify genetic diversity and population structure.

Isolates and DNA Extractions
A total of 82 B. cinerea isolates were collected from four Michigan vineyards, three southwest and one northwest Michigan locations (Supplementary Table 1). The "West" Michigan vineyard location represents B. cinerea isolates recovered from a Vitis interspecific hybrid (cv. Vignoles), and "Southwest" Michigan vineyards 1 and 2 are also hybrids recovered from cultivars Vignoles and Aurora, respectively. Finally, the "Northwest" Michigan samples were recovered from symptomatic Vitis vinifera (cv. Reisling).
These locations were sampled in both 2014 and 2018, where 42 and 40 samples were collected in 2014 and 2018, respectively, with at least eight isolates per location. Throughout the study, isolates were maintained on 20% clarified V8 agar media (100 ml of V8 juice, 1 g of CaCO 3 , 10 g of agar, and 400 ml of distilled water) at room temperature. All isolates were previously evaluated for fungicide resistance against eight fungicides with seven different chemical classes (Alzohairy et al., 2020). For each isolate, a multiple fungicide resistance (MFR) value was determined based on the total number of fungicides an isolate demonstrated resistance toward. For DNA extraction and quantification, we followed Alzohairy et al. (2020); in brief, mycelia were collected from 1-to 2-week-old cultures. Mycelia were lyophilized, then approximately 5 mg of tissue was ground using a tissuelyser (Qiagen, Valencia, CA). Automated DNA extraction was performed using a MagMax plant DNA isolation kit (ThermoFisher 192 Scientific, Waltham, MA) and processed on the KingFisher Flex purification system (ThermoFisher Scientific). For DNA quantification, we used two methods: first, DNA samples were processed with a Qubit 1X dsDNA HS assay kit (ThermoFisher Scientific), then DNA concentrations were quantified by a Qubit 4 fluorometer (ThermoFisher Scientific); second, DNA samples were processed using PicoGreen Quant-iT DNA reagent and kits (ThermoFisher Scientific), then DNA concentrations were determined using the Synergy HTX multi-mode reader (Biotec, VT). Insufficient DNA was available for all 82 isolates for both microsatellite and SNP evaluation, and only 74 isolates were consistent between the two datasets.

Design of Multiplex Primer Sets
Primers (microsatellite and candidate genes) were grouped in sets of five, six, or eight-plex for simultaneous amplification of Botrytis target genes. Grouping was based on the following parameters: (1) expected amplicon size, with at least 30-bp difference between each amplicon for clear validation by gel electrophoresis and (2) the same annealing temperature. Since the presence of multiple oligonucleotides in one PCR reaction could alter the efficiency of amplification, several combinations of primer concentrations (0.2-0.4 μM) were tested in parallel in a single-plex and in multiplex formats using the cycling conditions mentioned above. Amplicons were separated by electrophoresis at 50 V for 2 h on ethidium bromide-stained (10 μg/ml) 4% agarose gel and visualized under the UV light for validation of band presence.

Single-Nucleotide Polymorphism Processing and Alignments
Quality control, processing, read alignments, and SNP calling were completed using the Galaxy bioinformatics server. 1 Reads were trimmed using Trimmomatic (v0.38.0) with paired-end adapter trimming for Illumina MiSeq and HiSeq adapters. Reads shorter than 30 bp were removed, and read quality was assessed using a sliding window of 4 bp with a quality ≥20. Overrepresented sequences, GC content, and read quality were visualized using FastQC. Single reads (reads without a corresponding mate) were not retained for downstream analyses. Mapping was performed using Bowtie2, read groups (Sam/Bam format) were automatically set, and alignments were based on the sensitive local parameters to the B. cinerea reference genome (Van Kan et al., 2017). Unaligned reads were output to a separate file. Sam files were merged for downstream analysis. Single-nucleotide polymorphisms (SNPs) and insertions/deletions (indels) were identified for each isolate using FreeBayes. Indels were left aligned, and 60% of the observations were required for an alternate allele to be suggested within an individual. The variant call format (VCF) file was further filtered using TASSEL 5 (v20180222) to remove sites present in fewer than 10 individuals. The resulting VCF with 496 SNPs was used for all downstream analyses. An SNP was considered to be associated with a primer set, if the variant fell within the mapped boundaries of the forward and reverse primer visualized with Geneious Prime.

Microsatellite Genotyping
Published microsatellites were evaluated across 74 isolates using a three-primer method as previously described (Schuelke, 2000;

Genetic Diversity and Population Structure
For both the SNP and the microsatellite data sets, genetic diversity, and analysis of molecular variance (AMOVA) were calculated for the population based on location, year, and fungicide resistance using Poppr. Insufficient resistant or sensitive isolates were available for fludioxonil or pyraclostrobin, and thus, AMOVAs were not conducted for these fungicides. For the microsatellite data, MFRs for 1, 2, 3, 4, 5, and 6+ were compared due to an insufficient number of isolates in the 6 and 7 categories. Genetic distance trees were generated using a UPGMA tree based on Provesti's distance (Parada-Rojas and Quesada-Ocampo, 2018).

Single-Nucleotide Polymorphism-Based Diversity
The number of reads per isolate ranged from 14,000 to 32,000 and the percentage of reads mapped to the reference genome ranged from 95 to 99.92% (Supplementary Table 3).
Overrepresented sequences ranged from 7 to 11 sequences per isolate. On average, each isolate had 23,844 reads with nine overrepresented sequences and 97% of reads mapping to the reference genome. When BLASTed against the Botrytis genome, unaligned reads were able to be aligned or were determined to be Botrytis but not in the nuclear reference genome (e.g., Cyt b and Mat1-2). SNP distribution across the genome ranged from three SNPs on chromosome 11 to 101 SNPs on chromosome 1 with a total of 496 SNPs detected post filtering used for subsequent analyses ( Table 1). No SNPs were identified on chromosomes 16, 17, or 18, consistent with primer design (Figure 1). The greatest number of positional variants was detected for products associated with primer sets NEPO5, Bc_pop22, alpha 1, Flip3, Bc_pop16, Bc_pop59, Bc_pop82, and mrr1 ( Table 2). Multilocus statistics for the population grouped by year showed that the standardized index of association was significant at 0.0428 and 0.0189 for 2014 and 2018, respectively (Table 3; Figures 2A,B). A greater number of unique multilocus genotypes (MLGs) was detected in 2014 than in 2018, but the total number of MLGs between years was similar. A small (1.6%) yet significant (p = 0.049) difference was detected between isolates collected in 2014 and those collected in 2018 (Table 4)  When grouped by location (West, Southwest 2, Northwest, and Southwest 1), overall G'ST was 0.154. Both the genetic distance tree and AMOVA showed no significant differences based on location among isolates (p > 0.05; Figure 3). The standardized index of association for each location was significant, 0.025-0.074. The highest rbarD was detected at Southwest 1, and Southwest 2 had the highest number of MLGs ( Table 3). The greatest SNP contribution was from position 1,783,010 located on chromosome 1, which was not located within the aligned primer boundaries for any of the primer sets evaluated (Table 2; Supplementary Figure 2). When grouped by MFR, significant variation (p = 0.004, 7.76%) was explained by the populations (Table 5). When grouped by individual fungicide resistances, fluopyram, cyprodinil, and fenhexamid had no significant variation between resistant and sensitive populations. Thiabendazole, iprodione, and boscalid all had moderate  (Stoddart and Taylor, 1988). x Simpson's index (Simpson, 1949). y Nei's unbiased gene diversity (Nei, 1978). Frontiers in Microbiology | www.frontiersin.org 8 April 2021 | Volume 12 | Article 660874 variability (7.7-10%) explained between the resistant and sensitive populations ( Table 5).

Microsatellite-Based Diversity
Alleles identified in the population ranged from four to seven, with an average of 5.89 alleles per locus and a total of 53 alleles across all nine markers ( Table 6). After clone correction, 56 original multilocus genotypes were identified when grouped by year. Both years had a similar number of MLG identified (30 vs. 27), and the standardized index of association (rbarD) for the populations ranged from 0.040 to 0.029 in 2014 and 2018, respectively (Figures 2C,D). Hedrick's GST across all markers was 0.044 (Hedrick, 2005;Meirmans and Hedrick, 2011). AMOVA revealed no significant variation explained by grouping the population by year. When grouped by location, genetic diversity was similar across populations (H exp = 0.61-0.64), but rbarD varied widely with slightly negative values for Southwest 2 and West (−0.003 and −0.002, respectively) and a large positive value (0.15) for Southwest 1 (Table 7). However, AMOVA revealed no significant variation explained when grouped by location (p = 0.069). When grouped by fungicide resistance (individual or MFR), no significant variability was detected between resistant and sensitive populations for cyprodinil, iprodione, or fenhexamid. MFR was grouped into categories designated 3 (0-3), 4, and 5 (5-7) because insufficient numbers of isolates had resistance to 0, 1, 2, 6, or 7 fungicides. AMOVA results indicated significant variability in the population was explained when grouped by MFR (0-7) or MFR categories (3, 4, 5; 3.7% at p = 0.037 and 4.4% at p = 0.012, respectively; Table 8). When grouped into resistant or sensitive categories for thiabendazole or boscalid, 14% of the population variability (p = 0.001) was explained ( Table 8). When grouped by fluopyram resistance or sensitive groupings, 4.4% of the variability was explained (p = 0.026).

Comparison of Single-Nucleotide Polymorphisms and Microsatellite Markers
Both SNP and microsatellite markers were able to identify significant population differentiation and genetic diversity in Botrytis isolates from Michigan. A greater number of MLG groups was identified across the total population and within each location with SNP markers compared with microsatellites (Tables 3 and 7). Both marker systems demonstrated significant linkage among markers and genotypes across both years for rbarD suggesting that populations persist across years (Figure 2). For instance, MFR groupings were significant using both microsatellites and SNPs; however, SNP markers explained a greater proportion of the variation by this grouping than microsatellites (7.76% compared with 3.68%, respectively.) Thiabendazole and boscalid resistance groupings were consistently associated with significant population structure. Population structure associated with fluopyram resistance, a FRAC 7 similar boscalid, was only significant using the microsatellite markers (4.3% at p = 0.026).

DISCUSSION
In this study, we evaluated SNP and microsatellite markers for their ability to describe genetic diversity and population structure in B. cinerea using isolates collected from Michigan vineyards. B. cinerea is a globally distributed pathogen with a broad host range and widespread fungicide resistance. Thousands of studies across the globe have used morphological, genic, microsatellite, and other PCR-based markers to characterize genetic diversity and population structure of Botrytis spp. The pathogen is genetically diverse with closely related morphologically similar species that can be found on the same host or tissue complicate diversity studies (Walker et al., 2011;Li et al., 2012;Saito et al., 2016;Garfinkel et al., 2017;Rupp et al., 2017b;Hu et al., 2018;Harper et al., 2019). While it is widely accepted that population structure exists within Botrytis, the degree and factors by which populations can be differentiated has not been consistent. Depending on the populations and marker systems tested, population structures associated with continent, host, year, and fungicide resistance have all been reported at varying levels of differentiation and significance (Fournier and Giraud, 2008;Muñoz and Moret, 2010;Walker et al., 2015;Delong et al., 2020;Rupp et al., 2017a,b). This variability could be, in part, caused by regional differences in Botrytis population composition. Fungicide spray programs and the resulting resistance have repeatedly been shown to play a large role in Botrytis field population composition, but smaller differences caused by host, season, or marker resolution may also play a role (Wessels et al., 2016;Kozhar et al., 2020;Testempasis et al., 2020). Microsatellites have historically been used to assess genetic diversity and population structure because of their widespread accessibility, requiring no prior known sequence information, highly specialized equipment, or software to analyze. As sequencing technologies have become more affordable, SNP-based assessment of genetic diversity and population structure has become more prevalent (Loera-Sanchez et al., 2019;Sato et al., 2019;Li et al., 2020;Weldon et al., 2020). Using sequencing data has multiple advantages to microsatellites, primarily that no prior sequence information is required, marker transferability is not a concern, and sequence information can provide higher resolution for individuals (Helyar et al., 2011;Du et al., 2019). However, studies have shown that for organisms with larger genomes, well-designed microsatellites are often more effective than small (<400) numbers of SNPs at characterizing differentiation (Vali et al., 2008;Müller et al., 2015;Lemopoulos et al., 2018). However, Botrytis, like many fungal organisms, has a small genome (<45 Mb) and may not require large (>400) numbers of SNPs for accurately differentiating populations (Abbott et al., 2010;Tsykun et al., 2017;Van Kan et al., 2017). This could be achieved either through whole genome sequencing or reduced representation sequencing. However, whole genome sequencing may be cost prohibitive for large numbers of isolates. Amplicon sequencing, compared with GBS, has the added advantage of targeting known regions allowing for comparison of genes of interest for isolates across different studies. As more  Botrytis genomes are sequenced, better target sequences (semiconserved with high diversity across isolates, fungicide resistance genes, etc.) can be selected to improve detection. In this study, we identified 496 SNPs from 52 amplicons and nine microsatellites to characterize 82 and 74 isolates, respectively. Both types of markers were able to identify genetic diversity and population structure across the shared isolates. High levels of clonality were observed between years, suggesting that these populations can be persistent over multiple years in perennial crops. In our study, Botrytis isolates were collected in 2014 and 2018. When comparing genetic diversity between the two marker systems, the 496 SNPs were only able to distinguish an additional two to three MLGs for the locations with the highest number (Southwest 1 and Southwest 2), but with the two lower MLG locations (West and Northwest), more than six additional MLGs were detected. While both types of markers were able to differentiate population structure at the multiple fungicide resistance level, thiabendazole and boscalid resistance, neither marker system was able to detect significant differentiation based on location similar to other Botrytis studies (Hu et al., 2018;Delong et al., 2020;Testempasis et al., 2020). Minor differences associated with region were detected with microsatellite markers, but were not significant at the p = 0.05 level. Significant differentiation based on single fungicide resistances and year of collection were also detected, but differed, between the two marker systems. This could be, in part, due to specific known fungicide resistanceassociated loci being included with the SNP data set and not necessarily in the microsatellites. Yet microsatellite markers were able to differentiate population structure associated with a similar number of fungicides as the SNP markers. Resistance to multiple fungicides and iprodione each explained approximately 7% of the variability observed, but boscalid resistance explained the majority of the differentiation (10% with SNPs and 14% with microsatellites). This is similar to other studies showing that fungicide resistance may be a driving factor in Botrytis population structure in agricultural systems (Kozhar et al., 2020). In summary, microsatellites and SNP markers were both effective at identifying population structure associated with major factors (e.g., fungicide resistance) in Botrytis. However, as populations with greater numbers of individuals are evaluated, SNP markers will likely be more cost effective and useful for identifying closely related species and minor factors associated with population structure.  (Stoddart and Taylor, 1988). x Simpson's index (Simpson, 1949). y Nei's unbiased gene diversity (Nei, 1978

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/biosample/ 18170356).

AUTHOR CONTRIBUTIONS
RN conceived of experiment, analyzed the data, and contributed to writing. JD collected and analyzed the data and contributed to writing. NA and SS collected the data and contributed to the experimental design and writing. TM and SA contributed to the experiment design and writing. All authors contributed to the article and approved the submitted version.

FUNDING
Funding for this project was provided by the California Department of Food and Agriculture Specialty Crop Block grant Project #17-0275-048-SC and USDA ARS appropriate projects 2034-21220-007-00D.