Germplasm Screening Using DNA Markers and Genome-Wide Association Study for the Identification of Powdery Mildew Resistance Loci in Tomato

Powdery mildew (PM), caused by Oidium spp. in tomato, is a global concern that leads to diminished yield. We aimed to evaluate previously reported DNA markers linked to powdery mildew resistance (PMR) and identify novel quantitative trait loci (QTLs) for PMR through a genome-wide association study in tomato. Sequencing analysis of the internal transcribed spacer (ITS) of a PM strain (PNU_PM) isolated from Miryang, Gyeongnam, led to its identification as Oidium neolycopersici. Thereafter, a PM bioassay was conducted for a total of 295 tomato accessions, among which 24 accessions (4 S. lycopersicum accessions and 20 accessions of seven wild species) showed high levels of resistance to PNU_PM. Subsequently, we genotyped 11 markers previously linked to PMR in 56 accessions. PMR-specific banding patterns were detected in 15/22 PMR accessions, while no such bands were observed in the powdery mildew-susceptible accessions. The genome-wide association study was performed using TASSEL and GAPIT, based on the phenotypic data of 290 accessions and 11,912 single nucleotide polymorphisms (SNPs) obtained from the Axiom® Tomato SNP Chip Array. Nine significant SNPs in chromosomes 1, 4, 6, 8, and 12, were selected and five novel QTL regions distinct from previously known PMR-QTL regions were identified. Of these QTL regions, three putative candidate genes for PMR were selected from chromosomes 4 and 8, including two nucleotide binding site-leucine rich repeat class genes and a receptor-like kinase gene, all of which have been identified previously as causative genes for PMR in several crop species. The SNPs discovered in these genes provide useful information for understanding the molecular basis of PMR and developing DNA markers for marker-assisted selection of PMR in tomato.


Introduction
Powdery mildew (PM) in tomatoes (Solanum lycopersicum L., 2n = 2x = 24) is a common concern worldwide, especially in areas with cool and dry climates. It reduces production by more than 50%, depending on the growth stage of the infected plants and the surrounding environment. The major pathogens involved in PM are Oidium neolycopersici, Oidium lycopersici, and Leveillula taurica, which belong to Ascomycota [1]. Among them, O. neolycopersici is the main cause of PM in tomato plant cultivation. When infected, a small circular, opaque white powdery sign appears on the upper surface of the lower leaf of the plant, which then spreads rapidly to cover the entire leaf [2]. When the symptoms become severe, the leaves fall off prematurely, leading to yellowing and necrosis, which ultimately interferes with photosynthesis, thereby lowering the growth rate and fruit quality of crops [3].
For economical, sustainable, and eco-friendly cultivation, there is a need for use of PMresistant varieties. However, commercial varieties known to be PM-resistant have hardly been developed, in addition to which, only some related wild species showing potential yield, fruit quality, and abiotic stress [21]. To the best of our knowledge, however, no GWAS for PMR in tomato has been reported to date.
In order to develop PM-resistant tomato varieties, it is necessary to continuously discover breeding materials with high PMR. It is also necessary to identify PMR genes and molecular markers that can be selected by application of various analysis techniques that use not only bi-parental or multiparent segregation populations, but also natural populations. Therefore, this study aimed to (1) screen PMR from a population of diverse tomato genetic resources, including wild species, (2) search for PMR loci based on previously reported PMR-associated markers, and (3) identify novel PMR-associated QTL(s) and candidate genes through SNP genotyping and GWAS, based on the Axiom ® Tomato Genotyping Array.

Pathogen Identification
We sequenced the internal transcribed spacer (ITS) region of the PM pathogen (named PNU-PM) found on the leaves of the tomato plants growing in the greenhouse of Pusan National University's Miryang Campus, Korea. The length of the ITS sequence of PNU-PM amplified using PMITS1 and PMITS2 primers was found to be about 700 bp (Figure 1a). The ITS sequences for a total of 29 fungal species including PNU-PM, L. taurica, and pathogens used in Kiss et al. [24] were aligned using Clustal X, to create a phylogenetic tree.
The ITS sequence of PNU-PM (NCBI accession number: MW082786.1) matched 100% with that of the O. neolycopersici France isolate and 99.8% with that of the Netherland isolate and was most closely tied with a bootstrap value of 86% on the phylogenetic tree ( Figure 1b). The ITS sequences of O. lycopersici and L. taurica, the other PM-causing pathogens in tomatoes, showed a 79.5% and 72.6% match, respectively, with that of PNU-PM, but were distinct from that of O. neolycopersici (Figure 1b). In addition, BLASTn search of the ITS sequence of PNU-PM showed a 100% match (query cover = 100% and E-value = 0) with the ITS sequence of 42 O. neolycopersici isolates.
A total of 48 accessions were moderately resistant (MR, PDI = 10-30%), including LA2172 (S. peruvianum), which carries the resistance gene Ol-4. The remaining 223 accessions were susceptible (S, PDI > 30%), among which 204 lines were S. lycopersicum, including susceptibility controls, PNU-PMS (S. lycopersicum) and Moneymaker as well as one S. habrochaites accession (PI247087) carrying the resistance gene Ol-5 (Table S1).  [24]; (b) A phylogenetic tree of powdery mildew pathogens constructed using the neighbor-joining method based on the ITS sequences of ribosomal DNA. The duplicate PCR amplicons were loaded on an agarose gel, and the specificity of the ITS region of 'PNU-PM' was confirmed based on the size (700 bp) of the PCR amplicons. Lane 'M', 100 bp DNA size marker. The ITS of the 'PNU-PM' powdery mildew pathogen (red in b) isolated in this study was found to be identical to that of O. neolycopersici. region of the Oidium neolycopersici isolate 'PNU-PM' using the PMITS1 and PMITS2 primers [24]; (b) A phylogenetic tree of powdery mildew pathogens constructed using the neighbor-joining method based on the ITS sequences of ribosomal DNA. The duplicate PCR amplicons were loaded on an agarose gel, and the specificity of the ITS region of 'PNU-PM' was confirmed based on the size (700 bp) of the PCR amplicons. Lane 'M', 100 bp DNA size marker. The ITS of the 'PNU-PM' powdery mildew pathogen (red in b) isolated in this study was found to be identical to that of O. neolycopersici. A total of 48 accessions were moderately resistant (MR, PDI = 10-30%), inclu LA2172 (S. peruvianum), which carries the resistance gene Ol-4. The remaining 223 a sions were susceptible (S, PDI > 30%), among which 204 lines were S. lycopersicum, in ing susceptibility controls, PNU-PMS (S. lycopersicum) and Moneymaker as well as o habrochaites accession (PI247087) carrying the resistance gene Ol-5 (Table S1).

Validation of the PMR-Linked Markers
Eleven previously reported PMR-linked DNA markers were genotyped in a to 56 accessions, including 22 accessions that showed resistance (R, PDI < 10%), 28 acces that showed sensitivity (S, PDI = 100%), and 6 controls (Table 1). PMR-specific m genotypes were observed in 15 of the resistant accessions but in none of the suscep accessions. In the case of the recessive resistance gene Slmlo1 at the locus ol-2, for wh gene-based marker has been developed, there was no accession carrying the marker otype for the Slmlo1 allele, but the marker genotype specific to Slmlo1.1, another alle the gene, was observed in the resistant accessions KNU16 and KNU17 (Table 1).

Validation of the PMR-Linked Markers
Eleven previously reported PMR-linked DNA markers were genotyped in a total of 56 accessions, including 22 accessions that showed resistance (R, PDI < 10%), 28 accessions that showed sensitivity (S, PDI = 100%), and 6 controls ( Table 2). PMR-specific marker genotypes were observed in 15 of the resistant accessions but in none of the susceptible accessions. In the case of the recessive resistance gene Slmlo1 at the locus ol-2, for which a gene-based marker has been developed, there was no accession carrying the marker genotype for the Slmlo1 allele, but the marker genotype specific to Slmlo1.1, another allele of the gene, was observed in the resistant accessions KNU16 and KNU17 (Table 2).

Quality Evaluation of the Raw SNP Genotyping Data
A total of 290 accessions were genotyped using the Axiom ® Tomato Genotyping Array, resulting in the raw data for 51,214 SNPs. Following that, for SNP quality evaluation in the wild species, Dish QC (DQC), call rate, and heterozygosity rate were analyzed for 182 accessions, excluding the 160 cultivated tomato accessions (Table S2). The average DQC of the S. lycopersicum accession was as high as 0.79-0.99, while that of the wild species accessions was very low (0.04-0.15). In the wild species, the DQCs of S. cheesmaniae, S. galapagense, and S. pimpinellifolium, which are known to exhibit high crossing compatibility with S. lycopersicum [25], were 0.15, 0.14, and 0.10, respectively, while those of other wild species were in a lower range of 0.03-0.04 (Table 3). The heterozygosity rate was as low as 5.43-7.44% in the cultivated species, but as high as 7.74-33.44% in the wild species. For call rate, the minimum was 88.64% (SAL1875, S. pimpinellifolium), and the maximum was 98.85% (LA3871, S. lycopersicum), thereby indicating that the overall SNP chip array genotyping was performed properly; however; a relatively low SNP call rate in the wild species implied higher number of missing call signals in this set as compared to that in the cultivated species (Table S2). A correlation analysis between the three factors of DQC, call rate, and heterozygosity rate revealed a significant correlation coefficient (0.672, p < 0.01) of DQC and call rate, which indicated that the higher the DQC, the higher the call rate. However, the correlation coefficients between DQC and heterozygosity rate as well as call rate and heterozygosity rate were -0.807 and -0.881, respectively, indicating a strong negative correlation between these factors, i.e., the heterozygosity rate tended to be high when the DQC or call rate was low.

Homology Analysis of the SNP-Flanking Probe Sequence among Species
To identify the cause of low DQC, call rate, and high heterozygosity rate in the wild species, sequence homology analysis was performed between the probe sequence of 51,214 SNP loci used in the Axiom ® Tomato Genotyping Array and the reference genome assembly of 5 cultivated and wild species obtained from National Center for Biotechnology Information (NCBI). Upon alignment of the probe sequences on the reference genome of each species, all of the probe sequences of the 51,214 SNP loci were mapped in S. lycopersicum, while only 36,811 SNP loci (71.88% of the total) were mapped in the wild species S. chilense, 38,335 SNP loci (74.85%) in S. habrochaites, 38,117 SNP loci (74.43%) in S. pennellii, 49,738 SNP loci (97.12%) in S. pimpinellifolium, and 42,550 SNP loci (83.08%) in S. arcanum; thus, as compared to S. lycopersicum, there was low homology with the SNPflanking probe sequence used in the Axiom ® Tomato genotyping array in the wild species.
Further, we assessed the relationship of the number of nucleotide sequence variations between the reference genome and probe sequences for each species with the ratio of SNPs genotyped as missing or heterozygous on the SNP chip array (Figure 3). In the cultivar accessions, there was no or only one nucleotide sequence variation, and the proportion of SNPs that were genotyped as missing or heterozygous were 3.97-3.98% and 6.63-6.64%, respectively. In the wild species, the ratio of missing/heterozygous genotyping for SNPs with a probe sequence that perfectly matched the reference genome sequence was 10.48% in S. chilense, 8.66% in S. habrochaites, 8.27% in S. pennellii, 14.01% in S. pimpinellifolium, and 9.80% in S. arcanum. Whereas, among SNPs whose probe sequences did not map to the reference genome (aligned 0 times), the ratio of missing/heterozygosity was 56.61% in S. chilense, 61.64% in S. habrochaites, 58.96% in S. pennellii, 28.00% in S. pimpinellifolium, and 50.78% in S. arcanum. As the overall degree of nucleotide sequence variations between the reference genome and probe sequences increased, there was a proportional increase in the ratio of SNPs genotyped as missing or heterozygous.
To identify the cause of low DQC, call rate, and high heterozygosity rate in the wild species, sequence homology analysis was performed between the probe sequence of 51,214 SNP loci used in the Axiom ® Tomato Genotyping Array and the reference genome assembly of 5 cultivated and wild species obtained from National Center for Biotechnology Information (NCBI). Upon alignment of the probe sequences on the reference genome of each species, all of the probe sequences of the 51,214 SNP loci were mapped in S. lycopersicum, while only 36,811 SNP loci (71.88% of the total) were mapped in the wild species S. chilense, 38,335 SNP loci (74.85%) in S. habrochaites, 38,117 SNP loci (74.43%) in S. pennellii, 49,738 SNP loci (97.12%) in S. pimpinellifolium, and 42,550 SNP loci (83.08%) in S. arcanum; thus, as compared to S. lycopersicum, there was low homology with the SNP-flanking probe sequence used in the Axiom ® Tomato genotyping array in the wild species.
Further, we assessed the relationship of the number of nucleotide sequence variations between the reference genome and probe sequences for each species with the ratio of SNPs genotyped as missing or heterozygous on the SNP chip array (Figure 3). In the cultivar accessions, there was no or only one nucleotide sequence variation, and the proportion of SNPs that were genotyped as missing or heterozygous were 3.97-3.98% and 6.63-6.64%, respectively. In the wild species, the ratio of missing/heterozygous genotyping for SNPs with a probe sequence that perfectly matched the reference genome sequence was 10.48% in S. chilense, 8.66% in S. habrochaites, 8.27% in S. pennellii, 14.01% in S. pimpinellifolium, and 9.80% in S. arcanum. Whereas, among SNPs whose probe sequences did not map to the reference genome (aligned 0 times), the ratio of missing/heterozygosity was 56.61% in S. chilense, 61.64% in S. habrochaites, 58.96% in S. pennellii, 28.00% in S. pimpinellifolium, and 50.78% in S. arcanum. As the overall degree of nucleotide sequence variations between the reference genome and probe sequences increased, there was a proportional increase in the ratio of SNPs genotyped as missing or heterozygous.

SNP Validation Based on Cleaved Amplified Polymorphic Sequence (CAPS) and Derived Cleaved Amplified Polymorphic Sequence (dCAPS) Markers
To verify the reliability of the SNP chip genotyping results of wild species, six SNPs were converted into CAPS or dCAPS markers and the genotypes of the highly resistant and highly susceptible accessions were reanalyzed (Table 4). Of the 22 resistant accessions, SNP call in the SNP chip array was 'missing' in 3 accessions for the SNP AX-95781451, 1 for AX-95767557, 4 for AX-95784118, 3 for AX-95809776, and 2 for AX-95814666, while the SNPs were genotyped as heterozygous in 6 accessions for AX-95781451, 3 for AX-95767557, 4 for AX-95784118, 18 for AX-95773152, 8 for AX-95809776, and 10 for AX-95814666.
The PCR-based marker genotypes are presented in parentheses when they did not match the genotyping results from the SNP chip array. 2 NA: Missing SNP call signal.
However, all CAPS or dCAPS markers converted from the six SNPs showed homozygous genotypes, and this trend was mainly observed in wild species accessions, except for the case of AX-95773152 in three cultivars, A1161, A1162, and A1164. In 28 susceptible lines mainly composed of cultivars, CAPS or dCAPS markers for all 6 SNPs showed a showed homozygous genotype, and the test results between the SNP chip array and the CAPS or dCAPS markers matched, except for in case of the genotype of AX-95781451 in one accession, A1345.

SNP Filtering for Population Structure and GWAS
Filtering SNPs under the conditions of minor allele frequency >5%, missing proportion <5%, and heterozygosity proportion <10% led to the selection of 28,374 (55.4%) SNPs. Further, to avoid low genotyping error, we filtered SNPs with a heterozygosity rate lower than 10% in the wild species, which led to a final selection of 11,912 (23.2%) SNPs (Table S3).

Population Structure Analysis
Population structure analysis of 290 accessions was performed based on the filtered 11,912 SNPs, using the STRUCTURE analysis (Figure 4a). The optimal number of subgroups (delta K) of the group was 6. Most of the lines belonging to cultivar S. lycopersicum were divided into Groups 1, 2, 3, and 4, with 16 PGRBC lines and BT0902 lines in Group 1, 27 PGRBC lines in Group 2, and 97 PGRBC lines and lycopersicoides introgression in Group 3. There were 21 lycopersicoides IL lines, 7 hirsutum IL lines, and 7 other lineages. Group 4 included 9 accessions of PGRBC, 1 accession of lycopersicoides IL (LA3875), and 16 other accessions.

other accessions.
Cherry tomato accessions (var. cerasiforme) were mainly classified into Group 4, except for LA1338, which was classified into Group 1. Group E consisted of 30 wild accessions belonging to S. habrochaites, S. pennellii, S. cornerliomulleri, S. chilense, S. peruvianum, S. arcanum, S. huaylasense, and S. chmielewskii. Group F consisted of 18 accessions including most accessions (13) of S. pimpinellifolium (except LA2633, which was classified into Group 1), 3 of S. lycopersicum, 1 of S. galapagense, and 1 of S. cheesmaniae (Figure 4a). Figure 4. (a) Delta K graph calculated using the Evanno method (top) and optimal population structure (delta K = 6) cluster plot (bottom). The optimal population number has been indicated using a red circle; (b) Kinship heatmap plot and dendrogram for 290 accessions drawn using the genotype data of 11,912 SNPs. The genetic distance between accessions was calculated using the VanRaden algorithm. Figure 4. (a) Delta K graph calculated using the Evanno method (top) and optimal population structure (delta K = 6) cluster plot (bottom). The optimal population number has been indicated using a red circle; (b) Kinship heatmap plot and dendrogram for 290 accessions drawn using the genotype data of 11,912 SNPs. The genetic distance between accessions was calculated using the VanRaden algorithm.
The genetic relationships of the accessions obtained using Genome Association and Prediction Integrated Tool (GAPIT) software were further confirmed based on a kinship plot and dendrogram (Figure 4b). The 290 accessions were divided into three clusters; Clusters 1 and 3 mainly consisted of S. lycopersicum and cherry tomato (var. cerasiforme) accessions, while most wild species and S. pimpinellifolium accessions that belonged to Groups 5 and 6, respectively, in the STRUCTURE analysis, were grouped into Cluster 2. In particular, the kinship heatmap showed a relatively high genetic similarity among wild species accessions in Cluster 2. Also, Clusters 1 and 3 were further divided into 3 (1-1, 1-2, and 1-3) and 2 (3-1 and 3-2) Subgroups, respectively. Subgroups 1-3 consisted of 55 accessions out of 63 accessions in pop 4, and 1 and 3 accessions belonging to pop 1 and 3, respectively, while there were mixtures of the remaining accessions in pop 1 and accessions in pop 2, 3, and 4 in Subgroups 1-1 and 1-2, respectively. Subgroup 3-1 consisted of 12 accessions of pop 1 and 2 accessions of pop 6, while in Subgroup 3-2, the three accessions of S. pimpinellifolium, S. cheesmanii, and S. galapagense, as well as five cherry tomatoes (var. cerasiforme) accessions were included.

GWAS for PMR
GWAS was first performed using 11,912 SNPs and 290 accessions, based on the mixed linear model (PC matrix + kinship matrix), in Traits Analysis by aSSociation, Evolution and Linkage (TASSEL). In the quantile-quantile plot, the actual tested P-value of the SNPs corresponding to −log 10 (P) was smaller than the predicted p-value, indicating a traitgenotype association (Figure 5a). Additional GWAS was performed using GAPIT, with the same analysis model as that used in TASSEL. The GWAS results obtained using the software have been presented as a Manhattan plot (Figure 5b).   A total of 20 SNPs significantly associated with PMR at −log 10 (P) > 3 were detected on chromosomes 1, 4, 6, 8, and 12 in TASSEL, while 14 significant SNPs were found on chromosomes 1, 2, 4, 6, 8, 10, 11, and 12 in GAPIT (Table S4). Of these, 9 SNPs were commonly detected in both the analysis programs, which were located on chromosomes 1, 4, 6, and 8, and 12; these locations were designated as QTL-1, -2, -3, -4, and -5, respectively, in this study ( Figure 6 and Table S5). Unlike other QTLs in which one SNP of common significance was detected, five common SNPs (AX-95803738, AX-95799308, AX-95776428, AX-95771531, and AX-95783448) were located in the QTL-2 region. Among the significant common SNPs, the most significant one predicted by GAPIT, AX-95771531 [−log 10 (P) = 4.41], was located in QTL-2 of chromosome 4, while the most significant one predicted by TASSEL, AX-95810925 [−log 10 (P) = 4.57], was located in QTL-4 of chromosome 8 (Table S5).  Blue and green represent the phenotypic distribution for the homozygous genotypes, while red represents that for the heterozygous genotype.

Discussion
PM is a common disease that affects a wide range of host plants worldwide. To date, GWASes have been conducted on several crops such as barley, wheat, and oats, to explore PMR-related QTLs [26][27][28]. In tomato, PM is caused by O. neolycopersici.
Tomato disease resistance has been mostly found in wild species and has been utilized for breeding [29]. Genes exhibiting resistance to major diseases induced by various  (Table S6). Among those, two nucleotide-binding site-leucine rich repeat (NBS-LRR) genes (Solyc04g015630.2.1 and Solyc04g015640.1.1), which are known as genes involved in disease resistance (R gene), were identified at the distances of 48,380 bp and 47,495 bp from the SNPs AX-95771531 and AX-95810925, respectively. No SNPs of the Axiom ® Tomato Genotyping Array were located within these NBS-LRR class genes. There were nine genes in the QTL-4 region, and the gene harboring AX-95810925 was another R gene, receptor-like serine/threonine-protein kinase (STPK) gene (Solyc08g074980.4.1) (Table S6). Two SNPs, AX-95809404 (A/G) and AX-95781958 (T/C) were located in the STPK gene exon. These SNPs were converted to CAPS markers and the tomato accessions used in the GWAS were genotyped for it, to assess the association between the marker genotype and phenotype (Table S7). In the case of AX-95781958, the accession carrying the marker genotypes of TT, TC, and CC showed PDI averages of 35.68%, 50.66%, and 67.90%, respectively, while the same were 36.71%, 58.99%, and 66.74% for the AX-95809404 marker genotypes of AA, AG, GG, respectively, thereby indicating significant differences in the level of PMR depending on the marker genotypes of the STPK gene ( Figure 7). Figure 6. A schematic diagram of tomato chromosomes, with the genomic locations of the nine single nucleotide polymorphisms (SNPs, in green) significantly associated with powdery mildew resistance (PMR) in tomato depicted. The PMR loci and DNA markers identified in previous studies are shown in blue and red, respectively. Blue and green represent the phenotypic distribution for the homozygous genotypes, while red represents that for the heterozygous genotype.

Discussion
PM is a common disease that affects a wide range of host plants worldwide. To date, GWASes have been conducted on several crops such as barley, wheat, and oats, to explore PMR-related QTLs [26][27][28]. In tomato, PM is caused by O. neolycopersici.
Tomato disease resistance has been mostly found in wild species and has been utilized for breeding [29]. Genes exhibiting resistance to major diseases induced by various biotic stress factors such as tomato yellow leaf curl virus resistance genes Ty1,3, Ty2, Ty4, ty5, and Ty6; fusarium wilt resistance genes I-1, I-2, I-3, and I-7; late blight resistance genes Ph-1, ph-2, ph-3, and Ph-5; bacterial spot resistance genes Rx-4 and RXopJ4; and nematode root-knot resistance gene Mi-1 have been derived from wild species [30][31][32][33][34][35]. Except for the Ol-6 gene, which is found in the breeding line but whose origin of tomato PMR is unknown, all other PMR genes are derived from wild species (S. lycopersicum var. Blue and green represent the phenotypic distribution for the homozygous genotypes, while red represents that for the heterozygous genotype.

Discussion
PM is a common disease that affects a wide range of host plants worldwide. To date, GWASes have been conducted on several crops such as barley, wheat, and oats, to explore PMR-related QTLs [26][27][28]. In tomato, PM is caused by O. neolycopersici.
Tomato disease resistance has been mostly found in wild species and has been utilized for breeding [29]. Genes exhibiting resistance to major diseases induced by various biotic stress factors such as tomato yellow leaf curl virus resistance genes Ty1,3, Ty2, Ty4, ty5, and Ty6; fusarium wilt resistance genes I-1, I-2, I-3, and I-7; late blight resistance genes Ph-1, ph-2, ph-3, and Ph-5; bacterial spot resistance genes Rx-4 and RXopJ4; and nematode root-knot resistance gene Mi-1 have been derived from wild species [30][31][32][33][34][35]. Except for the Ol-6 gene, which is found in the breeding line but whose origin of tomato PMR is unknown, all other PMR genes are derived from wild species (S. lycopersicum var. cerasiforme, S. peruvianum, S. habrochaites, and S. neorickii). Our bioassay results also indicated that PMR is mainly found in wild species and appears very rarely in cultivated ones [5]. This is the first study to report PMR among the wild species accessions of S. chilense LA1963 and LA2932, S. chmielewskii LA2663, S. cornerliomulleri LA1274, S. galapagense LA1141, S. pennellii LA0716, LA1340, LA1272, LA1674, and LA1809, and S. pimpinellifolium LA2181 and LA2184, and thereby, highlights them as of important value for use as genetic resources for PMR breeding. According to Lindhout et al. [36], the wild species S. chilense and S. peruvianum, which belong to the peruvianum complex have a crossing barrier with cultivated species, whereas S. pimpinellifolium, S. pennellii, and S. habrochaites, which belong to the esculentum complex, have high cross-compatibility. Therefore, S. pimpinellifolium accessions LA2181 and LA2184, S. pennellii accessions LA0716, LA1272, LA1340, LA1674, and LA1809, and S. habrochaites accessions PI126445, PI134418, and PI308182, which have been identified as PM-resistant accessions in this study, can be used as a good material for crossbreeding with elite cultivars lines.
Similar to the previously reported tomato PMR gene (Ol)-associated marker test for the collected genetic resources, the resistance marker type was mostly found in the wild accessions that showed resistance, but not in the cultivars that showed sensitivity. This suggests that the previously reported PMR-linked markers can be utilized for MAS when introducing resistance genes from these wild species accessions into cultivars through crossbreeding. In addition, there was no resistance marker type found in the A1161, A1162, A1164, and A1216 lines, which are rare resistant cultivars. It can be speculated that the resistance of these cultivars is regulated by a new PMR gene(s) or allele that occurred during the domestication of tomato cultivars and has not been reported so far.
It can be observed from the results of this study that in some wild species lines, a band size different from the size expected for an existing PMR marker was not observed, or PCR amplification was not performed at all. This trend was observed regardless of species match with the accession used for resistance gene mapping. Cases of observation of amplification products at sizes different from the expected size have been reported in the study by Singh et al. [37], which screened diverse genetic resources for rice blast resistance gene(s) using existing resistance markers. It is presumed there was a nucleotide sequence variation in the primer sites or restriction enzyme recognition sites between the line used in this study and the line used for the existing resistance gene mapping. In addition, this non-specific amplification was mainly observed as most of the PMR-linked, non-genebased markers tested in this study were CAPS converted from restriction fragment length polymorphism markers [6][7][8]. However, the test results of the two PMR gene-based ol-2 markers (M/Slmlo1, dCAPS_Slmlo1.1) produced only the amplification product of the expected size from all accessions, which demonstrates that gene-based marker development is important for more accurate and consistent MAS targeting various genetic resources.
The SNP detection probe for the Axiom ® Tomato Genotyping Array was developed based on the genome sequencing of a commercial F1 cultivar (S. lycopersicum) [21]. To the best of our knowledge, this is the first study that has reported the use of the Axiom ® Tomato Genotyping for SNP genotyping of diverse wild species tomatoes. In this study, the marker genotype mismatch between the CAPS/dCAPS assay and the SNP chip was mainly observed in wild species, which could be due to a genotyping error caused by a weak hybridization reaction of the probe to the wild species genome, and thereby, lowering of the SNP calling signal intensity. This assumption is supported by the result that wild species accessions tended to have low DQC, low call rate, and high heterozygosity rate, and also the fact that the lower the homology between the reference genome assembly, the wild species, and the probe sequence, the higher probability that the SNP-chip genotype was called out as missing or heterozygosity. Nevertheless, the clear distinction between species in the population structure analysis shows that Axiom ® Tomato Genotyping Array can be fully utilized for genome-wide SNP genotyping of various wild tomato species and further genetic analysis based on that. However, for reliable analysis results, rigid SNP filtering, such as the removal of SNPs with a high heterozygosity rate, as performed in this study, is essential.
In this study, 5 QTL regions associated with tomato PMR and 9 novel SNP markers located therein were finally selected by means of GWAS. Among these 9 SNPs, the maximum −log 10 (P) value was between 4 and 5, which is significantly low in comparison to that observed in other disease resistance GWAS cases, including PMR in barley [−log 10 (P) = 9.69] [27], blight resistance in pepper (11.23) [38], and gray leaf spot resistance in tomatoes (19.96) [39]. In general, a statistically high correlation analysis (significant p-value threshold = 5 × 10 -8 ) using GWAS is possible when the gene of the target trait exists at a high frequency in the form of a common variant (allele) within a large population of core collection [40]. Two major factors lower the GWAS analysis ability of this study: the tomato PMR is a result of a rare allele with a low frequency in some wild species, and the resistance locus is also found at several genomic locations as QTLs [5]. Several recent studies have reported the strategies that could be adopted to set a genome-wide significant p-value threshold in GWAS, depending on the genetic structure and characteristics of the population (linkage disequilibrium, homogeneity, and ancestry) and the frequency of the target gene within the population [40,41].
In this study, genes annotated as NBS-LRR and receptor-like STPK were selected as potential candidate genes for PMR. The gene encoding the NBS-LRR protein is a representative disease resistance gene (R gene) that activates a pathogen-specific defense system together with receptor-like kinase [42][43][44]. In the QTL-2 region of chromosome 4 detected in this study, two NBS-LRR genes (PMR-NBS1 and PMR-NBS2) were clustered at very close positions, thereby indicating that this could be a potential genomic region involved in PMR. The previously reported PMR loci Ol-4 and Ol-6, which induced rapid resistance to O. neolycopersici through a single cell hypersensitive reaction, have been also presumed to be a homolog to the Mi-1 locus in which a number of R genes encoding NBS-LRR proteins form a cluster [5,6]. Other cases of the NBS-LRR gene clusters related to PMR include Mla in barley [45], REN1 in grape (Vitis vinifera) [46], and MtREP1 in alfalfa (Medicago sativa) [47].
In addition to NBS-LRR, receptor protein kinase, another protein taxon involved in plant disease resistance, and receptor-like STPK belonging to the receptor-like kinase family, are also involved in signal transduction processes required for disease resistance, through phosphorylation of serine and threonine residues [43]. In the case of the bacterial speck resistance gene Pto encoding serine/threonine kinase in tomatoes, Pto kinase recognizes the elicitor generated by the avrPto of the pathogen Pseudomonas syringae pv. tomato and activates another serine/threonine kinase, Pti1, in addition to interacting with the NBS-LRR gene, Prf, to induce resistance responses such as hypersensitivity reactions and oxidative bursts [48][49][50][51]. Furthermore, it has been shown that the receptor-like serine/threonineprotein kinase gene (Stpk-V) cloned from wheat (Dasypyrum vilosum) is located within the PMR locus Pm21 and reduces haustoria formation of the PM pathogen Blumeria graminis f. sp. tritici [52]. In our study, the AX-95810925 in chromosome 8, which showed the highest −log 10 (P) value in TASSEL, was located in the intron region of the gene encoding RSTK. Moreover, there was a clear distinction in the phenotypic distribution based on the genotypes of SNPs AX-95781958 and AX-95809404 found in the exons of this gene, thereby suggesting that the SL-STPK gene may be a putative candidate for PMR in tomato.
PMR gene cloning is essential for rapid breeding processes based on the plant transformation. The first cloned PMR susceptibility (S) gene, Slmlo1 for ol-2, is used for targeted mutagenesis such as RNAi-silencing and CRISPR/Cas9 technology [53][54][55]. These Slmlo1transgenic plants demonstrated improved PMR and can be useful breeding resources for expedited development of new PMR cultivars. In addition, a putative gene (ShORR-1) for Ol-1 was recently identified by complementary cDNA-amplified fragment length polymorphism (cDNA-AFLP) analysis of the resistant cultivar S. habrochiates G1.1560, carrying the Ol-1, which can be another target gene for genome editing mutagenesis [56].
In conclusion, PM bioassay and PMR-linked marker tests performed on a natural population with a diverse genetic background provided information about genetic resources that could be useful in MAS for PMR. Our GWAS identified nine significant SNPs from five novel QTL regions, which led to the selection of 3 PMR candidate genes from the two QTL regions where the most significant SNPs were located. The disease resistance locus information and resistance-related SNPs identified in this study will contribute to understanding the molecular basis of PMR and developing new PM-resistant cultivars. Further studies are needed to carry out gene-function analysis to verify the association of the putative candidate genes identified in this study with PMR.

Plant Materials
The genetic resources used for the PMR bioassay included a total of 295 accessions (Table S1)

Pathogen Isolation and Genomic DNA Extraction
A spore suspension was prepared for the isolation and identification of PM pathogens that occurred naturally in the greenhouse at Pusan National University in 2020. The PM spores were separated from those leaves of "Moneymaker" that showed symptoms of PM, by immersing them in double distilled water and then collecting the water after filtering it out using gauze. Thereafter, the suspension was placed in a 2.0 mL tube, centrifuged at 13,000× g for 1 min, following which the supernatant was discarded and the spores collected at the bottom were collected.
Thereafter, 400 µL of sodium dodecyl sulfate (SDS) was added to the spore pellet and incubated with it for 50 min, in a water bath maintained at 65 • C, to release the genomic DNA. To separate the cell debris and SDS solution from the DNA, 200 µL of 7.5 M ammonium acetate was added to the solution, mixed evenly, left on ice for 15 min, and centrifuged at 13,000× g for 10 min. Following that, 400 µL of the supernatant was collected and transferred to a new 1.5 mL tube. For DNA precipitation, 2.5 µL of glycogen (10 mg/mL) and 600 µL of isopropanol were added to the supernatant, and it was centrifuged at 13,000× g for 10 min. The obtained supernatant was discarded, and 300 µL of 70% ethanol was added to the generated DNA pellet, to remove the residual salt. The obtained solution was again centrifuged at 13,000× g for 10 min. After the DNA pellet was completely dried, it was resuspended in 60 µL of Tris-EDTA buffer, and quantified with a NanoDrop™ 1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).
PCR was carried out using the Veriti™ 96-Well Thermal Cycler (Thermo Fisher Scientific) and the following PCR conditions; 1 cycle of pre-denaturation at 94 • C for 1 min, 35 cycles of denaturation at 94 • C for 1 min, annealing at 65 • C for 1 min, and extension at 72 • C for 2 min, followed by 1 cycle of final extension at 72 • C for 7 min. The PCR amplicons were electrophoresed on a 1% agarose gel, at 200 V for 45 min, and purified from the gel using the Expin™ Gel SV Gel Extraction Kit (GeneAll, Seoul, Korea). Direct sequencing of the purified PCR amplicons was performed using the dye-terminator method at Genotech (Daejeon, Korea).

Pathogen Identification
The ITS region sequence of PNU-PM was analyzed for similarity to 27 ITS sequences of different fungal species used for phylogenic analysis in Kiss et al. [24] and the L. taurica ITS sequence (accession number: MH698492.1) obtained from NCBI (https://www.ncbi. nlm.nih.gov/ accessed on 15 August 2022), by means of pairwise/multiple alignments using Clustal W [60] in Molecular Evolutionary Genetics Analysis version X (MEGA version X) [61]. Phylogenetic analysis of the PM pathogens was based on the neighborjoining method [62], with the Jukes-Cantor model [63] set as the nucleotide substitution model. The phylogenetic tree branch was created through 1000 iterations of the bootstrap method. In addition, the ITS sequence of PNU-PM was subjected to a BLASTn search (https://blast.ncbi.nlm.nih.gov/Blast.cgi accessed on 15 August 2022), to confirm the pathogen species with the highest concordance rate.

PM Bioassay
PM bioassay was performed in March 2020 in the greenhouse of Pusan National University. For a total of 348 accessions, 10 seeds per accession were sown into 50-cell trays filled with seedling medium (horticultural medium no. 2, Nongwoo Bio, Suwon, Korea) and the seedlings were grown under the conditions of 26 ± 5 • C and 40-70% relative humidity. The spore suspension for pathogen inoculation was prepared as described in Section 4.2.1. From 4 weeks after sowing, a spore suspension of 2.4 × 10 4 spores/mL concentration was sprayed 3 times, at intervals of 2 weeks, using a manual sprayer, to the extent that the plants were soaked.
The disease severity index (DSI) was estimated by observing PM symptoms from 1 to 6 true leaves, in the 3rd week after completion of the last inoculation. The DSI was established based on the following criteria. 0: leaf area symptoms in up to 6 true leaves = 0%, 1: <10% leaf area symptoms, 2: 10-30% leaf area symptoms, 3: >30% leaf area symptoms ( Figure 8)  PCR was performed using a Veriti™ 96-Well Thermal Cycler and the following conditions: 1 cycle of pre-denaturation at 95 °C for 1 min, 35 cycles of denaturation at 95 °C for 30 s, annealing at Tm (°C) in Table 4 for 30 s, and extension at 72 °C for 1 min, followed by a final extension at 72 °C for 7 min. In case of touchdown PCR, after pre-denaturation at 95 °C for 2 min, there were 10 cycles of denaturation at 94 °C for 15 s, annealing at 60 °C (-0.5 °C/cycle) for 30 s, and extension at 72 °C for 1 min, followed by 30 cycles of denaturation at 94 °C for 15 s, annealing at 55 °C for 30 s, and extension at 72 °C for 1 min, and final extension at 72 °C for 3 min. For sequence-characterized amplified region markers, the amplicons were electrophoresed immediately after PCR, while for CAPS markers, the amplicons were treated with the restriction enzyme, according to the manufacturer's manual, and then electrophoresed. Electrophoresis was performed on a 2% agarose gel, at 190 V for 2 h, and after gel staining with ethidium bromide, the PCR amplicons were visualized under UV light using a gel image analysis system (CoreBio i-MAX™, Davinch-K, Seoul, Korea). The marker genotypes of the accessions used in the marker assay were determined by comparing the resistance-specific and sensitivity-specific band size information reported for each marker with the band size observed in this experiment (Table 4). Table 4. List of DNA markers for powdery mildew resistance in tomato reported in previous studies. The markers used for the genotyping tomato germplasm were highly resistant or susceptible to For each accession, while referring to the PDI values of MM and PMS as controls, PDI < 10% was scored as resistance (R), PDI in the range of 10-30% as moderate resistance (MR), and PDI > 30% as susceptibility (S).

Evaluation of PMR-Linked Molecular Markers
A total of 56 accessions, including 22 that were highly resistant (R), 28 highly sensitive (S), and the PMR and PMS controls were genotyped for a total of 10 previously reported PMR-linked markers (Table 5). For PCR amplification, a PCR mixture of 10 µL total volume, containing 1 µL of template DNA (20 ng/µL), 6.7 µL of double distilled water, 0.2 µL of 10 mM dNTPs, 1 µL of 10× buffer, 0.5 µL each of the 10 pmol forward/reverse primers, and 0.1 µL Taq polymerase (5 U/µL), was used per sample. PCR was performed using a Veriti™ 96-Well Thermal Cycler and the following conditions: 1 cycle of pre-denaturation at 95 • C for 1 min, 35 cycles of denaturation at 95 • C for 30 s, annealing at Tm ( • C) in Table 5 for 30 s, and extension at 72 • C for 1 min, followed by a final extension at 72 • C for 7 min. In case of touchdown PCR, after pre-denaturation at 95 • C for 2 min, there were 10 cycles of denaturation at 94 • C for 15 s, annealing at 60 • C (-0.5 • C/cycle) for 30 s, and extension at 72 • C for 1 min, followed by 30 cycles of denaturation at 94 • C for 15 s, annealing at 55 • C for 30 s, and extension at 72 • C for 1 min, and final extension at 72 • C for 3 min. For sequence-characterized amplified region markers, the amplicons were electrophoresed immediately after PCR, while for CAPS markers, the amplicons were treated with the restriction enzyme, according to the manufacturer's manual, and then electrophoresed. Electrophoresis was performed on a 2% agarose gel, at 190 V for 2 h, and after gel staining with ethidium bromide, the PCR amplicons were visualized under UV light using a gel image analysis system (CoreBio i-MAX™, Davinch-K, Seoul, Korea). The marker genotypes of the accessions used in the marker assay were determined by comparing the resistance-specific and sensitivity-specific band size information reported for each marker with the band size observed in this experiment (Table 5).

SNP Genotyping Using the Axiom ® Tomato Genotyping Array
For SNP genotyping, for use in the GWAS analysis, genomic DNA was extracted for the tomato accessions used in the PM bioassay. Fresh leaves of 3rd-4th true leave stage were collected from three individual seedlings per accession and pooled. Genomic DNA extraction was conducted using the SDS extraction method, with a partial modification, following which the quantity and quality of the purified DNA were measured as described in Section 4.2.1. All DNA samples were diluted to the final concentration of 20 ng/µL, and genotyping of a total of 51,214 SNP loci using the 55K Axiom ® Tomato Genotyping Array (Thermo Fisher Scientific) was performed by DNA Link Inc. (Seoul, Korea).

SNP Genotyping Data QC and SNP Filtering
In order to check whether the wild species accessions used in this study are well applied to the Axiom ® Tomato Genotyping Array, 125 cultivated species and 57 wild species accessions were selected and the DQC, call rate, and heterozygosity rate for the 51,214 raw SNPs in them were analyzed. Correlation analysis between these three factors was carried out by calculating the Pearson correlation coefficient and significance probability (p-value) using SPSS Statistics 25.0.0 (IBM Inc., Chicago, IL, USA).
To confirm sequence homology, DNA sequences corresponding to the 70 bp probe sequence for each of the 51,214 SNP loci of the 55K Axiom ® Tomato Genotyping Array were obtained from the reference genome assembly of S. lycopersicum and five wild species in NCBI (Table S8). Sequence alignment was performed using the Bowtie2 program [65], under the following conditions: maximum mismatches in seed alignment(-N) = 1, length of seed substrings(-L) = 20, and the interval between seed substrings w/r/t read len(i) = S, 1, 0.50. Then, for each SNP locus in the Axiom ® Tomato Genotyping Array, the correlation of the ratio of missing data and heterozygosity with the degree of variation between probe and reference genome sequence was calculated in the 169 accessions (125 from S. lycopersicum, 8 from S. chilense, 8 from S. habrochaites, 10 from S. pennellii, 18 from S. pimpinellifolium, and 3 from S. arcanum 3) that were used in SNP genotyping data QC.

SNP Conversion to CAPS and dCAPS Markers
To verify the reliability of the SNP genotyping results, six of the SNPs genotyped on the Axiom ® Tomato Genotyping Array were selected and converted into CAPS or dCAPS markers ( Table 6). The physical location of the selected SNPs was searched for in the SL4.0 assembly of SGN, to extract the flanking sequence. For CAPS and dCAPS marker conversion, restriction enzymes capable of distinguishing polymorphisms were confirmed using NEB cutter v2.0 (http://nc2.neb.com/NEBcutter2/index.php accessed on 15 August 2022) and dCAPS finder 2.0 (http://helix.wustl.edu/dcaps/ accessed on 15 August 2022), respectively. The primers were designed using Primer3 v0.4.0. (https: //bioinfo.ut.ee/primer3-0.4.0/ accessed on 15 August 2022). PCR, restriction enzyme treatment and electrophoresis were performed in the same manner as that for the CAPS marker, as described in Section 4.4, except that PCR was performed with annealing at 58 • C, in a non-touchdown method. GWAS for PMR was performed using the PM bioassay and SNP genotype data from 290 accessions (244 accessions of S. lycopersicum and 46 accessions of wild species) that showed relatively low standard deviation in the bioassay (Table S1). SNP filtering from the 51,214 raw SNPs was performed by applying the conditions of minor allele frequency > 5%, missing proportion <5%, and heterozygosity proportion <10% (only for wild species).

Population Structure Analysis
Population structure analysis was performed based on the admixture model of STRUC-TURE v2.3.4 [66] using the 290 accessions and the filtered SNPs selected for GWAS. In the K (number of population) range of 1 to 1000 burn-in periods and 10,000 Markov Chain Monte-Carlo (MCMC) for each K were repeated 10 times. To determine the optimal size of the population, the value of ad hoc statistics (∆K) was calculated using the Evanno method [67] with STRUCTURE HARVESTER (http://taylor0.biology.ucla.edu/ structureHarvester/ accessed on 15 August 2022). The genetic relationship among the accessions was analyzed using the VanRaden kinship algorithm [68] on the R-based package GAPIT, and confirmed through a Kinship heatmap plot and dendrogram.

GWAS Analysis
GWAS was performed based on the mixed linear model (PC matrix + kinship matrix) and applied to the analysis using two programs, Traits Analysis by aSSociation, Evolution and Linkage (TASSEL) 5.0 [69] and R-based package Genome Association and Prediction Integrated Tool (GAPIT) ver. 3 [70]. The results from both programs were visualized and confirmed with a Manhattan plot. From the analysis result, SNPs with a p-value < 0.001 were considered significant SNPs associated with PMR, and only SNPs that were common to both TASSEL and GAPIT software were reselected and used for subsequent analysis.

Candidate Gene Selection
For the identification of putative candidate genes for PMR, genes located within 50 kb on each flaking region of the SNPs with the highest −log 10 (P) value in TASSEL and GAPIT were selected from the SL4.0 assembly (ITAG4.0) in SGN. The genes related to disease resistance were then selected as candidate genes based on the gene function annotation information.
Alle experimental procedures described in the section '4. Materials and Methods' were schematically represented in Figure S1.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms232113610/s1, Table S1: Phenotypic data of 295 tomato accessions obtained from powdery mildew (PM) (Oidium neolycopersici) disease assay; Table S2: Dish QC (DQC), call rate, heterozygosity rate information for each sample of the 182 accessions used for assessing the quality of results for the single nucleotide polymorphism (SNP) genotyping by the 55K Axiom ® Tomato Genotyping Array; Table S3: Genotypes of 11,912 single nucleotide polymorphisms (SNPs) for the 290 tomato accessions, revealed using the 55K Axiom ® Tomato Genotyping Array; Table S4: Single nucleotide polymorphisms (SNPs) selected for significant association [−log 10 (P) > 3] with powdery mildew resistance (PMR) in tomato, upon carrying out a genome-wide association study (GWAS); Table S5: Genotypes in 290 tomato accessions for nine single nucleotide polymorphisms (SNPs) identified to be significantly associated with powdery mildew resistance (PMR) in the genome-wide association study (GWAS); Table S6: List of annotated genes in the tomato powdery mildew resistance (PMR)-associated QTL-2 (chromosome 4, 5.87-6.07 Mb) and QTL-4 (chromosome 8, 57.17-57.27 Mb) regions containing two significant SNPs, AX-95771531 and AX-95810925, respectively; Table S7: The genotyping results reconfirmed by converting chip single nucleotide polymorphisms (SNPs) AX-95809404 and AX-95781958 into cleaved amplified polymorphic sequence (CAPS) markers to the 290 accessions used in the genome-wide association study (GWAS); Table S8: Tomato reference genome assembly information used for similarity analysis with the probe sequences of the 55K Axiom ® Tomato Genotyping Array; Figure S1: A schematic representation of experimental procedure.