In a previous study (Shu et al., 2021), 490 inbred lines collected from maize breeders at China Summer Corn Belt were genotyped via genotyping-by-sequencing (GBS); among these, 140 inbred lines that can grow and mature normally in both the summer corn region and Sanya winter nursery in Hainan Province were scored for resistance to southern corn rust (SCR). All 140 of these inbred maize lines were planted in a complete randomized design at Nanbin Station (N18°21′7′′, E109°10′20′′), Sanya, Hainan province, China, in November 2013, in two-row plots of 5 m length and spaced 0.60 m apart, with 25 plants per row. Standard crop management practices were followed throughout the growing season to ensure normal growth without weed or insect damage. Natural infection with airborne uredospores of P. polysora was employed to induce the disease phenotype. To provide an abundant supply of the SCR pathogen, highly SCR-susceptible inbred lines Ye478 and Lx9801 were planted in a rectangular lattice across the field surrounding the SCR-phenotyping plots; this controlled field test plot design has been confirmed as an effective way to ensure efficient spread of the pathogen to neighboring plants in fields (Jines et al., 2007; Zhao et al., 2013) and has been adopted as the standard method of SCR phenotyping in maize breeding and genetic research in China (Zhao et al., 2013; Deng et al., 2020; Lv et al., 2020; Wang et al., 2020b). Disease scores were recorded at 30 days post-pollination. The degree of infection on each plant was visually scored on a scale of 1 to 9 using the Stakman infection type scale with 1 = HR (highly resistant, almost free of rust pustules), 3 = R (resistant, infection covering 6% to 10% of leaf area), 5 = MR (moderately resistant, infection covering 11% to 30% of leaf area), 7 = S (susceptible, infection covering 31% to 70% of leaf area), and 9 = HS (highly susceptible, nearly completely covered with rust pustules) (Wang et al., 2014; Deng et al., 2019; Deng et al., 2020; Li et al., 2021). The scores for each plot, based on the average of all plants in the plot, were used as the phenotype scores for each inbred line.

A leaf sample from each inbred line was used for DNA extraction via a CTAB procedure. The protocol reported by Elshire et al. (2011) was followed for DNA sequencing. Briefly, genomic DNA was digested with the restriction enzyme ApeK1, and genotyping-by-sequencing (GBS) libraries were then constructed in 96-plex and sequenced on an Illumina HiSeq 2000. SNP calling was performed using the TASSEL-GBS pipeline (Glaubitz et al., 2014), with B73 RefGen V2.0 as the reference genome. Initially, 87,6297 SNPs were filtered with minor allele frequency (MAF) > 5%, missing rate< 20%, and heterozygosity rate< 25% (Shu et al., 2021); data on 68,768 high-quality SNP loci were retained for entry into all analyses conducted in this study.

ADMIXTURE 1.3.0 (http://dalexander.github.io/admixture/download.html) was used to determine the population structure among all 140 inbred maize lines using 68,768 SNPs (Alexander and Lange, 2011), and the “Expectation Maximization” clustering algorithm was run with cluster numbers ranging from K = 1 to K = 12 in order to determine the optimum number of clusters K producing the smallest cross-validation error. Linkage disequilibrium (LD) analyses were carried out using TASSEL 5.2.25 (https://www.maizegenetics.net/tassel) (Bradbury et al., 2007; Glaubitz et al., 2014). The following six multi-locus GWAS methods from the R software package mrMLM.GUI V4.0.2 (Zhang et al., 2019; Zhang et al., 2020) were used to detect significant QTNs for SCR resistance: (1) mrMLM (https://cran.r-project.org/web/packages/mrMLM/index.html, Wang et al., 2016); (2) FASTmrEMMA (Wen et al., 2018); (3) ISIS EM BLASSO (Tamba et al., 2017); (4) pLARmEB (Zhang et al., 2017a); (5) pKWmEB (Ren et al., 2018); and (6) FASTmrMLM (Tamba & Zhang, 2018, bioRxiv). Q-matrices generated using ADMIXTURE 1.3.0 and K-matrices generated using mrMLM.GUI V4.0.2 were applied to correct the population structure and polygenic backgrounds. An intermediate result file and a final result file were generated. Manhattan plots were generated based on the intermediate and final results files using mrMLM.GUI V4.0.2.

Case–control sampling is a common design in human genetic mapping, in which cases and controls are defined based on the presence or absence of a certain phenotype, respectively (Uffelmann et al., 2021). In this study, case–control sampling was employed to verify the phenotype–SNP association detected in GWAS. Briefly, 10 inbred lines with the highest SCR disease phenotype scores (all with phenotype score 9) were sampled as cases, and 10 inbred lines with the lowest SCR disease phenotype scores (all with phenotype score 1) were sampled as controls; the SNP allele genotype and the mean phenotypic value for the corresponding allele for each QTN identified in multi-locus GWAS were compared between the case sample and the control sample for consistency. The phenotypes are visualized on a color scale in (with phenotype values from 1–9 expressed on a color scale from blue to red).

FASTA sequences containing significant QTNs were re-aligned to the B73 v4 reference genome in order to obtain a more accurate physical localization for improved gene annotations (https://www.maizegdb.org/gbrowse). For each statistically significant QTN identified via multi-locus GWAS, the annotated genes within ±300kb of the chromosome region (based on B73 Reference V4) were considered as candidate genes, as described in Ahmed et al. (2022). The location of a QTN in a gene and its effect were analyzed using the ANNOVAR software package (Wang et al., 2010). The polymorphic SNPs surrounding key significant QTNs, their SCR phenotype associations, and their association with gene structures were explored using scatter plots and gene structure diagrams. The expression profiles for each candidate gene in different organs and tissues at different developmental stages and under exposure to different biotic and abiotic stresses were extracted from the transcriptomic databases at MaizeGDB (https://www.maizegdb.org/). The relationship between expression of each candidate gene and resistance to southern corn rust was analyzed.

Using the control procedure of the Tassel pipline (see Materials and Methods section), the quality of the genotype data of the set of 140 inbred lines was improved by filtering out SNPs with minor allele frequency (MAF)< 5% and SNPs with more than 20% of the genotyping data missing ( Figure S1 ). The filtered genotype data of 68,769 SNP loci were evaluated for LD decay on all 10 maize chromosomes ( Figure 1A ); the mean LD decay distance across all chromosomes for this diversity panel was r² = 0.1, or approximately 200~300 kb ( Figure 1A ). The 140 inbred lines were evaluated for population structure using the SNP data (see Material and Methods section). ADMIXTURE 1.3.0 was used to detect population structure; the cross-validation error was smallest at K=6 ( Figure 1B ), indicating there were 6 subgroups. This population structure is visualized in Figure 1C .

Genome-wide association study (GWAS) was conducted using SCR phenotype data, SNP genotype data from the diversity panel of 140 inbred lines, and the population structure information through use of a multi-locus GWAS R package that offers six different statistical and computational models (see Material and Methods section). All statistically significant QTNs, including 21 QTNs identified by GWAS (defined as a LOD score greater than 3.00), are reported in Supplementary Table S1 ; among these, 13 QTN loci associated with putative candidate genes are listed in Table 1 . Pairwise LD values among 10 QTNs are shown in Supplementary Figure S2 . Among the 21 QTNs identified by GWAS ( Supplementary Table S1 ), eight QTNs were located on exons, six in intergenic regions, three in the upstream regions of genes, two on introns, one in the UTR3 region, and one in the UTR5 region. Among the eight QTNs located on exons, three variants located on Zm00001d032244, Zm00001d033259, and Zm00001d052781 were identified as nonsynonymous SNVs causing an amino acid change, whereas 5 variants were identified as synonymous SNVs not causing an amino acid change ( Table S1 ).

The map positions of the 13 QTNs reported in Table 1 were also labeled on a Manhattan plot ( Figure 2 ) produced using the R software package mrMLM.GUI V4.0.2 (see Materials and Methods section). As shown in Table 1 ; Supplementary Table S1 ; Figure 2 , the 13 significant QTNs were located on Chromosomes 1, 2, 4, 5, 6, and 8; some chromosomal regions were found to contain multiple tightly linked QTNs, such as bin1.11 (QTNs S1_299a and S1_299b (nbcs4)), bin5.06 (QTNs S5_210 (Fcf2) and S5_211(cel25)), and bin6.07 (QTN S6_164a (hsp101) and S6_164b (hsp101)). Three QTNs that were found to explain more than 10% of the phenotypic variation (r²%) are shown in Table 1 ; these are S1_299b (nbcs4), S5_145 (GTE4), and S5_211 (cel25).

As shown in Table 2 , the two alleles from each QTN were associated with significantly different phenotypic effects on their carriers. Based on measures of the difference in allelic effect, the 13 QTNs fell into four different categories or allele effect types: (1) type A, or HR/S type: one allele is associated with high resistance (HR) and the other with susceptibility (S); (2) type B, or MR/S type: the two alleles are associated with moderate resistance (MR) and susceptibility (S), respectively; (3) type C, or MR/HS type: the two alleles are associated with moderate resistance (MR) and high susceptibility (HS), respectively; and (4) type D, or S/HS type: the two alleles are associated with susceptibility (S) and high susceptibility (HS), respectively. From Table 2 , it can be seen that most of the highly resistant (HR) inbred lines, such as Qi319, X178, P138, CT01, and CT1251, carry not only HR alleles of type A QTNs (S5_145, S5_210), but also either MR alleles of type B and type C QTNs or S alleles of type D QTNs; in contrast, most of the highly susceptible (HS) inbred lines, such as Zong3, Zheng58, 478, and Xun926, carriers of the HS alleles of type D QTNs ( Tables 2 , 3 ).

The allele phenotypic effect of each single QTN ( Table 3 ) and the joint phenotypic effects of haplotypes from two QTNs (known as a QTN combo) were estimated via allele/haplotype effect analysis ( Table 3 ); these effects were also visualized using boxplots, two examples of which are shown in Figure 3 . QTN S6_164a and S6_164b ( Table 3 , combo 1) were two tightly linked QTN loci with D′= 0.93 ( Supplementary Figure S1 ), and physically were only 3036 bp apart on Chr6 (see B73 RefGen V4 position in Table 2 ); both were type D QTN loci (S/HS allele effect type), and both were found to have AA and GG alleles ( Figures 3 -A2, 3-A3). In order to evaluate their joint phenotype effects, the haplotype effects for all four observed haplotypes (GGGG, AAGG, GGAA, AAAA), in the form of the average SCR score of their carrier inbred lines, were plotted as shown in Figure 3 -A1. The haplotypes GGGG and AAAA were found have an almost identical SCR phenotypic effect to that of the GG and AA alleles of a single QTN alone ( Figures 3 -A2, 3-A3); the haplotype AAGG and GGAA had very few counts, indicating that these are likely recombinant types, and the GGGG and AAAA are likely parental types. The joint phenotypic effects of combo 2 (S6_165.S8_123) are shown in Table 3 and in Figure 3 -B1, 3-B2, 3-B3. S6_165 (MR/S) and S8_123 (MR/HS) were found to be two almost completely independent QTNs (D′= 0.28, Figure S1 ), producing four observed haplotypes (AATT, GGTT, AAGG, GGGG); as shown in Figure 3 -B1 and Table 3 , the haplotype AATT carriers were found to have an SCR phenotypic effect of 2.4, and were more likely to be resistant (R) inbred lines than either the AA allele carriers of QTN S6_165 (5.1, Figure 3 -B2; Table 3 ) or the TT allele carrier of S8_123 (4.4, Figure 3 -B3; Table 3 ), both of which were found to have an SCR phenotypic effect of around 5.0. Similarly, the GGGG haplotype was found to be a better predictor of an HS phenotype in its carriers than the GG allele of S6_165 or the GG allele of S8_123 alone ( Figure 3 -B1, 3-B2, 3-B3, Table 3 ). As shown in Table 3 , the maximum haplotype effect of combo 2 (S6_165.S8_123, AATT vs GGGG) was -6.0, much larger than the allele effect of QTN S8_123 (AA vs GG, -3.6) or QTN S6_165 (TT vs GG, -2.7) alone by absolute value. In comparing the two combos, it is clear that combo 2 has much better predictive value, because the two QTNs have low LD and they can freely shuffle to generate a more advantageous allele combination or haplotype than combo 1.

To verify the significant association between the SCR phenotype and QTN genotype that we had observed in multi-locus GWAS, a case–control sampling study was conducted (see Materials and Methods section). Specifically, 10 case inbred lines and 10 control inbred lines were genotyped for 10 QTN loci; the results are summarized in Figure 4 , where rows represent inbred carrier, columns represent QTNs, and the letters and the color of each cell C_ij represent the SNP genotype and the allelic phenotypic value of QTN j of inbred carrier i, respectively. The allelic effect type (A, B, C, or D, from Table 2 ) for each QTN are indicated in the top margin of Figure 4 . Based on Figure 4 ; Table 2 , our findings were as follows: (1) at type A QTNs, all case inbred lines carried S alleles in type A and type B QTNs and had a susceptible allelic effect phenotype (shown by the orange color), whereas all control inbred lines carried HR alleles and had a highly resistant allelic effect phenotype (shown by the blue color) at type A QTNs, with the exception of Xun248, which was found to carry only S alleles of type A, and CT019-21, for which genotype data were missing; (2) all case inbred lines carried HS alleles and had highly susceptible allelic effect phenotype (shown by the red color) at two type C QTNs and at one or more type D QTNs, while none of the control inbred lines carried HS alleles of type D QTNs. As shown in Figure 4 , the associations between the allelic genotypes carried by the case and control inbred lines and the phenotypes observed were highly consistent across different QTNs.

For each of the 13 QTNs and the 8 combos of two QTNs, the allelic (haplotype) effect, expressed as the difference between its resistant alleles (or haplotypes) and its susceptible alleles (or haplotypes) in terms of average SCR phenotype score is reported in Table 3 ; a negative value indicates that the resistant alleles (or haplotypes) can reduce the SCR score of their carriers, thereby enhancing the carrier’s disease resistance performance. Among single QTNs, the two type A QTNs had the largest allele effect by absolute value, with a value of -5.8. Type D QTNs had the smallest allelic effect by absolute value (ranging from -2.1 to -1.1; Table 3 ). Among combos of two QTNs, the combo of two type A QTNs (combo 3) had an effect of -6.5, larger than any single type A QTN, while the combo of one type D QTN and one type A QTN (combo 8, S2_12.S5.210) had the largest allele effect, -7.2. For each QTN, the relative performance of its resistant allele over its susceptible allele was measured in the form of RR_R, and the relative performance of its susceptible allele over its resistant allele was measured in the form of RR_S. As shown in Table 3 , for any QTN, its resistant allele significantly enhanced disease resistance performance over its susceptible allele, with RR_R values ranging from 2.89- to 9.65-fold, and its susceptible alleles significantly increased disease risk over its resistant allele, with RR_S values ranging from 1.25- to 17.64-fold. For combos of two QTNs, RR_R values ranged from 6.09- to 56.45-fold, significantly greater than values for a single QTN, and RR_S values ranged from 1.39- to 19.38-fold, only slightly greater than values for a single QTN. The above results suggest that the disease resistant haplotype arising from a combo of two QTNs could significantly increase resistance performance without significantly increasing the disease risk of its carriers.

According to the QTN–QTN interaction estimates calculated using Plink v1.07 (see Supplementary Table S2 ), among the 78 QTN–QTN pairs generated from the 13 QTNs shown in Table 3 (1/2(13×12)), six pairs showed very significant epistatic effects (p<0.01) and nine pairs showed significant epistatic effects (0.01≤p<0.05). Among the eight combos (QTN–QTN pairs) listed in Table 3 , three showed epistatic effects (0.01≤p<0.05), and the absolute SCR_R value was little greater than that of the more resistant QTN in the pair in each case, with differences of -0.1 (SCR_R of S5_211.S8_123 minus S8_123), 0.3 (SCR_R of S6_165.S6_164a minus S6_165), and 0.5 (SCR_R of S5_145.S5_210 minus S5_145 or S5_210). Combo 2 (S6_165.S8_123), with SCRR=2.4, showed no significant epistatic effect (p = 0.1273). These findings indicate that epistasis or gene modification is quite common in the gene regulation network, but other types of interaction, such as additive effects, likely also play important roles in enhancing SCR resistance.

The associations between resistance to southern corn rust (SCR) and the polymorphic SNPs surrounding key significant QTNs were further examined using the intermediate results file via a scatter diagram. The relationships between associations and gene structures were studied using a gene structure diagram ( Figure 5 ). The transcriptomic profiles of the candidate genes were obtained from maize transcriptomic databases via MaizeGDB (see Material and Methods section, Supplementary Figures S3A-L , S4A-L ). Supplementary Figure S2 shows the expression profiles of these candidate genes in different tissues, as obtained from MaizeGBD, and Supplementary Figure S3 shows an organ-specific and stress-induced gene expression atlas for them, also retrieved from MaizeGDB. Here, we only list a select set of QTNs and putative candidate genes for further scrutiny.

(1) QTN S1_218 ( Table 1 , V4: chr1:218619398, bin 1.07) was identified as significant by two GWAS models, with an LOD score of 4.58 and phenotype contribution of 8.1%. S1_218 was located in the CDS region of Zm00001d032244 (probable inactive receptor kinase) and A MYB-transcription factor 146 (Zm00001d032240) was located 110kb upstream; Zm00001d032240 is notably upregulated after infection by Cercospora zeina ( Supplementary Figure S4A ) (Hoopes et al., 2019). MYB transcription factors have been reported to function as negative regulators of plant immunity (Wang et al., 2021) and positive regulators of the defensive response to stripe rust in wheat (Zhu et al., 2021), and have also been found to be associated with resistance to rust disease in foxtail millet (Bai et al., 2019).

(2) QTN S1_299b ( Figure 5A , V4: chr1:299623487, bin 1.11): with an LOD score of 8.58 and phenotype contribution of 18.13%, this encodes a permease I (Zm00001d034678, nbcs4, nucleobase-cation symporter4). According to the results of genome scanning of the genome region of Zm00001d034678 (nbcs4) ( Figure 5A ), three SNPs (Chr1: 299623551, 299623552, 299623554) with -log₁₀ (Pvalue.Q)< 2 were located on the 5’ UTR region of Zm00001d034678 (nbcs4); a high level of gene expression in leaves was also reported ( Supplementary Figure S3B ). One of the nucleobase-cation symporters (NBCS), named MdNAT7, has been reported to transport xanthine and uric acid, which enable more efficient scavenging of ROS by enhancing H₂O₂-scavenging enzymes in apple (Sun et al., 2021a).

(3) QTN S2_12 ( Table 1 , V4: chr2:12916090, bin 2.01) was identified by two GWAS models, with an LOD score of 5.07 and phenotype contribution of 4.06%; an rlk12 gene (Zm00001d002447, receptor-like protein kinase 12) was located 11kb upstream, and has been found to be a fungal resistance-related gene (Zhang et al., 2017b; Yang et al., 2019), with Zm00001d002447 having been found to express in mature leaves and its level of expression to be elevated by fungal infection ( Supplementary Figure S3C ; Supplementary Figure S4C ) (Hoopes et al., 2019).

(4) QTN S5_145 ( Figure 5B ,V4: chr5:145816043, bin 5.04): with an LOD score of 3.47 and phenotype contribution of 11.16%, two SNPs (Chr5:145816043, 145817495) with -log₁₀ (Pvalue.Q) > 2 were located on the CDS regions (the SNP chr5:145816043 is a synonymous SNV, and the SNP chr5:145817495 is a nonsynonymous SNV; Supplementary Table S3 ) of Zm00001d016131 (Transcription factor GTE4); Zm00001d016131 has been found to express in leaves, and its expression is elevated by fungal infection ( Supplementary Figure S3F , Figure S4F ) (Hoopes et al., 2019).

(5) QTN S5_211 ( Figure 5C , V4: chr5:211183684, bin 5.06) was identified by three GWAS models, with an LOD score 6.81 and phenotype contribution 11.69%. S5_211 was located in the CDS region (synonymous SNV, Supplementary Table S3 ) of Zm00001d017978 (cellulase25). An immediate-to-early fungal elicitor protein (CMPG1) encoding gene, Zm00001d017976, was located 280kb upstream; expression of Zm00001d017976 was significantly elevated by Colletotrichum graminicola infection ( Supplementary Figure S4H ) (Hoopes et al., 2019).

(6) QTN S6_164a (V4: chr6:164808768, bin 6.07) and QTN S6_164b (V4: chr6:164811804, bin 6.07) were two QTNs found to be only 3036bp apart and tightly linked with high LD ( Supplementary Figure S1 ). The two significant QTNs S6_164a (Chr6:164808768) and S6_164b (Chr6:164811804) were located in the CDS (synonymous SNV, Supplementary Table S3 ) and 3′UTR regions of Zm00001d038806 (hsp101), respectively ( Figure 5D ); the expression of Zm00001d038806 (hsp101) is notably upregulated when the temperature increases ( Figure S4K ) (Hoopes et al., 2019). Another gene named rlk10 (Zm00001d038791, receptor-like protein kinase 10) was located 410kb upstream; the expression of Zm00001d038791(rlk10) has been found to be significantly elevated by Colletotrichum graminicola infection ( Supplementary Figure S4J ) (Hoopes et al., 2019). Overall, Zm00001d038806 (hsp101) and Zm00001d038791(rlk10) are two potential candidate genes for these two QTNs.

(7) QTN S8_123 (V4:chr8:123503579, bin 8.06) was identified by two GWAS models; two genes encoding Pgk2- phosphoglycerate kinase 2 (also known as Metacaspase type II, Zm00001d010672, Zm00001d010673) were located 300kb downstream and might be related to SCR resistance (Luan et al., 2020). Zm00001d010672 has been found to be expressed in mature leaves and its expression is elevated by fungal infection ( Supplementary Figures S3K, L ; Supplementary Figure S4L ) (Hoopes et al., 2019).

Genetic mapping, either by linkage-based mapping or association-based mapping, mainly provides information on associations between DNA markers and trait phenotypic variation; further verification is required to establish that the association represents a causal relationship between underlying genetic variants and phenotypic variation. Traditional QTL mapping using linkage-based F₂-derived populations is limited by a low frequency of chromosome recombination and low density of molecular markers; the QTLs identified are far away from the causal DNA variants, and intensive fine mapping is required to further pinpoint the precise chromosome location. A large number of SCR-related genetic loci have been reported since the 1960s, including some major effect QTLs on chromosome 10S: RppM (Wang et al., 2020a; Wang et al., 2022), RppS/RppK (Wu et al., 2015; Chen et al., 2022), RppQ (Chen et al., 2004), RppC (Deng et al., 2022), RppP25 (Liu et al., 2003), RppD (Zhang et al., 2009), and RppCML496 (Lv et al., 2020). Only a small number of QTLs, such as RppC, RppK, and RppM, have been verified by further genetic fine mapping and molecular studies (Chen et al., 2022; Deng et al., 2022; Wang et al., 2022). A set of genes/QTLs including ZmREM1.3 (Wang et al., 2019b), qSCR4.01 (Deng et al., 2020), qSCR4.08 (Deng et al., 2019), qSCR6.01(Lu et al., 2020), ZmMM1(Wang et al., 2021), and qSCR9.04 (Deng et al., 2019; Deng et al., 2020; Li et al., 2023b) have also been reported on chromosomes 2, 4, 6, 7, 8, and 9, but most of these have been identified based on preliminary linkage-based mapping data or comparative proteomic analysis. Association-based mapping strategies, such as GWAS, with a large mapping population and high-density SNP markers, have been shown to be very powerful in locating causal chromosome variation and precise DNA variants without further fine mapping, but are also more likely to generate a large number of false positive signals, and thus require extensive post-GWAS verification (Wen et al., 2018; Xiao et al., 2017; Wang et al., 2019a; Zhang et al., 2020; Li and Ritchie, 2021; Shikha et al., 2021). GWAS has been employed to identify SCR-related QTNs in several studies (Zhou et al., 2017; De Souza Camacho et al., 2019; Sun et al., 2022; Li et al., 2023b). The QTNs Chr4:173,863,109 and Chr6:169,030,253, identified by Sun et al. (2022), are approximately 4Mb away from the QTNs reported in this GWAS study: QTNs S4_170 (Chr4:170324384) and S6_164a (Chr6:164808768), and QTNs S6_164b(Chr6:164811804) and S6_165 (Chr6:165682422), respectively. Zm00001d029980 (mca1 - metacaspase 1, Chr1:97331766-97337093) and Zm00001d047078 (mca2 - metacaspase 2, Chr9:117546061-117551091) are genes related to disease resistance in the NLR (Rp1-D21)-mediated defense response (Luan et al., 2020). The NLR protein Rp1-D21 is coded by Zm00001d023317 (rp1-resistance to Puccinia sorghi 1, Chr10:2861471-2865816). Zm00001d010672 and Zm00001d010673, as candidate genes of QTN S8_123 in our study, encode phosphoglycerate kinase 2, which is also known as metacaspase.

Most of the QTNs identified on chromosomes 1, 2, 4, 5, 6, and 8 were found to be type B, type C, and type D QTNs with small effects. To minimize the risk of reporting false positive QTNs, six different multi-locus GWAS models from Zhang et al. (2020) were applied to generate a shortlist of significant QTNs; two association-based post-GWAS procedures (case–control sampling and allele/haplotype effect analysis) were then followed for verification. For ultimate confirmation that a statistical association represents a true causal relationship between the QTNs identified and the SCR-related phenotypes observed, additional molecular genetic tools, such as transgenic studies and genome-editing, or signal transduction cascade verification, are needed (Zhu et al., 2021; Chen et al., 2022; Deng et al., 2022; Wang et al., 2022).

One critical component for successful application of GWAS is genotyping. DNA markers with sufficient density and low LD are important for association mapping. Due to multiple occurrences of historical chromosomal recombination, the LD block in diversity panel of inbred lines is expected to be much shorter than in linage-based F₂-derived populations; low LD together with high marker density render GWAS with a diversity panel an effective approach to genetic mapping. In the past, the trade-off between genotyping cost and marker density had always been an issue: although SSR and SNP-chip platforms are more cost-effective than complete sequencing methods such as NGS, they yield fewer markers and a lower marker density on chromosomes, and thus have lower mapping resolution. In this study, we solved this dilemma by using an NGS platform known as Genotyping by Sequencing (GBS), which uses the DNA methylation-sensitive enzyme ApeK1 to cut maize chromosomes. Because ApeK1 cuts more frequently at less methylated, gene-rich regions of the maize chromosomes, it generates more short DNA fragments and thus more DNA sequencing reads and more SNP calls. In this project, even though the GBS method only achieves 5X genome coverage on average, the real genome coverage can reach over 30X in gene-rich chromosomal regions and may yield over 700K SNP loci (Bradbury et al., 2007; Elshire et al., 2011; Glaubitz et al., 2014; Shu et al., 2021).

The selection of appropriate statistical models to detect and measure association is critical to the success of GWAS. The models should be able to deal with various features of phenotypic and genotype data, such as continuity and normality of phenotypic data, population structure and kinship information in genotype data, and confounding from other covariables. The R software package provided by Zhang’s group, mrMLM.GUI V4.0.2 (Zhang et al., 2020), includes six multi-locus GWAS statistical models and a multiple-step algorithm for dealing with different data types and reducing false positives. Under the framework of multi-locus random-SNP-effect mixed linear modeling, each marker on the maize chromosome was first scanned for statistical significance and a less stringent Bonforroni correction was adopted in the statistical test. The significant marker loci identified were then incorporated into a new multi-locus genetic model; their effects were estimated by an empirical Bayes method, and all non-zero effects were further evaluated by the likelihood ratio test. Comparative studies have shown that this R package is more powerful than other methods both in terms of the detection of both strong and weak signals and in terms of lower output of false-positive signals (Wang et al., 2016; Tamba et al., 2017; Wen et al., 2018; Zhang et al., 2017a; Ren et al., 2018; Tamba and Zhang, 2018; Wang et al., 2019a; Zhang et al., 2020; Uffelmann et al., 2021; Park et al., 2022; Li et al., 2023a). As shown in Table 1 ; Supplementary Table S1 , a large portion of QTNs were detected by more than one models, indicating that these QTNs have strong associations with the SCR trait.

As reported in other GWAS projects in human and maize, our GWAS also identified a large number of statistically significant SNP variants or QTNs ( Table 1 ; Supplementary Table S1 ), and the majority of the QTNs were located in non-coding regions of a putative gene or maize genome regions of unknown biological function. Verification of QTNs identified by GWAS has always represented a major challenge, but great progress has recently been made, with various post-GWAS methods of analysis having been developed and some important genes for complex human diseases having been successfully identified (Li and Ritchie, 2021; Uffelmann et al., 2021). In this study, two association-based post-GWAS procedures were conducted: first, case–control sampling was performed to verify the QTNs identified by GWAS and to check the expected QTN genotype predicted by the phenotype for consistency; and second, the SNPs near the QTNs were collected, their haplotypes were generated, and haplotype–trait associations were examined (see Figures 3 – 5 ). Additionally, putative genes within a 300kb neighborhood of the QTN were taken as candidate genes, and their transcriptomic profiles and biological functions were examined. The genome-wide transcriptomic association study (TWAS) approach has gained in popularity over the years due to its distinct ability to connect SNP variants with expressed genes and complex trait phenotypes (Li and Ritchie, 2021; Uffelmann et al., 2021). However, the majority of QTNs identified in this study were located in non-coding regions of putative genes or non-genic regions of chromosomes; similar findings have also been reported in a number of studies of humans, rice, and maize, where the majority of trait-associated DNA markers are located in non-genic regions, and their genetic effects are generally greater than those in genic regions (Li et al., 2012; Rodgers-Melnick et al., 2016; Wei et al., 2020). Many studies have shown that complex cis- and trans-transcriptional regulatory mechanisms do not take effect within the coding regions of genes and more likely reside in the non-coding genomic neighborhood of identified QTNs (Li and Ritchie, 2021); thus additional approaches for post-GWAS analysis, such as transgenic validation, genome editing, and molecular function analysis, are necessary.

In this study, we attempted to connect the QTNs identified by GWAS with molecular breeding practice. Firstly, the allele of each QTN or a haplotype of a QTN combo was assigned to one of the five categories HR, R, MR, S, or HS based on the value of its allelic (or haplotype) effect; subsequently, thirteen QTNs were grouped into four allelic effect types (A, B, C, and D) according to their combination of allele types. The two type A QTNs, S5_145 and S5_210, both located on chromosome 5, were found to have large allelic effects; each contributed approximately 10% of SCR-related phenotypic variation (r^2%), and their highly resistant (HR) alleles were present in all major SCR-resistant inbred lines in the diversity panel. The two type B, two type C, and 7 type D loci, located on chromosomes 1, 2, 4, 5, 6, and 8, were found to have small allelic effects, and likely provide a broad spectrum of SCR resistance to resistant inbred lines. Following this analysis, the phenotypic impact of each QTN (and its alleles) on its carrier was measured in terms of two statistics, relative resistance (RR_R) and relative susceptibility (RR_S; see Materials and Methods for definitions); these measure the resistance efficacy and disease risk conferred on the carrier by a resistant allele or a susceptible allele, enabling maize breeders to easily understand the benefits and risks of integrating a particular allele or QTN into their breeding materials or into new maize varieties. GWAS has been shown to be effective in identifying QTLs with small effects or conferring broad resistance (Xu et al., 2014; Zhou et al., 2017; Xu et al., 2018; De Souza Camacho et al., 2019; Zhang et al., 2019; Sun et al., 2022). Most published SCR-related QTL-mapping projects have used F₂-derived segregated populations or their permanent versions, such as DHs and RILs from a genetic cross of a highly resistant (HR) parent × a highly susceptible (HS) parent; such studies have identified a number of QTLs associated with major effects on the short arm of chromosome 10, but very few small-effect QTLs have been reported and verified (Zhao et al., 2013; Lv et al., 2020; Wang et al., 2020a), indicating that linkage-based QTL-mapping using the HR × HS genetic cross method is less effective in identifying small-effect resistance genetic loci for broad resistance. Large-scale GWAS using a large, well-designed diversity panel of breeding lines and high-density SNP markers is likely to be a more effective way of discovering small-effect, broad SCR-resistance genetic loci.