Genetic diversity, QoI fungicide resistance, and mating type distribution of Cercospora sojina—Implications for the disease dynamics of frogeye leaf spot on soybean

Frogeye leaf spot (FLS), caused by Cercospora sojina, causes significant damage to soybean in the U.S. One control strategy is the use of quinone outside inhibitor (QoI) fungicides. QoI resistant isolates were first reported in Tennessee (TN) in 2010. To investigate the disease dynamics of C. sojina, we collected 437 C. sojina isolates in 2015 from Jackson and Milan, TN and used 40 historical isolates collected from 2006–2009 from TN and ten additional states for comparison. A subset of 186 isolates, including historical isolates, were genotyped for 49 single nucleotide polymorphism (SNP) markers and the QoI resistance locus, revealing 35 unique genotypes. The genotypes clustered into three groups with two groups containing only sensitive isolates and the remaining group containing all resistant isolates and a dominant clonal lineage of 130 isolates. All 477 C. sojina isolates were genotyped for the QoI locus revealing 344 resistant and 133 sensitive isolates. All isolates collected prior to 2015 were QoI sensitive. Both mating type alleles (MAT1-1-1 and MAT1-2) were found in Jackson and Milan, TN and recovered from single lesions suggesting sexual recombination may play a role in the epidemiology of field populations. Analysis of C. sojina isolates using SNP markers proved useful to investigate population diversity and to elaborate on diversity as it relates to QoI resistance and mating type.


Introduction
Frogeye leaf spot (FLS) of soybean (Glycine max Merr.), caused by the fungal pathogen C. sojina Hara, was first identified in Japan in 1915 and South Carolina, United States in 1924 [1,2]. FLS is an important foliar disease of soybean in the US and symptoms can appear on stems, pods, and seeds [3]. Initial symptoms appear as small, light brown circular spots which develop a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 a darkish brown to reddish margin [4]. As the infected leaf area reaches 50%, the leaves blight, wither and fall prematurely. Conidia, produced on conidiophore are primary and secondary sources of inoculum and are produced on infected leaves, stems, and pods [2]. Warm temperature and frequent rainfall promote severe disease development and fully expanded leaves have smaller lesions compared to younger leaves [2].
According to the food and agriculture organization (FAO), the United States is the world's leading producer of soybean (second only to corn) and produced 106.9 million metric tons of soybeans in 2015 [5]. Yield losses are due to reduced photosynthetic area and the premature defoliation of leaves [3,6]. In the US, FLS has been predominately a disease of southern warm and humid regions [3,7], but has been reported in northern states, such as Iowa in 1999, Wisconsin in 2000 [8] and Ohio in 2006 [9]. Damage depends on cultivar and location, with yield losses ranging from 10% to more than 60% [4,[10][11][12].
FLS is a polycyclic disease and remains active throughout the growing season [11,13]. Spores are dispersed by the wind and water splashes [11]. C. sojina can overwinter in plant debris and may survive for two years in northern environments [9,14]. Greenhouse trials suggest C. sojina may survive on alternate hosts in the absence of soybean [15]. Control is accomplished through the use of cultural practices, fungicides and planting of resistant cultivars. Cultivars with genetic resistance to FLS can be effective and three resistance genes, Rcs (Resistant to C. sojina), have been identified including Rcs1 [16], Rcs2 [17] and Rcs3 [18]. The Rcs3 gene confers resistance against race 5 and all known races of C. sojina present in the U.S. [18,19]. Crop rotation for two years is thought to reduce infection and limit disease severity [14,20]. In addition, the use of pathogen-free seeds and application of fungicides are used to decrease disease severity [20]. While these management practices can still be successful, they have placed selection pressures on C. sojina populations and resulted in isolates that have overcome genetic resistance, namely the Rcs1 and Rcs2 genes [16,17,21]. Similarly, C. sojina isolates have developed resistance to the quinone outside inhibitor (QoI) fungicide group [22].
Previous studies characterized population diversity of C. sojina using AFLP and simple sequence repeat (SSR) markers and there is evidence for sexual outcrossing in field populations [13,23]. Currently, there are no universally accepted soybean differentials and reports of race diversity include 12 races of C. sojina in the US, 22 races in Brazil and 14 races in China [20]. Use of the same 12 soybean differentials produced differing numbers of proposed races in two separate studies [3,9].
Our primary objective was to investigate the dynamics of QoI resistance on eight cultivars of soybean that were fungicide treated and untreated at two locations in Tennessee in 2015. Supporting objectives included the development of novel SNP markers to characterize genotypic diversity, comparison of the 2015 isolates to isolates collected previously and from surrounding states, and an assessment of mating type loci at Tennessee locations and within individual lesions.

Materials and methods
Sample collection, single-lesion isolation, and DNA extraction This study was carried out on private land and the owner gave permission to collect all samples. In 2015, soybean leaves with typical FLS lesions were collected from research plots at two locations in Tennessee (Milan and Jackson). In total, 437 isolates, 203 from Jackson and 234 from Milan, were collected from eight fungicide treated and non-treated Maturity group III soybean cultivars (Table 1). Cultivars were planted in four rows with 76.2 cm row spacing (30-inch row spacing) in 7.6 m (25 ft) long plots in a randomized complete block design with four replications. Plots were split. Two rows were untreated, and two rows were treated at the R3 growth stage (beginning pod) with Quadris Top SB at 8 fl oz/a (0.12 kg a.i./ha Azoxystrobin and 0.07 kg a.i./ha Difenoconazole, Syngenta Corp., Basel, Switzerland). A single isolate of C. sojina was obtained from a single lesion from an individual leaflet. Sporulation was induced by incubating leaves in a plastic bag with moist towels at room temperature (approximately 24˚C/ 72˚F). Spores were harvested with a flame-sterilized needle using a dissecting microscope and 8-10 spores transferred to RA-V8 agar media (rifampicin 25 ppm, ampicillin 100 ppm, 160 mL unfiltered V8 juice, 3 gm calcium carbonate and 840 mL water). Observations were made daily and contaminated sectors removed. After seven days, single-lesion isolates of C. sojina were transferred to new plates. In addition, a set of 40 isolates from 10 different states, collected between 2006 and 2009, were included (Table 2).
For DNA extraction, single-leision isolates were grown in 24-deepwell plates (Fisher Scientific) with 1 mL RA-V8 liquid broth (same as above, minus the agar) per well. DNA was extracted as described by Lamour and Finley (2006). Briefly, this includes harvesting mycelium from the broth cultures into a 96-well 2 mL deepwell plate pre-loaded with 3-5 sterilized 3 mm glass beads. The plates are freeze dried and the dried mycelium powdered using a Mixer-Mill bead beating device (Qiagen). The powdered mycelium was then lysed and a standard glass fiber spin-column DNA extraction completed. The resulting genomic DNA was visualized on a 1% gel and quantified using a Qubit device (Thermo Fisher Scientific Inc.) using the high-sensitivity assay.

SNP marker discovery and Targeted-sequencing based genotyping
Whole genome sequencing was accomplished for three FLS isolates, two from the collection of Dr. Dan Philips at the Univesity of Georgia and one from Tennessee. These include isolate FLS19 (TN10) recovered from the Georgia Experiment Station in 1994, FLS21 (TN85) recovered from Charleston Mississippi in 1994, and FLS11 (CS10117) recovered from Milan, Tennessee in 2010. Genomic DNA was extracted from freeze-dried and powdered mycelium using a standard phenol-chloroform approach and the resulting DNA was submitted to the Beijing Genomics Institute in China for 2x100 paired-end sequencing on an Illumina HiSeq2000 device. De novo assembly, read mapping and SNP discovery were accomplished with CLC Genomics workbench 7 (Qiagen). As there was no public reference genome available at the time, FLS21 was de novo assembled using the default settings in CLC and the resulting contigs used as a reference genome. All open reading frames (ORFs) longer than 300 amino acids were predicted using CLC and annotated onto the FLS21 contigs. The raw reads from FLS11 and FLS19 were then mapped to the draft reference (separately), and putative single nucleotide variants (SNVs) identified at sites with at least 20X coverage and an alternate allele frequency greater than 90%. A subset of the SNVs were chosen from the largest contigs for further genotyping using a targeted sequencing approach. Custom Perl scripts (www.github.com/sandeshsth) were used to extract the flanking sequences for the panel of SNPs and primers were designed using Batch-Primer3 v1.0 (http://probes.pw.usda.gov/batchprimer3/) to amplify targets between 80 and 120bp in length. Primers for 50 SNPs, including the mitochondrial QoI resistance locus, are summarized in Table 3. Primer sequences and genomic DNA were sent to Floodlight Genomics (Knoxville, TN) for processing as part of a non-profit Educational and Research Outreach Program (EROP) that provides targeted-sequencing services at cost for academic researchers. Floodlight Genomics uses an optimized Hi-Plex approach to amplify targets in multiplex PCR reactions and then prepares libraries for sequencing on an NGS device [24]. The services include testing of primers to determine optimal multiplex mixtures followed by standard PCR amplification that includes the addition of sample-specific barcode sequences. The resulting amplicons are prepared for sequencing using PCR-free library construction and the samples presented here were sequenced on an Illumina HiSeq3000 device running a 2x150bp configuration per the manufacturer's directions (Illumina). The resulting sequences were demultiplexed based on the sample-specific barcodes and then mapped to the reference contigs and genotypes assigned for loci with at least 6X coverage and an alternate allele frequency greater than 90%.

QoI resistant locus genotyping and analysis
A single nucleotide polymorphism (G/C) in the Cytochrome b gene at position 143 of the C. sojina mitochondrial genome confers resistance to QoI fungicides. A custom TaqMan SNP genotyping assay was designed using the online design tools from Applied Biosystems (Thermo Scientific) and include the forward primer GGGTTATGTTTTACCTTACGGACAAATG and reverse primer GTCCTACTCATGGTATTGCACTCA and two probes to discriminate resistant and sensitive isolates: ACTGTGGCAGCTCATAA with VIC for the "G" resistance allele and ACTGTGGCACCTCATAA with FAM for the "C" sensitive allele [21]. Each 5 μl reactions consisted of 2.5 μl 2X Taqman master mix, 0.25 μl of 20X Assay working stock (primers + probes) and 2.25 μl of genomic DNA or water as negative control. Reactions were conducted in a Table 3. Summary data for 49 nuclear SNP loci and the mitochondrial QoI-resistance locus.

Genetic analysis
Only SNP loci with no missing data were retained for analyses. Isolates with identical multilocus genotype were considered members of a clonal lineage and assigned a unique identifier (G1-G35) ( Table 4). All subsequent analyses were conducted on the clone-corrected data except for analysis of molecular variance (AMOVA). For each marker, the allele frequencies, effective number of alleles, Nei's gene diversity, and Shannon's Information index were measured using POPGENE [26]. Population structure was assessed using Bayesian clustering in Structure 2.3.4 [27] with the following settings: no prior population information, admixture model, and allele frequency correlated. Structure was run for each possible cluster (K = 1 to 35) with 30 replications or independent runs at 1,000,000 burn-in followed by 1,000,000 Markov Chain Monte Carlo simulations. Structure Harvester was used to find the most optimal value of K (using Evanno's method) from the results obtained from Structure [28,29]. Distance based multivariate analysis (multiple loci and multiple genotypes) was performed to observe the relationship among different genotypes. A pairwise genetic distance matrix was computed and used for covariance-standardization and principle coordinate analysis (PCoA) in GENALEX 6.502 [30]. A phylogenetic tree of the unique genotypes was constructed using the maximum likelihood method based on the General Time Reversible model with 1000 bootstrap replications using Mega 6.06 [31,32]. Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 4.4125)). Minimum spanning networks [33] were constructed using PopART at epsilon zero (http://popart.otago.ac.nz/). This network shows the relationship among closely related intraspecific individuals with alternative potential evolutionary paths in the form of cycles. The partitioning of genetic variance was performed with AMOVA using GENALEX 6.502. The AMOVA analysis was used to estimate molecular variance for four different groupings of the isolates including: JTN and MTN from 2015; isolates from 2015 and historical TN isolates;    (Table 3).

Population structure
Clone-correction of the 186 multi-locus genotypes revealed 35 unique genotypes ( Table 4). Most of the samples (130 of 186) from Jackson and Milan in TN were clonal and belonged to  genotype G22, although some isolates belonged to other genotypes ( Table 4). The discriminatory power of these SNPs for identifying clonal lineages is illustrated in Fig 1. In total, 35 unique genotypes were captured by 15 SNPs. Diversity statistics for the 50 SNP markers are presented in Table 5. Each SNP locus had two alleles. The minor allele frequency (MAF) ranged from 0.057 to 0.486, with 88% of the markers having MAF > 0.1. The effective number of alleles ranged from 1.121 to 1.998, with a mean of 1.653. Gene diversity, the probability that two randomly chosen alleles from the population are different, ranged from 0.108 to 0.5, with a mean of 0.377. Shannon's information index varied from 0.219 to 0.693, with average of 0.557. Bayesian clustering of the 35 unique genotypes predicted three as the most probable K as shown in Fig 2. Grouping of the 35 unique genotypes into three clusters is shown in Fig 3A. A phylogenetic tree constructed using the maximum likelihood approach showed three distinct clades identical to the three clusters generated by Structure (Fig 3B). The three groups are also supported by the principle coordinate analysis (Fig 4). The Minimum spanning network provides an alternative way to visualize how the 186 samples are distributed among the 35 unique genotypes. Fig 5A is colored to indicate sample locations. Fig 5B is the same network colored to indicate the distribution of QoI resistant and sensitive isolates among the 35 unique genotypes. All 130 isolates in G22 clonal lineage were resistant with three other genotypes having resistant isolates: G08, G23, and G30. The remaining isolates fell into 31 genotypes, and all were QoI sensitive.  (Table 6). Similarly, pairwise Nei genetic identity ranged from 0.534-0.999 with maximum identity between Jackson and Milan (Table 6). AMOVA showed high genetic variance within the population and no genetic variance between JTN and MTN (Table 7). However, AMOVA showed 64.8% variance between isolates from Tennessee in 2015 and 2007, which should be assessed with caution due to the low number of isolates analyzed from 2007. Overall, the seven locations had a variance of 40.6% among populations. PhiPT analyses also indicates high population differentiation for isolates collected between 2007 and 2015 (Table 8). Fungicide treated and untreated populations had a non-significant PhiPT value (p = 0.025) and only 3% variation.

QoI resistant isolates
The Taqman assay successfully discriminated the resistant and sensitive alleles present in the mitochondrial Cytochrome b gene based on our positive controls. The number of QoI resistant and sensitive isolates of C. sojina recovered from treated and non-treated cultivars from Jackson and Milan are given in Table 9. Resistant isolates dominate the collection from Tennessee in 2015 (344 out of the 437 isolates tested) while all historical isolates were sensitive ( Fig  6A and 6B). The Chi-square Fisher's exact two tailed test indicates a significant increase of resistant isolates in the field in 2015 compared to isolates collected prior to 2015 (P < 0.0001). Jackson and Milan were dominated by resistant isolates, although Milan had a significantly greater proportion compared to Jackson (85% vs. 72%, respectively p<0.0001) (Fig 6C). The  number of resistant isolates recovered from fungicide treated cultivars was significantly greater (91% resistant, 191 out of 209 isolates) than those collected from non-treated cultivars (67% resistant, 153 out of 228 isolatesz) (P < 0.0001) (Fig 6D). Excluding cultivar C3-VAR Beck's 393R4 (due to only 3 resistant isolates being recovered from it), the percentage of resistant isolates ranged from 68 to 90% and was not significantly different across cultivars (Fig 6E).

Mating-type distribution
Both mating types (Mat1-1-1 and Mat1-2) were indentified (Table 10). Mating types of representative samples from the 35 unique genotypes are given in Table 4. Both mating types were identified in all three of the genetic groups described and were present in single-lesion isolates.

Discussion
Our primary objective was to investigate the dynamics of FLS QoI resistance on eight cultivars, fungicide treated and untreated, at two locations in Tennessee in 2015 and to gain perspective by comparison to isolates from previous years and surrounding states. The development of novel SNP markers proved useful and it appears a subset of 15 markers was sufficient to characterize the overall genotypic diversity. Tennessee isolates from 2015 were dominated by a single QoI resistant clone, found on both fungicide treated and untreated cultivars. Although a greater proportion of resistant isolates were recovered from fungicide treated plants, there were no differences in the recovery of resistant vs. sensitive isolates among seven cultivars. Interestingly, a combination fungicide (e.g. Quadris Top SB) containing active ingredients from two different fungicide groups (Azoxystrobin from the QoI fungicide group and Difenoconazole from the DeMethylation Inhibitor (DMI) fungicide group) resulted in a higher proportion of resistant isolates, suggesting fungicide mixes can still exert selection pressure for QoI resistance. This may explain why C. sojina populations in TN have continued to increase since the first report in 2010; even when producers are using fungicides that contain active ingredients in different (non-QoI) fungicide groups (unpublished data).
Although the high levels of genetic diversity reported within fields for populations of FLS in Arkansas were not found at our two locations [13], this may be because our sampling was conducted late in the season, allowing FLS an extended period for polycyclic reproduction and dissemination. It will be useful to conduct sampling at earlier time points in the growing season to capture the full complement of genetic diversity within fields [13,34]. The finding of both mating types in individual lesions and across all three genetic groups suggests ample opportunity for sexual recombination to play a role in shaping the population structure and to transfer the maternally inherited QoI resistance into diverse nuclear genetic backgrounds. Overall, the low level of genetic differentiation among sites was similar to the findings with C. zeina, which also had a lack of regional population differentiation [35].
Interestingly, Tennessee isolates from 2007 and 2015 shared four unique genotypes, including a second dominant genotype differing from the modern dominant clonal lineage by two SNP markers. This suggests clonal lineages may survive many years. This is a reasonable conclusion considering our data and previous studies suggesting C. sojina can remain viable in plant residues for more than two years [9,14].  Most soybean cultivars, fungicide treated or untreated, had a high frequency of resistant isolates compared to the sensitive isolates. However, only three isolates were recovered from the FLS resistant cultivar VAR Beck's 393R4 and only from untreated plants. This suggests deployment of resistant cultivars can reduce FLS development and the further spread of QoI resistant C. sojina isolates and warrants further investigation.
Although a limited number of historical isolates were available for analyses and the population analyses are correspondingly weak, these isolates provide a useful perspective on the current situation in Tennessee. Clearly there has been a shift in the proportion of QoI resistant isolates in the last decade. Additional fine-scale sampling over time and over a wider geographical area will be useful to measure the longevity of the dominant QoI resistant clonal lineage (G22) and to investigate the role outcrossing and sexual recombination may play in the epidemiology of FLS in the face of multiple selection pressures.

Author Contributions
Conceptualization: SS AM HY.