A total of 384 DH lines belonging to four DH populations were developed from four elite inbred lines ( Table 1 ) in BeiJing (N40°08’ E116°10’) in 2017. The founders of the four DH populations were C783 × C229, C783 × UH306, C783 × EH, C229 × UH306, respectively, and we named them as POP1-4. The quantities of DH lines in POP1-4 are 66, 107, 127, and 77, respectively. Then, we testcross each population with three testers, yielding a total of 903 hybrids ( Table 1 ). Thus, the hybrid population is divided into 12 subgroups with quantities ranging from 40 to 119.

To test the SCR infection levels of the accessions, we planted the DH lines and hybrids in Huang-Huai-Hai summer corn planting region in China for phenotypic identification. This region is the main area where SCR occurs in China due to high temperature and rainy summer. The DHs were planted in Jinan (N37°42’ E117°27’) in 2021, Xinxiang (N35°9’ E113°47’) in 2018 and 2021; the hybrids were planted in Jinan in 2021, Jining (N35°6’ E116°31’) in 2020 and 2021, Xinxiang in 2020. We used the augmented experimental design, setting every 20 accessions as a block. Each block consisted of 19 rows and standard accessions were planted in random order. For the DHs, the standard accession was a susceptible inbred line C116A. For the hybrids, the standard accession was a susceptible commercial hybrids ZhengDan958. In the field, each accession was planted in a one-row plot for DHs and a two-row plot for hybrids at a spacing with 0.6 x 0.25 m spacing (66,000 plants per hectare).

SCRRS for each accession was visually collected from the leaf area covered by lesions at 4 weeks after flowering ( Figure 1A ). A rating scale of 1 corresponds to severe infection covering > 75% of the leaf surface, 3 corresponds to moderate-to-severe infection covering 50–75% of the leaf surface, 5 corresponds to moderate infection covering 25–50% of the leaf surface, 7 corresponds to weak to moderate infection covering 10–25% of the leaf surface, and 9 corresponds to high resistance covering 0–10% of the leaf surface (Ren et al., 2021).

The raw phenotypic data were analyzed using the linear mixed model with an R add-on package “lme4” (Bates et al., 2014). Best linear unbiased predictors (BLUPs) were calculated for DHs and hybrids. In the model,

where, y_ij is the mean phenotypic value of the ith DH or hybrid in the jth environment;

μ is the overall mean of the trait; g_i is the random effect of the ith accession; l_j is the random effect of the jth environment; g_i ×l_j is the random interaction effect between the ith accession and the jth environment; and ϵ is the random error.

Heritability was calculated using variance components estimated from the above model. The following equation was used to estimate heritability on an individual plot basis,

Where V_g is the genotypic variance component, V_e is the error variance, and l is the number of environments.

Young leaves of all the DHs and the tester lines were sampled for DNA extraction using the CTAB method (Porebski et al., 1997). Then, genotyping was conducted using the Maize-6H-60K SNP chip (Tian et al., 2021). SNPs with minor allele frequency (MAF) > 0.05 and per locus missing rate< 0.1 were filtered out using plink 1.90 (http://www.cog-genomics.org/plink2/). The genotypes of hybrids were obtained with the cleaned SNPs (N=34,037) of DHs and testers using TASSEL V5.2 software (Bradbury et al., 2007). Pairwise measures of linkage disequilibrium (LD) were performed to analyze the squared allele‐frequency correlation coefficient (r²) between two loci using plink software. Only SNPs with a MAF > 0.05 and less than 0.1 missing data were used to estimate LD. Principal component analysis (PCA) was used to assess the level of genetic structure using TASSEL software.

Genome wide association analysis was performed with the BLUPs obtained from the combined analysis for the DHs and hybrids. A Fixed and Random Model Circulating Probability Unification (FarmCPU) method, as proposed by Liu et al. (Liu et al., 2016) was applied in GAPIT V3 software (Wang and Zhang, 2021). Two genetic models, additive and dominant, are used for the hybrids panel, and only the additive model was used for the DH panel. Under the additive model, homozygous genotypes with recessive allele combinations were coded as 0, homozygous genotypes with dominant allele combinations are coded as 2, and heterozygous genotypes were coded as 1. Under the dominant model, both types of homozygous genotypes are coded as 0 and heterozygous genotypes were coded as 1. The Bonferroni testing was used to determine the genome-wide significance thresholds (0.05/34,034 = 1.47 × 10⁻⁶), where 34,037 is the total number of SNP markers (Holm, 1979). Markers whose P-values passed the threshold were identified as candidate loci. Unlike natural material populations, such as the artificial DH population or testcross hybrid population, which had a high LD level, our candidate intervals were selected according to LD decay and LD block. Makers with a physical distance of<20 Mb and in high LD (r² ≥ 0.8) were considered to mark the same genomic region (Mayer et al., 2020). The corresponding candidate region was described by the positions of the first and last maker, respectively.

The Genomic prediction was performed for the hybrid panel with three conditions, including 1) extended GBLUP models, 2) maker based prediction methods, and 3) prediction models with QTL calculated by GWAS as fixed effects.

For extended GBLUP models, which comprised additive (G_a ), dominant (G_d ) and epistatic relationship matrices. G_a and G_d matrices were calculated using the “sommer” package in R (Covarrubias-Pazaran, 2016). The epistatic matrices terms were computed using Hadamard products (i.e., cell-by-cell product denoted as “#”) of the following form: (i) additive-by-additive interactions (G_a#G_a ); (ii) dominance-by-dominance interactions (G_d#G_d ); and (iii) additive-by-dominance interactions (G_a#G_d ), respectively (Muñoz et al., 2014). In total, six GBLUP models were used in this study ( Table 1 ). The programs were implemented in the “BGLR” package in R (Pérez and De Los Campos, 2014). The extended GBLUP models can be described as

where y is the vector of phenotypic data; 1_n is the n-dimensional vector of ones; μ is the overall mean; u_a, u_d, u_aa, u_ad, u_dd are the vectors of random effects for additive, dominant, additive-by-additive, additive-by-dominance and dominance-by-dominance effects assumed to obey the normal distributions N(0, G a σ a 2 ), N(0, G d σ d 2 ), N(0, G aa σ aa 2 ), N(0, G ad σ ad 2 ) and N(0, G dd σ dd 2 ), respectively; G_a , G_d , G_aa , G_ad and G_dd are the genomic relationship matrices corresponding to additive, dominance, additive-by-additive, additive-by-dominance and dominance-by-dominance genotypic values, respectively.

We also performed maker based prediction models including RRBLUP (Whittaker et al., 2000), BRR (Pérez and De Los Campos, 2014), BL (Park and Casella, 2008), BayesA-C (Meuwissen et al., 2001b). The RRBLUP method is based on a restricted maximum likelihood (REML) approach to ridge regression, we performed it by R package “rrBLUP” (Endelman, 2011). Meanwhile, we also used Bayes-based methods to fit models, containing different prior densities, i.e., Gaussian (BRR), Double exponential (BL), Scaled-t (BayesA), Scaled-t mixture (BayesB), Gaussian mixture (BayesC) in BGLR package (Pérez and De Los Campos, 2014). The basic model is,

where y is the vector of phenotypes; 1_n is the n-dimensional vector of ones; μ is the overall mean,; α is a vector of random regression coefficients of all the marker effects; Z is an genotypic matrix for markers; and ϵ is a vector of residuals. The alternative methods discussed here differ primarily in their specific prior used for α . For RRBLUP, α~N(0, I σ a 2 ) and σ a 2 has a scaled inverse chi-square distribution. For BayesA, the unconditional distributions of the marker effects follow identical and independent univariate t distributions, each with mean zero. BayesB employs a mixture distribution that includes a point of mass at zero and a univariate scaled t distribution. The assumption of BayesC is that each marker effect is zero with probability π and follows a univariate normal distribution with probability (1 − π) having mean zero and variance σ j 2 , which has a scaled inverse chi-square distribution.

To further improve the prediction ability, we added QTL into the mixed linear model as fixed factors. Two representative models were selected, namely G_A and BayesB. We added additive localization maker and dominant localization maker obtained by GWAS into the model separately or together, including G_A_qa (G_A with additive GWAS SNPs), G_A_qd (G_A with dominant GWAS SNPs), G_A_qad (G_A with additive and dominant GWAS SNPs), BayesB_qa (BayesB with additive GWAS SNPs), BayesB_qd (BayesB with dominant GWAS SNPs), BayesB_qad (BayesB with additive and dominant GWAS SNPs). When the prediction was performed with additive QTL, homozygous genotypes with recessive allele combinations were coded as 0, homozygous genotypes with dominant allele combinations were coded as 2, and heterozygous genotypes were coded as 1. When the prediction was performed with dominant QTL, both types of homozygous genotypes were coded as 0 and heterozygous genotypes were coded as 1. The models can be described as,

where y is the vector of phenotypes; X_QTLa, X_QTLd, X_QTLad are incidence matrices of additive localization makers, dominant localization makers and both, respectively; β_a, β_d and β_ad are vectors of fixed effects for X_QTLa, X_QTLd and X_QTLad , respectively; G_a is the genomic relationship matrix corresponding to additive genotypic values; Z is a genotypic matrix for all markers; α is a vector of random regression coefficients of all the marker effects; and ϵ is a vector of residuals.

In this study, we used a five-fold cross validation approach to assess the ability of the tested GP models. Prediction accuracy was quantified using two methods, 1) the Pearson correlation between the input trait values and the genomic estimated breeding values (GEBVs) predicted from a given GS model evaluated in the test set, 2) the number of top 20% accessions intersections selected by GEBVs and true values derived by the total number of accessions in the test set. The process was repeated 100 times to eliminate the prediction error.

Data analysis was carried out with R software (Version 3.6.2). Microsoft Excel for Mac (Version 16.50) was used for collation of phenotypic data. Tukey’s test and Students’ t-test were performed to assess the significance of differences between values, and P < 0.05 was considered to be statistically significant.

We evaluated the SCRRS in 384 DH lines and 903 hybrids under three and four environments, respectively. The results indicated that there were abundant phenotypic variations within each panel ( Figures 1B, C ). The descriptive statistics for each population are presented in Table 2 . For DH founders, C229 and UH306 showed the highest and lowest SCRRSs, which were 5.90 and 4.89 respectively. In the DH panel, the SCRRS ranged from 4.10–7.06, and POP1 showed significantly high resistance to SCR (Tukey-test, P<0.05), with the mean SCRRS was 6.02. In the hybrid panel, the scores ranged from 3.65 to 6.08, with a mean of 4.97. The most resistant subgroup was POP2/C229, with the mean SCRRS was 5.39. In particular, the DHs and hybrids were planted at different locations, so we didn’t make a comparison between the two panels. The broad-sense heritability (H² ) analysis revealed that the H² in the DH panel and the hybrid panel were 0.64 and 0.54, respectively, suggesting that the phenotypic variation in the two panels was genetically controlled.

After marker quality control (see Materials and Methods), 34,037 SNP markers for 384 DH genotypes were available for further analysis. The 903 hybrid genotypes were imputed by their parents. The molecular diversities among the DH lines and hybrids were examined by applying principal coordinate analysis ( Figures 2A, B ). There were 4 subgroups in the DH population, among which POP2 and POP3 were relatively close, possibly because they share a common parent C783 and another parent EH were closely related to UH306. In the hybrid panel, three subgroups were observed, including hybrids using C116A as the tester, Z58 as the tester, DH founders or J2416 as the testers. The LD was estimated for the two panels using SNPs. The LD rapidly decreased with increasing the physical distance between SNPs ( Figure 2C ), but the decay rate varied among the two panels. At r²= 0.2, the mean LD decay was about 20 Mb and 5 Mb for the DH panel and the hybrid panel.

Three GWAS processes were performed using the FarmCPU method, including additive GWAS in the DH panel, additive GWAS in the hybrid panel and dominant GWAS in the hybrid panel. The quantile–quantile (q–q) plot implied that the population structure and family relatedness were well controlled in the three GWAS jobs ( Figures 3B, D, F ). One SNP (AX-107958879) on chromosome 10 significantly associated with the SCRRC trait was identified at P< 1.47 × 10⁻⁶ in the DH panel, with effect value was -0.25 ( Figure 3A and Table 3 ). LD analysis suggested candidate region was 1,150,363–3,990,150 bp, which overlapped the previously reported gene RPPC or RPPK ( Supplementary Figure 1 ). For the additive GWAS in hybrids, three significant SNPs (AX-90698604, AX-108029030, AX-108089672) were detected, with -log10 (P) ranging from 6.30 to 40.61. These SNPs were distributed on chromosomes 1, 8 and 10, with the candidate regions Chr1: 181,330,348-188,255,567, Chr8: 13,140,413-18,429,572, Chr10: 2,656,837-4,990,741, respectively ( Figure 3C and Table 3 ). The effects of them were 0.24, 0.22 and -0.54. For dominant GWAS in the hybrid panel, three SNPs (AX-107981937, AX-108109448, AX-108089672) on chromosomes 7, 8 and 10 significantly associated, with -log10 (P) were ranged 7.17-37.12, the effects were -0.16, 0.14 and -0.5 ( Figure 3E and Table 3 ). Their candidate regions were suggested as Chr7: 13,581,102-23,774,017, Chr8: 167,766,262-168,856,337, Chr10: 2,656,837-4,990,741. The QTL on chromosome 10 obtained by the three GWAS processes were identified as the same region using LD analysis ( Supplementary Figure 1 ).

The effect of marker density and training population size on the GP accuracy is shown in Figure 4 . For marker density, the prediction accuracy increased as the number of markers increased. The prediction accuracy increased rapidly when the number of markers increased from 10 to 5,000. Then, the prediction accuracy increased slightly when the number of markers kept increasing. For training population size, prediction accuracy increased as the size increased, and no slowdown in the growth rate was observed.

To meet the breeding needs of SCR-resistant hybrid selection, different GP methods were implemented to improve the prediction accuracy. Firstly, six extended GBLUP models with combinations of additive, dominant, epistatic matrices were tested ( Figure 5 and Supplementary Table 1 ). For test set correlation, the G_A model which only used the additive matrix was found significantly better than the G_D model which used the dominant matrix, with accuracy were 0.60 and 0.57, respectively. Another less accurate model than G_A, but not significant, was the G_A_D_E model, which had a mean accuracy of 0.59. The accuracy of the G_A_AA model was higher than that of the G_A model, suggesting that the epistatic effect was beneficial to GP in this study. Other models (G_A_D, G_A_AD, G_A_DD) performed as well or slightly better than G_A, with the accuracy of 0.60, 0.61 and 0.60, respectively. For top selection accuracy, the overall accuracy was lower than that of the test set correlation. The correlation test shows a significant correlation between the two accuracy evaluation methods, with R=0.77 (P value<0.05). Interestingly, G_A is better than other models for top selection, which is different from previous reports (Muñoz et al., 2014). The accuracy of the six models ranged from 0.45 to 0.48, indicating that further improvement is needed.

Then, given that our hybrid panel had several significant resistance QTL, six maker based prediction models were performed ( Figure 6 and Supplementary Table 1 ). The results showed that RRBLUP, BRR, BL and BayesC were at the same level with an accuracy of 0.60 for test set correlation. BayesA and BayesB were significantly better than other models with an accuracy of 0.65. The top selection accuracy showed the same trend, the accuracy of BayesA and BayesB were 0.53 and 0.52, respectively, which were significantly higher than other models. Meanwhile, a more significant correlation than extend GBLUP models was detected between the two accuracy evaluation methods (R = 0.98).

Two representative methods (G_A, BayesB) were selected to test the effect of adding QTL as fixation factors ( Figure 7 and Supplementary Table 1 ). For the two models, the GP results showed that the test set correlations were significantly improved whether additive GWAS QTL, dominant GWAS QTL or both were added. The test set correlations for G_A, G_A_qa, G_A_qd, G_A_qad, BayesB, BayesB_qa, BayesB_qd, BayesB_qad were 0.60, 0.66, 0.67, 0.67, 0.65, 0.66, 0.67, 0.66. For top selection using G_A models, the addition of QTL significantly improved the accuracy, among which the G_A_qa model performed best, with an accuracy of 0.55. In contrast, the BayesB model was not significantly changed the accuracy after the QTL addition for top selection accuracy, in which BayesB_qa (0.53) was slightly improved, while BayesB_qd (0.51) and BayesB_qad (0.51) slightly decreased compared with BayesB (0.52).