Genome-Wide Associations with Resistance to Bipolaris Leaf Spot (Bipolaris oryzae (Breda de Haan) Shoemaker) in a Northern Switchgrass Population (Panicum virgatum L.)

Switchgrass (Panicum virgatum L.), a northern native perennial grass, suffers from yield reduction from Bipolaris leaf spot caused by Bipolaris oryzae (Breda de Haan) Shoemaker. This study aimed to determine the resistant populations via multiple phenotyping approaches and identify potential resistance genes from genome-wide association studies (GWAS) in the switchgrass northern association panel. The disease resistance was evaluated from both natural (field evaluations in Ithaca, New York and Phillipsburg, Philadelphia) and artificial inoculations (detached leaf and leaf disk assays). The most resistant populations based on a combination of three phenotyping approaches—detached leaf, leaf disk, and mean from two locations—were ‘SW788’, ‘SW806’, ‘SW802’, ‘SW793’, ‘SW781’, ‘SW797’, ‘SW798’, ‘SW803’, ‘SW795’, ‘SW805’. The GWAS from the association panel showed 27 significant SNPs on 12 chromosomes: 1K, 2K, 2N, 3K, 3N, 4N, 5K, 5N, 6N, 7K, 7N, and 9N. These markers accumulatively explained the phenotypic variance of the resistance ranging from 3.28 to 26.52%. Within linkage disequilibrium of 20 kb, these SNP markers linked with the potential resistance genes included the genes encoding for NBS-LRR, PPR, cell-wall related proteins, homeostatic proteins, anti-apoptotic proteins, and ABC transporter.


Introduction
Switchgrass (Panicum virgatum L.) is a perennial biomass crop native to North America. Biomass and other agronomic traits have been the focus of breeding programs [1][2][3][4]. Although many diseases have been reported to cause deleterious effects on yield, research on disease resistance, especially breeding for the resistance, is scarce. Only 47 of the 2328 research articles directly related to diseases in switchgrass [5]. Bipolaris oryzae (Breda de Haan) Shoemaker (teleomorph: Cochliobolus miyabeanus) is one of the major fungi causing Bipolaris leaf spot (BLS) in switchgrass that can reduce biomass by 70% [6]. To manage the disease, breeding is one of the most economical approaches [7].
The natural distribution of switchgrass stretches latitudinally across North America east of the Rocky Mountains. By simple morphological differentiation, switchgrass can be separated into two groups of upland and lowland ecotypes. Upland ecotypes are well adapted to higher latitudes and have higher drought tolerance, whereas lowland ecotypes provide higher yield and require more water [8]. The more in-depth genetic diversity was revealed by the Network-Based single nucleotide polymorphism (SNP) discovery

Phenotyping and Correlation between Traits
The association panel of 478 genotypes from 66 populations was evaluated for the severity of BLS from the mean of Bipolaris lesion percentage from detached leaf assay by vision (DTVI), by image analysis (DTIA), the severity from the leaf disk assay by vision (DSVI), by image analysis (DSIA), mean of two locations (Ithaca, New York and Phillipsburg, Philadelphia), the highest scores between two locations (MTL), mean in New York (NY), the highest score in NY (MNY), mean in Philadelphia (PA), and highest score in PA (MPA) based on population (rank). 'Population' in this study meant a seed source with a specific origin and, in some cases, breeding history. 'Genotype' meant the individual plants propagated from populations.
In the field evaluation based on the mean of the severity score (0 to 5) from two locations, 'Shelter' showed the lowest severity (0.44) and 'SW787' showed the highest severity. However, when considering each location, different populations performed differently. In comparison between NY and PA, Wilcoxson's rank test of the severity showed significant difference between the two locations (p-value = 0.042). The re-ranking incidence across locations indicates the significant GxE effect [32]. In NY, SW803 showed the lowest severity (0.22), and 'SW787' showed the highest severity (3.0). Whereas in PA, 'Shelter' showed the lowest severity (0.33), and 'Pathfinder' showed the highest severity (2.9). Between the two locations, 'SW123', 'SW33', 'SW793', 'SW781', 'High Tide', and 'Timber' had the lowest MTL at 3 (resistant), and 41 populations had the highest MTL at 5 (susceptible). In NY, 'SW123', 'SW128' and 'ECS.6' had the lowest MNY at 2, and 32 populations had the highest MNY at 5. In PA, 'SW115', 'SW802', 'SW31', and 'SW43' had the lowest MPA at 2, and 19 populations had the highest MPA at 5. Although each phenotyping approach cannot determine a single resistant population, a ranking of detached leaf, leaf disk, and mean from two locations showed that 'SW788', 'SW806', 'SW802', 'SW793', 'SW781', 'SW797', 'SW798', 'SW803', 'SW795', 'SW805' were always in the lowest severity group and can be the source of resistance (Table S1). When comparing ecotypes, the field evaluation of means between two locations showed a trend of more resistance in lowland than upland ecotypes significantly ( Figure S1 and Table S1). The difference between lowland and upland ecotypes also occurred in NY and PA. The severities based on genotypes in each phenotyping approach were listed (Table S2).
Broad-sense heritabilities (H 2 ) in different groups of genotypes were similar across groups with some variations (Table 1). For example, H 2 from DSIA in upland was only 0.14 while other groups had around 0.40. Although the severity from different phenotyping approaches in different groups had significantly moderate to high H 2 , H 2 of severity from two locations was zero. This was supported by the significantly different mean of severity from two locations, NY and PA ( Figure S1 and Wilcoxson p-value = 0.042). The different rankings did not only show in the field evaluation, but the other phenotyping approaches also gave different rankings of resistance from DTIA, DSIA, and mean from two locations and correlation among the traits ( Figure S2 and Wilcoxson's p-value < 2.2 × 10 −16 ). For example, from 478 genotypes, 'KY1625_07' ranked 27th from DTIA but ranked 329th in DSIA and 100th in mean from two locations. Such differences led to a low correlation among approaches. The r 2 between DTIA and DSIA was only 0.1 and between DTIA and the mean of two locations (TL) was 0.03 ( Figure 1). In contrast to these low correlations, there were high correlations in DTVI-DTIA (r 2 = 0.84) and DSVI-DSIA (r 2 = 0.85). There was no correlation between BLS resistance and other agronomic or biomass quality traits ( Figure S3). 329th in DSIA and 100th in mean from two locations. Such differences led to a low correlation among approaches. The r 2 between DTIA and DSIA was only 0.1 and between DTIA and the mean of two locations (TL) was 0.03 ( Figure 1). In contrast to these low correlations, there were high correlations in DTVI-DTIA (r 2 = 0.84) and DSVI-DSIA (r 2 = 0.85). There was no correlation between BLS resistance and other agronomic or biomass quality traits ( Figure S3).  , the highest score between two locations (MTL), mean from the two locations (TL), minerals total ash, in vitro dry matter digestibility, cell wall concentration, and ethanol conversion per day from a previous study [33]. , the highest score between two locations (MTL), mean from the two locations (TL), minerals total ash, in vitro dry matter digestibility, cell wall concentration, and ethanol conversion per day from a previous study [33].
Additionally, Qst-Fst comparisons were used to investigate whether the mean difference of each phenotype among populations is greater than the expected mean difference under genetic drift. None of the traits (DTVI, DTIA, DSVI, and DSIA) has Qst-Fst significantly more than zero ( Figure S4), showing no positive selection for resistance to BLS in populations.

Linkage Disequilibrium, Principal Component Analysis and Admixture
To determine the intervals that were potentially linked to the significant markers from GWAS, linkage disequilibrium (LD) decay was estimated by plotting the allele frequency correlations (R2) against the physical distance in base pairs ( Figure 2). In the full set of 478 genotypes, the LD decayed sharply within 20 kb and reached a plateau background within 50 kb. Therefore, in this study, we focused on genes within the 20 kb from the significant SNPs.
Additionally, Qst-Fst comparisons were used to investigate whether the mean difference of each phenotype among populations is greater than the expected mean difference under genetic drift. None of the traits (DTVI, DTIA, DSVI, and DSIA) has Qst-Fst significantly more than zero ( Figure S4), showing no positive selection for resistance to BLS in populations.

Linkage Disequilibrium, Principal Component Analysis and Admixture
To determine the intervals that were potentially linked to the significant markers from GWAS, linkage disequilibrium (LD) decay was estimated by plotting the allele frequency correlations (R2) against the physical distance in base pairs ( Figure 2). In the full set of 478 genotypes, the LD decayed sharply within 20 kb and reached a plateau background within 50 kb. Therefore, in this study, we focused on genes within the 20 kb from the significant SNPs. Based on the principal component (PC) analysis and admixture, there was a population structure in this association panel ( Figure 3a). From PC1 and PC2, upland and lowland ecotypes separated from each group. Moreover, within the lowland ecotype, the latitude of the lowland can be distinguished to lowland north and south groups. The separation of the upland ecotype can also be noticed in PC1 and PC3 that the upland north was grouped apart from the other upland. Besides ecotypes, ploidy levels can be differentiated into groups. Both lowland north and south, and upland north genotypes were tetraploid (4X), while the upland east and upland west genotypes were octaploid (8X). In total, PC1, PC2, and PC3 explained the variance due to the population of 49.35%. The admixture (Figure 3b) also showed the distinction among five gene pools, including lowland North, lowland South, upland East, upland North, and upland West. Based on the principal component (PC) analysis and admixture, there was a population structure in this association panel ( Figure 3a). From PC1 and PC2, upland and lowland ecotypes separated from each group. Moreover, within the lowland ecotype, the latitude of the lowland can be distinguished to lowland north and south groups. The separation of the upland ecotype can also be noticed in PC1 and PC3 that the upland north was grouped apart from the other upland. Besides ecotypes, ploidy levels can be differentiated into groups. Both lowland north and south, and upland north genotypes were tetraploid (4X), while the upland east and upland west genotypes were octaploid (8X). In total, PC1, PC2, and PC3 explained the variance due to the population of 49.35%. The admixture (Figure 3b) also showed the distinction among five gene pools, including lowland North, lowland South, upland East, upland North, and upland West.

Genome-Wide Association and Candidate BLS Resistance Genes
We conducted GWAS with five subsets of genotypes, including 478, 4X, 8X, lowland, and upland. Most of the single phenotypes did not show significant Manhattan peaks or theoretically expected QQ-plots. To improve the association analysis, multiple combinations of BLS severity approaches were implemented. The full set of 478 genotypes conducted the GWAS with the combination of three traits of DTVI, DTIA, and DSIA, which showed 14 significant peaks over the FDR threshold ( Figure 4a and Table 2).

Genome-Wide Association and Candidate BLS Resistance Genes
We conducted GWAS with five subsets of genotypes, including 478, 4X, 8X, lowland, and upland. Most of the single phenotypes did not show significant Manhattan peaks or theoretically expected QQ-plots. To improve the association analysis, multiple combinations of BLS severity approaches were implemented. The full set of 478 genotypes conducted the GWAS with the combination of three traits of DTVI, DTIA, and DSIA, which showed 14 significant peaks over the FDR threshold ( Figure 4a and Table 2).
In a tetraploid group, the two-approached combination of DSVI and MNY provided three significant peaks with the expected QQ-plot ( Figure 4b). Interestingly, one peak of SNP Chr07N_61547213 on chromosome 7N from tetraploid group overlapped with the same peak in GWAS from the 478 group. In the octaploid group, the two-approached combination of DTVI-DSIA showed one significant peak on chromosome 9N (Figure 4c). In the lowland group, the best GWAS was from the single trait mean of two locations yielding    In this study, therefore, we focused on five selected GWAS analyses, including 478-DTVI-DTIA-DSIA, 4X-DSVI-MNY, 8X-DTVI-DSIA, lowland-mean TL, and upland-DTIA-MTL ( Figure 4). In total, there were 27 significant peaks (one overlapped between 478 and the lowland group) across 12 chromosomes: 1K, 2K, 2N, 3K, 3N, 4N, 5K, 5N, 6N, 7K, 7N, and 9N. Accumulatively, when we considered phenotypic variance explained (PVE) by all 27 significant SNP markers for each trait, PVE of DSIA was the highest at 26.52% and PVE of DTIA was the lowest at 3.28% (Table 3).  (Table S6). In 4X-DSVI-MNY, there were three SNPs on chromosome 3K, 3N, and 7N, accumulatively explaining 9.78% of DSVI and 11.31% of MNY. In 8X-DTVI-DSIA, only one significant marker showed on chromosome 9N, explaining 0.64% of DTVI and 11.31% of MNY. In lowland-mean TL, the eight significant markers were on chromosomes 1K, 5N, 6N, 7N, and 9N, explaining 45% of mean TL. In up-DTIA-MTL, the significant markers were on chromosome 1K and 7K explaining 5.03% of DTIA and 7.77% of MTL. The PVEs of each SNP in each subgroup and each trait were reported in Table S7.
In addition to the effect size of each SNP on each trait via PVE, the direction and inheritance of genotypic effects were shown in the format of the allele combination effect of each SNP on phenotypes ( Figure 5). Of 27 SNPs, there were 11 SNPs that showed dominant effects, such as Chr05K.5450256 (Figure 5a), Chr02N.79645558, Chr07k.50656672, to name a few. Five SNPs showed overdominance, including Chr07N.61547213 (Figure 5b) (Table S6). In 4X-DSVI-MNY, there were three SNPs on chromosome 3K, 3N, and 7N, accumulatively explaining 9.78% of DSVI and 11.31% of MNY. In 8X-DTVI-DSIA, only one significant marker showed on chromosome 9N, explaining 0.64% of DTVI and 11.31% of MNY. In lowland-mean TL, the eight significant markers were on chromosomes 1K, 5N, 6N, 7N, and 9N, explaining 45% of mean TL. In up-DTIA-MTL, the significant markers were on chromosome 1K and 7K explaining 5.03% of DTIA and 7.77% of MTL. The PVEs of each SNP in each subgroup and each trait were reported in Table S7.
In addition to the effect size of each SNP on each trait via PVE, the direction and inheritance of genotypic effects were shown in the format of the allele combination effect of each SNP on phenotypes ( Figure 5). Of 27 SNPs, there were 11 SNPs that showed dominant effects, such as Chr05K.5450256 (Figure 5a), Chr02N.79645558, Chr07k.50656672, to name a few. Five SNPs showed overdominance, including Chr07N.61547213 (Figure 5b  The effects of three different allele combinations were tested at alpha = 0.05. "ns" means non-significant, "*" means significant at alpha 0.05, "**" means significant at alpha 0.01, "***" means significant at alpha 0.001, and "****" means significant at alpha 0.0001. When we compared these significant SNP markers from GWAS to the potential candidate genes for resistance to BLS in rice, Chr07N_48546809, on chromosome 7N from 478-DTVI-DTIA-DSIA linking with Pavir.7NG240300 overlapped with the BLAST hit from LOC_Os04g39430.1 on rice chromosome 4 from AE11005627 marker (Tables S4 and S5). Moreover, Chr07_54502620, on chromosome 7K from up-DTIA-MTL linking with Pavir.7KG238800 overlapped with the BLAST hit from LOC_Os04g43730 on rice chromosome 4 from ad04009558 marker. showed an additive effect. The effects of three different allele combinations were tested at alpha = 0.05. "ns" means non-significant, "*" means significant at alpha 0.05, "**" means significant at alpha 0.01, "***" means significant at alpha 0.001, and "****" means significant at alpha 0.0001.

Discussion
When we compared these significant SNP markers from GWAS to the potential candidate genes for resistance to BLS in rice, Chr07N_48546809, on chromosome 7N from 478-DTVI-DTIA-DSIA linking with Pavir.7NG240300 overlapped with the BLAST hit from LOC_Os04g39430.1 on rice chromosome 4 from AE11005627 marker (Tables S4 and  S5). Moreover, Chr07_54502620, on chromosome 7K from up-DTIA-MTL linking with Pavir.7KG238800 overlapped with the BLAST hit from LOC_Os04g43730 on rice chromosome 4 from ad04009558 marker.

The Most BLS Resistant Population
Since the various phenotyping approaches yielded various resistant populations, it was challenging to determine the most resistant candidates for further breeding. Based on the field evaluation, BLUPs from two locations were zero due to zero broad-sense heritability. The different trend for BLS resistance between two locations was confirmed by the significant rank test. This re-ranking incidence suggested that the resistance was under effects between two locations [32,34]. If the mean of field scores were used, the upland 'Shelter' should be considered one of the most resistant populations across two locations with a mean of 0.44. The high resistance of 'Shelter' from field evaluation explained the reason in the prior study that it cannot be improved for resistance to BLS via recurrent phenotypic selection in seedlings. However, the one of most susceptible 'SW787' had a mean score of 2.33 out of 5. This suggested the low natural inoculation in 2017. When comparing each location, 'SW803' showed the lowest severity in NY but had the highest score in PA at 2.11. Therefore, this location variation needs to be reduced.
To lessen the GxE effect on the BLS, artificial inoculation was conducted in the laboratory. To handle the variation among sets and unbalanced design [35][36][37][38], BLUPs were used by fitting the effect from the set as a fixed effect. With the high correlation between vision and image analysis, the latter approach was recommended for the future BLS evaluation due to repeatability and reanalysis [39]. For the sake of comparison, the mean of percentage lesion from image analysis was considered. In the detached leaf approach, the lowest severity population from DTIA was ECS.6 at 8%, but the population appeared very susceptible to DSIA at 94%. These supported the low correlation among all three disease evaluation approaches. The low correlation between laboratory and field evaluations could be explained by the different disease pressure between the two conditions. The natural inoculation in the field varied by the weather condition of each year, whereas the laboratory inoculation was conducted with a high inoculum concentration to maximize the disease pressure that each genotype can respond to. Another counterintuition was the low correlation between detached leaf assay and leaf disk assay. This can be the result of the different ratios of leaf area to inoculum. The detached leaf assay used a 5-cm long whole leaf applied by spraying inoculum, whereas the leaf disk assay used only an 8-mm bored leaf section and a 2-microliter inoculum droplet on a single location. The difference between spraying and droplet was also shown in the assessment of early blight (Alternaria solani) in tomatoes [40].
The Qst-Fst of all phenotypes showed that none of them has gone under positive selection across five ecotypes. This neutrality was not expected as most of the resistance genes were hypothesized to be under natural selection [41][42][43][44][45]. However, neutral evolution of resistance genes like NBS-LRRs can be undergone in a relaxed selective constraint. For example, when pathogens are absent or at low frequency, the resistance genes may experience much weaker or no cost of resistance [46]. Such a low disease pressure can be noticed by means of field evaluation that the highest severity was only 2.33 from 5 (the most severe).

Dissecting BLS Resistance
There were only two significant SNP markers that overlap with the BLAST hits of candidate resistance genes from rice Chr07N_48546809 from the 478 GWAS and Chr07K_54502620 from upland GWAS ( Table 2 and Table S6). Chr07N_48546809 linked to Pavir.7NG240300 encoding for cytochromes P450, which involved diverse oxidation reactions and triterpene synthesis [47]. The most well-known triterpene against the fungal disease was Avenacins in oats (Avena spp.). Chr07K_54502620 linked to Pavir.7KG238800 encoding for EGF_CA, which has been proved to be the main structure of ZmWAK conferring resistance to maize head smut by controlling the galacturonan-binding for cell wall mediation [48]. In addition to genes overlapping with the BLAST hits of the resistance rice genes, Chr07N_61547213 showed a significant peak in both 478 and 4X GWAS. The SNPs linked with Pavir.7NG322400 encoding for PPR repeat, which was one of the potential resistance genes in rice. These PPR repeat families were commonly known for disease resistance genes [49]. Pavir.2NG397900 encoding GDP-fucose protein O-fucosyltransferase relating to cell wall development and disease resistance [50].
Although the rest of the significant peaks are unique within each group of GWAS, we considered the biological functions for the resistance over the fixation of the alleles in each population for a better resistance mechanism. In the full set of 478, on chromosome 2N, Pavir.2KG226200 was predicted to function as a zinc-binding RING finger. The zinc finger domains had the broad-spectrum nature of rice blast resistance gene Pi54 classed as NBS-LRR, which played an important role in effector-triggered immunity (ETI) [51]. Moreover, this gene also encoded glycosyltransferase, which was responsible for cell wall synthesis and modification, leading to disease resistance [52]. Pavir.2NG424700 encoded for pyruvate dehydrogenase, which was a key protein in a tricarboxylic acid cycle. This suggested the high energy-intensive resistance response as happened in wheat leaf rust (Puccinia triticina Eriks) [53]. Pavir.2NG500000 encoded for casein kinase II. The protein itself did not confer the resistance, but it was often found as many motifs in Mlo [54], which was identified to provide broad resistance to powdery mildew in barley [55]. On chromosome 3K, Pavir.3KG329700 encoded DNA polymerase III, mainly known for DNA replications and repairing DNA damaged by hydrogen peroxide [56]. The DNA repair process was a part of the plant immune response [57], especially in the infection of B. oryzae that hydrogen peroxide was shown to accumulate in a leaf [23]. Pavir.3KG358400 encoded for serine aminopeptidase that conferred no resistance but played a major role in necrosis [58]. On chromosome 4N, Chr04N_47494808 was more closely linked to Pavir.4NG267200, which was predicted to have a hydrolase activity. The further gene, Pavir.4NG267500, was encoded for Villin, which was tissue-specific actin modifying protein responsible for anti-apoptotic activity suppressing necrosis [59,60]. On chromosome 5N, Pavir.5KG028600 encoded for tetratricopeptide repeat, which was one of five domains of the SGT1 resistance gene in Arabidopsis [61]. Pavir.5KG482000 encoded for ATP-binding cassette (ABC) transporter. There are many subgroups of ABC. In Arabidopsis, the mutant lacking the ABC transporter of penetration3/pleotropic drug resistance8 (PEN3/PDR8) conferred nonhost susceptibility to Blumeria graminis, suggesting that the ABC transporter was involved in exporting secreted toxins from the fungus and fungal suppressor from the host [62]. Moreover, in wheat, the ABC suggested the exportation of mycotoxin, conferring resistance to Fusarium head blight [63]. On chromosome 7K, Pavir.7KG196600 encoded for interleukin-1 receptor-associated kinase, which was the pathogen recognition resulting in plant response [64]. Pavir.7KG255000 encoded for oligopeptide transporter protein, which was also shown in rice QTL resistance to rice blast (Manaporthe grisea) [65]. Lastly, on chromosome 9N, Pavir.9NG844200 encoded for Mitochondrial Fe2+ transporter MMT1. Despite no direct report of resistance, mitochondria were proven to control redox homeostasis under stress [66].
In GWAS from 4X-DSVI-MNY, on chromosome 3K, Pavir.3KG095200 encoded for both peroxisomal membrane protein and ABC transporter. Due to the high involvement of reactive oxygen species (ROS) in the infection of B. oryzae [67], peroxidase played an important role in the pathogen-associated molecular pattern-triggered immunity (PTI) [68]. Additionally, Pavir.3NG101800 on chromosome 3N encoded ABC transporter.
The GWAS from 8X-DTVI-DSIA yielded only one significant SNP on chromosome 9N. Pavir.9NG061100 encoded for DNA topoisomerase and, more importantly, dirigent. Dirigent protein conferred resistance by the lignin and lignan synthesis [69]. The lignin synthesis was expected to relate to BLS resistance by strengthening cell walls [70].
In GWAS from lowland mean of two locations, on chromosome 1K, the further linked Pavir.1KG225100 encoded for MORC family-ATPases. This protein was the key component compromising the recognition of Turnip Crinkle Virus (CRT1) in Arabidopsis that was required for PTI, basal resistance, no-host resistance, and systemic acquired resistance [71]. On chromosome 5N, Pavir.5NG192000 encoded for COP9 signalosome complex. Interestingly, this signalosome complex interacts with SGT1, indicating the disease resistance [61]. Pavir.5NG248900 encoded basic leucine zipper (bZIP) transcription factor, which was a key modulator in hypersensitive response and salicylic acid [72]. Pavir.5NG476900 encoded for abscisic-acid-induced TB2/DP1, which was responsible for membrane turnover and reduced unnecessary secretion against Southern corn leaf blight by Cochliobolus heterostrophus (Drechs.) Drechs. in maize [73]. On chromosome 6N, Pavir.6NG083100 encoded homeobox-leucine zipper protein playing an important role in responses to abscisic acid under stress [74], transcriptional regulation for disease resistance in rice [75], and programmed cell death [76]. On chromosome 7N, Pavir.7NG091000 encoded sucrose synthase. The hexose sensing was important to activate defense-related genes as well as repression of photosynthetic genes [77]. On chromosome 9N, Pavir.9NG173600 encoded Phytochromeinteracting factor 4. The photoreceptors in plants played an essential role in hydrogen peroxide and ROS in the cell resulting in disease response [78]. Lastly, Pavir.9NG499500 encoded for alpha/beta hydrolase, which was a diversely biochemical protein family involving the activation of hydrogen peroxide [79].
In the last GWAS from upland DTIA-MTL, Chr01K_61692577 was linked to two interesting genes. First, Pavir.1KG372200 encoded Exo70 exocyst complex subunit, which was required for recognizing the avirulent effector AVR-Pii from rice blast resulting in ETI [80]. The other linked gene was Pavir.1KG372200 encoding for Nicotinamidase, which showed resistance to tomato mosaic virus [81].
In conclusion, the resistance to BLS from the different phenotyping approaches with low correlations suggested that resistance was sensitive to phenotyping methods in different conditions. Although it is difficult to standardize the phenotyping approach that can take into account all conditions, 'SW788', 'SW806', 'SW802', 'SW793', 'SW781', 'SW797', 'SW798', 'SW803', 'SW795', 'SW805 should be considered as candidates for resistance. To increase the statistical power in SNP detection, the five subgroups of populations were conducted with multi-trait GWAS. The full population of 478 genotypes showed the most SNP markers that can be considered as potential candidate genes related to resistance to BLS with reasonable biological functions. Additionally, the SNP markers from the other subgroups provided a broader resistance mechanism. Nevertheless, one needs to be prudent regarding the utilization of these markers. The dependency of multi-trait GWAS to yield significant SNP peaks brought the challenge for genomic-assisted breeding for resistance to BLS.
When we considered biological functions from the GWAS regardless of subgroups and trait combinations, the defense mechanisms against BLS were complicated and multifaceted, probably controlled by multiple small-effect alleles. The defense mechanism started from basal defense via cell-wall mediated proteins and PTI via peroxidase. Then, ETI depending on NBS-LRR and PPR repeat took place. The responses to necrosis were suppressed by homeostasis and anti-apoptotic proteins. To our knowledge, this was the first research to dissect the resistance to BLS in switchgrass. With the current power of detection within the population, multi-trait GWAS can provide functional insight into the resistance within multiple subgroups of genotypes.

Switchgrass in the Northern Association Panel
The switchgrass used in this study has been developed for accelerating breeding progress, especially for bioenergy traits at northern latitudes [9,33]. The association panel consisted of 478 genotypes from 66 populations representing the mostly upland northern population and some southern lowland population (Table S8). 'Population' in this study means a seed source with a specific origin and, in some cases, breeding history. Six to ten seeds from each population were planted into individual plants called 'genotypes.' The association panel was initially planted in Ithaca, NY, in 2008, and then vegetatively cloned and planted in a randomized, complete block design with three replicates in 0.9 m spaced planting in Ithaca, NY, USA and Phillipsburg, PA, USA, in 2016. No fungicide has been applied to the field.

Disease Evaluation and Phenotype Processing
Since the recurrent phenotypic selection for resistance to BLS cannot improve the trait from screening in the seedling stage [26], the phenotypic evaluations in this study were explored in a mature stage under different conditions. Bipolaris leaf spot has been evaluated in the fields to determine the resistance in natural incidence under different environments and in the laboratory to minimize environmental effects on the resistance.
Field evaluation was conducted in the same week in both locations in August 2017. Each plant was visually evaluated by a single person to reduce the variance among persons, the severity ranging from 0 to 5 as 0 = 0% BLS, 1 = 1-10% BLS, 2 = 11-25% BLS, 3 = 26-50% 134 BLS, 4 = 51-75% BLS, and 5 = 76-100% of the leaf area of the whole plant with BLS. In addition to field evaluation, to minimize environmental effects, an artificial inoculation in the laboratory was conducted. Disease severity was evaluated by detached leaf and leaf disk assay by taking leaf samples from one replicate of genotypes in Ithaca, NY. In the detached leaf assay, three healthy leaves were randomly selected from 30 cm below the top of the shoot to control the same stage of leaf development. The leaves were cut into five centimeters, washed with deionized water, placed in a pre-wet petri dish bedded with filter paper, and kept cold in a cooler in the field, then refrigerated overnight.
The inoculum was prepared from a subculture of a single conidium isolated from a 'Carthage' switchgrass leaf in the warm season biofuels field experiment in Ithaca, NY (Waxman, 2011). The day after collecting the detached leaf, the three-week-old plate was flooded, filtered via gauze, and the concentration adjusted to 10 5 conidia·mL −1 by hemocytometer with two drops of Tween-20 per 100 mL. An airbrush pressuring at ten psi was used to spray inoculum on each dry-surfaced adaxial side of the leaf on each plate. After letting the droplets of inoculum dry and firmly attached to the leaf surface, the plates were then sealed with paraffin tape to maintain high moisture conditions. Each plate contained three detached leaves from one genotype. The plates were placed on the table randomly and kept at room temperature with 12-h light for seven days showing necrotic lesions ( Figure 6). In each set of the experiment, 85 to 200 genotypes were sampled with some overlapping genotypes in one to three duplicated plates and additional noninoculated plates as controls. The leaf samplings were done weekly from May to June 2017 and controlled the similar stage of a leaf by sampling from the same height from the top. A total of seven sets of the experiment eventually covered all 478 genotypes. The disease evaluation was conducted by vision and image analysis software. Both evaluations were based on a percentage of leaves covered with lesions. The evaluation by vision was from 0 to 100% with a 10% increment. Then, each leaf was dried and scanned with a flatbed scanner (Canon CanoScan LiDE 700F manufactured by Canon Inc, Tokyo, Japan) at a resolution of 1200 dpi, and images were saved in .tiff format. The best linear unbiased predictions (BLUPs) were applied to extract random effects of the genotypes fitting the experiment factor as a fixed effect [35][36][37][38]. In the leaf disk evaluation, all 478 genotypes were collected on the same day in August 2017 by cutting the same age leaves from each genotype and keeping them refrigerated overnight. On the next day, the leaves were washed with deionized water, surface dried, bored by keeping the midrib in the middle into an 8-mm disk, and placed on a water agar plate. There were 28 leaf disks in each plate, including four check disks (control and inoculated disks of resistant 'SW43_09' and susceptible 'SW122_02') and three repli- In the leaf disk evaluation, all 478 genotypes were collected on the same day in August 2017 by cutting the same age leaves from each genotype and keeping them refrigerated overnight. On the next day, the leaves were washed with deionized water, surface dried, bored by keeping the midrib in the middle into an 8-mm disk, and placed on a water agar plate. There were 28 leaf disks in each plate, including four check disks (control and inoculated disks of resistant 'SW43_09' and susceptible 'SW122_02') and three replicated disks of eight genotypes (Figure 6d). The inoculum was prepared the same way described above. For each leaf, 2 µL of inoculum was dropped in the middle of the right side of the leaf. After letting the droplet dry, all 64 plates were sealed with paraffin tape to keep high moisture, randomly placed on a table, and kept under 12-h light at room temperature. The disease evaluation was conducted at seven dpi both by vision and image analysis software. The percent of leaf disk covered with lesions and necrotic tissue was assessed by a vision from 0% to 100% with 10% increments. Moreover, on the same day, a photo was taken of each plate with a digital DSLR Canon Rebel T6 DSLR camera with a resolution of 2644 × 3084 in .JPG.
Images from both detached leaf and leaf disk assays were analyzed using ImageJ macro (the settings available upon request). In brief, each detached leaf or leaf disk was measured for total leaf area and total necrosis lesion area (yellow and black area). The percentage of the area covered with lesions was computed by dividing the total necrotic lesion area by the total leaf area.
Package 'LME4' [82] in R was used to calculate BLUPs from various disease evaluations. For the field evaluation, BLUPs were computed based on each location separately and on the two locations combined using genotypes nested in populations and replicates as random effects in each location model and populations as a fixed effect. BLUPs for the two locations combined were fitted genotype nested in population, replicates, locations, the interaction between genotypes and location, as random effects, and populations as a fixed effect. In addition to generating BLUPs from field evaluation, the highest score, which was the most severe symptoms observed from each location and two locations over replicates, were used due to successful QTLs of resistance to foliar symptoms caused by potato virus Y in autotetraploid potato (Solanum tuberosum L.) [83]. The highest scores were used on the assumption that in the field, there was a chance that some replicates may expose or avoid the pathogens differently. The maximum scores were considered the worst response of each genotype. For disease evaluations under laboratory conditions, since all phenotypes were expressed as percentages to leaf area, log transformation was performed before computing BLUPs. For the detached leaf assay, BLUPs were computed with experiment, genotype nested in populations, the interaction between genotypes and experiment, plates within the experiment, and leaf replicates of genotypes as random effects and populations as a fixed effect. For the leaf disk assay, BLUPs were calculated with plates as a fixed effect and genotype, replicates within plates, and genotype within plates as random effects and populations as a fixed effect. Broad-sense heritability was estimated from variances in each model, and standard error was computed by the bootstrap approach 1000 times (Supplement file). Moreover, phenotypic correlations were computed between resistance to Bipolaris leaf spot to the other 20 morphological and biomass quality traits from a previous study [33], such as plant height, anthesis date, acid detergent lignin, minerals, ethanol·g −1 dry weight, etc.
Therefore, the resistance to BLS was evaluated via three approaches-detached leaf, leaf disk, and field evaluation-providing phenotyped traits including BLUPs of detached leaf percent lesion via vision and image analysis (DTVI and DTIA), BLUPs of leaf disk percent lesion via vision and image analysis (DSVI and DSIA), the highest score from two locations, NY and PA (MTL, MNY, and MPA).

Genotyping, Linkage Disequilibrium Analysis and Population Structure
The exome-capture sequencing was conducted as previously described [84]. In short, the DNA of each genotype was processed via exome-capture using the Roche-Nimblegen exome-capture probe set [85] and DNA sequencing. The sequences were aligned to the P. virgatum genome assembly v.4.1 (P. virgatum v.4.1 [86]) for SNP discovery [87]. The reads were sorted with PicardTools version 2.1.1 [88] via functions of SortSam and MarkDuplicates. Samtools version 1.3.1 [89] was used to pile up files with base alignment quality disabled and map quality adjustment disabled. In this calling, SNPs were required to be biallelic, sequenced in all samples, be monomorphic in at least two samples, and filtered with a minimum read mapping quality score of 30 and more than 5X coverage in more than 95% of the samples. At each SNP, the genotype dosages range from zero to two copies of minor alleles and can be nonintegers by using the EM algorithm [90]. Naturally, the switchgrass association panel includes allopolyploid switchgrass: tetraploids (4×), octoploids (8x), and hexaploids (6x). It was challenging to perform GWAS with polyploid models. In this study, we modeled them under the assumption of disomic inheritance, similar to the study of GWAS of flowering time in the previous study [29]. This was because tetraploid switchgrass was confirmed to have disomic inheritance [91]. Although octaploid switchgrass had four copies of each homologous chromosome, which made it difficult to precisely determine heterozygous genotypes, the disomic segregation was used for the model with the caution of the increasing standard error of estimates in GWAS. Moreover, linkage disequilibrium (LD) was computed via the '-r2' tag in PLINK as r 2 between all SNPs within 1 MB and recorded all values [92]. Principle component analysis was used to evaluate the population structure via PCA in the FactoMineR package [93]. Additionally, ADMIXTURE was used to evaluate population structure with five-fold cross-validation [94].
We were also interested in testing whether the variations in phenotypes were from the high differentiation across five ecotypes (lowland north, lowland south, upland east, upland north, upland west) based on Qst-Fst analysis. Estimated Qst-Fsts were computed by R package "QstFstComp" by bootstrapping 1000 times by sampling random genotypes [95].

Candidate Genes Based on Resistance in Rice
Since B. oryzae is one of the major pathogens in rice (Oryza sativa), major resistant QTLs have been identified in recombinant inbred lines (RILs) and doubled haploids (DH) using restriction fragment length polymorphism (RFLP), simple sequence repeat (SSR) markers, and sequence-tagged site (STS) [30,31,94]. Three major QTLs in Chromosomes 1, 4, and 11 were then studied further in near-isogenic lines (NILs) and 13 SNP markers linked with the resistance [30]. The NILs are BC3F5 between indica 'Tadukan' (resistance) and temperate japonica 'Koshihikari' (susceptible). According to the previous GWAS in rice, linkage disequilibrium (LD) was~100 kb in indica and~200 kb temperate in japonica [96]. Thus, in this study, the candidate resistance genes were screened from 200 kb around the 13 SNPs from Oryza sativa genome assembly v.7.0 (Oryza sativa v.7.0 [86]). The candidate genes from rice are listed in Table S6. The sequences of these candidate genes were used to identify potential homologs in switchgrass using the BLASTP tool on Phytozome and kept the hits with E-value < 10 −10 [86].

Genome-Wide Association Studies
Genome-wide Efficient Mixed Model Association (GEMMA) [97] was used to implement a multivariate linear mixed model for GWAS of single-and multiple-phenotyping approaches for resistance to BLS to analyze each phenotype and combined phenotypes for improving statistical power [98,99]. Kinship was included as a random effect. Moreover, the first three principal components (PCs) were included as fixed covariates. The GWAS was conducted in five groups, including all 478 genotypes, tetraploid (4×), octaploid (8X), lowland, and upland genotypes, to determine if there were any SNPs linked to specific switchgrass population. Minor allele frequency (MAF) at 0.05 was used to filter the minor alleles. In total, there are 135,791 SNP markers used in association analysis. To correct multiple testing, the false discovery rate (FDR) was calculated as where P is the p-value tested at 0.1; k is the number of traits in the combination; A is the number of SNP that were significant at the p-value tested; T is the total number of SNP tested. The significant markers then were used to determine candidate genes linked to them within the range of LD and search in JBrowse in P. virgatum v.4.1 in Phytozome (e.g., if the position of the marker was 4019500 on chromosome 9K and LD was 20 kb, "Chr09K:4019500..4039500" was used). In case the candidate gene was not fully aligned within that LD range, we considered the gene as a candidate if 80% of its length was within the LD range. Since some of the SNP positions could not be aligned to the P. virgatum genome assembly v.4.1, they were excluded from determining the candidate genes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11101362/s1, Figure S1: Bar chart of mean severity score from 0 to 5 from two location (top), mean from NY (middle) and mean from PA with standard error from each 66 population; Figure S2: Bar chart of mean severity from the leaf detachment assay via image analysis (DTIA) (top), from the leaf disk assay via image analysis (DSIA) (middle) and mean field severity from two locations with standard error from each of 66 populations; Figure S3: Correlation plot among BLUPs of severity from the leaf detachment via vision (DTVI), via image analysis (DTIA), from leaf disk assay via vision (DSVI), via image analysis (DSIA), the highest score between two locations (MTL), mean from two locations (mean TL), highest score in NY (MNY), highest score in PA (MPA) and BLUPs of other agronomic and biomass quality traits; Figure S4: The bootstrapped distributions of Qst-Fst for each phenotypes; Table S1: Phenotype summary of the severity of Bipolaris leaf spot from mean of Bipolaris lesion percentage from detached leaf assay by vision (DTVI), by image analysis (DTIA), the severity from the leaf disk assay by vision (DSVI), by image analysis (DSIA), mean of two locations, highest scores between two locations (MTL), mean in NY, highest score in NY (MNY), mean in PA, and highest score in PA (MPA) based on population (rank); Table S2: Phenotype summary of DTI, DTIA, DSVI, DSIA, mean TL, NY and PA, and MTL, MNY, and MPA based on each genotype (rank) ordered based on DTIA; Table S3: Summary of P-value from ANOVA tests of different groups of ecotypes in each trait: detached leaf assay by vision (DTVI), by image analysis (DTIA), the severity from the leaf disk assay by vision (DSVI), by image analysis (DSIA), mean of two locations, highest scores between two locations (MTL), mean in NY, highest score in NY (MNY), mean in PA, and highest score in PA (MPA) and their ranks; Table  S4: Candidate genes and their annotated functions for resistance to BLS in rice (Oryza sativa) from QTLs from Sato et al. (2015) within LD of 200 kb; Table S5: BLAST results from candidate resistance genes from rice to switchgrass; Table S6: Phenotypic variance explained (adjusted PVE) by set of significant SNP markers in each subgroup for traits used in each subgroup; Table S7: Phenotype variance explained (PVE) by each significant SNP markers in each subgroup for each trait; Table S8: Switchgrass association panel sample information.