Two tomato germplasm collections, TGC1 (n = 162) and TGC2 (n = 191), were used in this study. The TGC1 accessions included 119 determinate breeding lines, 42 semi-determinate breeding lines, and one undetermined breeding line from a private breeding program, originating from seven different countries ( Supplementary Table S1 ). For TGC2 representing indeterminate accessions, 98 breeding lines were derived from the National Institute of Horticultural and Herbal Science (NIHHS) in Rural Development Administration (RDA), the Republic of Korea (ROK). The other 93 accessions were collected from the National Agrobiodiversity Center (NAC) in RDA, the Germplasm Resources Information Network (GRIN) in the U.S. Department of Agriculture, the C. M. Rick Tomato Genetics Resource Center (TGRC), and Sejong University ( Supplementary Table S1 ).

These collections were independently evaluated for BW resistance at the seedling stage in a greenhouse. Three seedling assays for each collection were conducted with artificial inoculation of a virulent strain WR-1 (race 1, biovar 3, and phylotype I) in different seasons (spring, summer, and fall). For inoculum, this strain was cultured in the Difco™ nutrient broth medium containing 3 g/L of beef extract and 5 g/L of peptone (BD, Sparks, MD, USA) at 28°C for 48 h. Bacterial cells were collected and resuspended in sterile, double-distilled water, and the resulting suspension was standardized to OD₆₀₀ = 0.3 (10⁸ CFU/mL) using the NanoDrop™ One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The roots of 6- to 8-week-old seedlings were wounded by cutting for inoculation and then dipped in the bacterial suspension for 30 min (Carmeille et al., 2006). For each assay, five to seven seedlings per tomato accession were inoculated and then placed in the greenhouse with a randomized complete block design. Wilting symptoms were scored 10 or 14 days after inoculation using a 1–5 scale, where 1 = no wilting symptom, 2 = one or two leaves wilted, 3 = most of the leaves wilted, 4 = all the leaves wilted, and 5 = plant died (Kelman, 1953). To correct environmental effects between assays, the BLUP values were calculated for the phenotypic data of BW resistance using the R package “lme4” (Bates et al., 2015) and used for further analysis.

Genomic DNA of each accession was isolated from fresh and young leaf tissue of 3- to 4-week-old seedlings using a modified cetyl trimethyl ammonium bromide (CTAB) method (Kabelka et al., 2002). The isolated DNA pellets were resuspended by T_1/10E buffer (10 mM of Tris–HCl pH 8.0, 0.1 mM of EDTA), and their concentrations were estimated using the NanoDrop™ One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The final concentration was adjusted to 50 ng/µL for genotyping with the Axiom^® tomato array containing 51,912 SNPs (Yamamoto et al., 2016). For this SNP array-based genotyping, 200 ng of genomic DNA for each accession was amplified and then fragmented into 25–125 bp using the Axiom^® 2.0 reagent kit (Thermo Fisher Scientific, Waltham, MA, USA). The resulting DNA fragments were hybridized to the array in the Affymetrix^® GeneTitan system according to the manufacturer’s instructions. SNP calling was conducted using the Affymetirx^® Power Tools software package v1.18, and high-quality SNPs were filtered based on <10% of missing data and >5% of minor allele frequency. The remaining missing data of these SNPs were imputed using BEAGLE v5 with default parameter settings (Browning, 2008).

Six GS models were used to predict GEBVs for BW resistance in three training populations (TGC1, TGC2, and combined). Of these, the parametric models included ridge regression-best linear unbiased prediction (RR-BLUP), BayesA, and Bayesian LASSO. The non-parametric models were reproducing kernel Hilbert space (RKHS), support vector machine (SVM), and random forest. The SNP genotypes were used to train the GS models along with the BLUP values for BW resistance in each population. The estimates of GEBVs were obtained in the GS models implemented in several R packages: “rrBLUP” for RR-BLUP and RKHS (Endelman, 2011), “BGLR” for BayesA and Bayesian LASSO (Pérez and De Los Campos, 2014), “e1071” for SVM (Meyer et al., 2023), and “randomForest” for random forest (Liaw and Wiener, 2002). Cross-validation was conducted using the leave-one-out cross-validation (LOOCV) method (Molinaro et al., 2005). Prediction accuracy was determined based on the Pearson correlation coefficients between GEBVs and observed phenotypes (BLUP values) in the training populations.

The prediction accuracy of low-density markers was investigated using the selected models, which showed the best performance in each of the two groups (parametric and non-parametric). For this analysis, several subsets were generated by filtering the total marker set of confident 31,142 SNPs. The first subset was obtained based on the marker effect values of 31,142 SNPs that were determined by estimating the effects of SNP allele substitution in a parametric model. With a non-parametric model, the second subset was generated based on variable importance (VIM) values that were determined as the percentage of increased mean squared error (MSE) after this marker was randomly permuted in a new sample (Nicodemus et al., 2010). The marker effect and VIM values were calculated in the combined population and 1,557 of 31,142 SNPs were selected as the top 5% of marker effects or VIM values. The third subset consisted of SNPs that were significantly associated with BW resistance. For this subset, a GWAS was conducted in the combined population. Population structure was inferred to determine the best K (number of clusters) using the STRUCTURE v.2.3.4 program (Pritchard et al., 2000). The 10 Ks (1–10) were tested in 10 independent simulations for each K with a burn-in of 20,000 iterations and a run length of 100,000 iterations. With the resulting log-likelihood values, the best K was found in the delta K method (Evanno et al., 2005), and the corresponding membership coefficients of tomato accessions were used as the Q matrix. In addition, a kinship matrix was generated using the VanRaden algorithm (Vanraden, 2008) for association analysis. Marker–trait association was detected using a multilocus mixed model (MLMM) implemented in the genomic association and prediction integrated tool (GAPIT) (Lipka et al., 2012). In this model, both Q and kinship matrices were used as covariates to reduce false-positive associations due to population structure and familial relatedness (Yu and Buckler, 2006). Significant associations were determined using two thresholds, P <0.05 and P <0.005. Two additional subsets (fourth and fifth) were also generated using the first to third subsets. For the fourth subset, all of the SNPs in the three subsets were combined, while the common SNPs were used for the fifth subset. These five subsets and a total marker set (31,142 SNPs) were used to estimate GEBVs for BW resistance, and their prediction accuracies were evaluated in the combined population.

Two collections (TGC1 and TGC2) were evaluated for BW resistance in three independent greenhouse trials, respectively. For TGC1 (n = 162), we found asymmetric distributions skewed toward resistant responses in all the seedling assays ( Figure 1A ). The first and third assays, which were performed in spring and fall, showed the five rating scales (1 = no wilting symptom to 5 = plant died) with means of 1.76 and 1.98, respectively. The Pearson correlation coefficient was 0.69 between these assays ( Table 1 ). The second assay conducted in summer showed a mean of 2.80 and relatively lower correlations with the first assay (0.16) and the third assay (0.36). The BLUP data of all the seedling assays showed less asymmetric distribution with a mean of 2.18 and correlation coefficients of 0.70 (vs. the second assay) to 0.87 (vs. the third assay) ( Figure 1A , Table 1 ). For TGC2, the first and third assays showed skewed distributions to susceptible responses with means of 3.93 and 3.57, while the second assay showed a relatively symmetric distribution with a mean of 2.99 ( Figure 1B ). Interestingly, the first assay showed a higher correlation with the second assay (0.64) compared with the third assay (0.35) ( Table 1 ). The first and second assays were conducted in two fall seasons, but the third assay was done in summer. In addition, the BLUP data of TGC2 showed a mean of 3.50 and correlation coefficients of 0.74 (vs. the third assay) to 0.84 (vs. the second assay) ( Table 1 ). For the combined population, additional BLUP data were generated with six phenotypic data sets of both TGC1 and TGC2 and showed normal distributions with a mean of 3.09 ( Figure 1C ). The BLUP data of the three populations (TGC1, TGC2, and combined) were used to assess the prediction accuracy of GEBVs for BW resistance.

The 51K SNP array generated common 31,142 SNPs with reliable polymorphisms in both tomato collections. Their genotypic data were used to train six GS models along with the BLUP data in three training populations (TGC1, TGC2, and combined). The parametric models, RR-BLUP, BayesA, and Bayesian LASSO, showed higher prediction accuracies in TGC2 (0.672–0.680) relative to the TGC1 (0.518–0.544) and combined (0.609–0.625) populations ( Table 2 ). Of these models, Bayesian LASSO showed the highest prediction accuracy in both TGC1 (0.544) and combined (0.625) populations. The prediction accuracy of Bayesian LASSO in TGC2 was 0.672, which was slightly lower than 0.680 for RR-BLUP. Similarly, three non-parametric models (RKHS, SVM, and random forest) revealed higher prediction accuracies in TGC2 (0.595–0.683) relative to the TGC1 (0.451–0.526) and combined (0.549–0.614) populations ( Table 2 ). In addition, random forest provided the highest prediction accuracy in both TGC2 (0.683) and combined (0.614) populations, while SVM resulted in the lowest prediction accuracy in all the populations. With these results, Bayesian LASSO and random forest were selected as the best models to predict GEBVs for BW resistance and used to identify low-density SNP sets for genomic selection.

The marker effects of 31,142 SNPs were estimated using Bayesian LASSO in the combined population and 1,557 SNPs in the top 5% were selected to generate a subset of markers ( Figure 2A ). Similarly, the second subset of 1,557 SNPs was produced using the VIM values, which were calculated in random forest ( Figure 2B ). For the third subset, GWAS for BW resistance was conducted in the combined population. Population structure analysis revealed that 353 tomato accessions were separated into seven clusters, and the number of accessions per cluster ranged from 10 (cluster 7) to 102 (cluster 6) ( Supplementary Table S1 ). The majority of accessions in clusters 1 (88.0%) and 2 (89.7%) were derived from TGC1, while clusters 6 and 7 were represented by the TGC2 accessions (89.2% and 100.0%). In the other clusters, we found relatively high levels of mixture (31.9%–68.1% for each collection). The multilocus mixed model (MLMM) identified 63 SNPs significantly associated with BW resistance at P <0.005 on nine chromosomes (1, 2, 4, 6–9, 11, 12) including three major QTL (Bwr-4, Bwr-6, and Bwr-12) ( Supplementary Figure S1 , Supplementary Table S2 ). With the threshold of P <0.05, a total of 1,357 SNPs showed significant associations across 12 chromosomes ( Figure 2C ). This number of SNPs could be large enough to capture QTL with minor effects and, thus, was used to generate the third subset.

Each subset included not only subset-specific but also common SNPs, which were found in more than two subsets ( Figure 3 ). For the first subset, 845 (54.3%) of 1,557 SNPs were subset-specific and other SNPs were also present in the second (344 SNPs) and third (503 SNPs) subsets. Of the 1,557 SNPs, 1,118 (71.8%) in the second subset were subset-specific, while 230 SNPs were also found in the third subset. Moreover, a total of 135 SNPs were found in all the subsets ( Figure 3A ). With only 63 significant SNPs at P <0.005 in the third subset, we found much smaller numbers of common SNPs in all pairwise comparisons, ranging from 24 to 344 ( Figure 3B ). Therefore, the significant 1,357 SNPs at P <0.05 were used as the third subset to generate two additional subsets. The fourth subset consisted of 3,529 that were obtained by combining all the SNPs in the first to third subsets, while the fifth subset included 135 common SNPs ( Figure 3A ).

Five SNP subsets were assessed for prediction accuracy using two GS models, Bayesian LASSO and random forest, in the combined population. With Bayesian LASSO, all the subsets showed higher prediction accuracies compared with 0.625 in the total set of 31,142 SNPs ( Figure 4 ). Of these, the first (marker effect-based) and third (GWAS-based) subsets resulted in 0.753 and 0.756 that were higher than 0.699 of the second subset (VIM-based). Interestingly, the fourth and fifth subsets, which consisted of 3,529 and 135 SNPs, showed a little difference in prediction accuracy (0.740 vs. 0.743) ( Figure 4 ). The random forest model also provided higher prediction accuracies in the five subsets (0.661 to 0.702) compared with the total set (0.614). However, the first to fourth subsets showed similar levels of prediction accuracy ranging from 0.661 to 0.682. In addition, the highest prediction accuracy was found in the fifth subset ( Figure 4 ).