The data for the 26 maize pan-genomes were obtained from the study by Hufford et al. (2021). Protein domains were identified using HMMER v3.1 with an E-value threshold of 1e−15 to search against the following Pfam profiles: PF00931 (NB-ARC), PF01582 (TIR), PF05659 (RPW8), and PF00560/PF07725/PF12799 (LRR) (Ma et al., 2021; Sun et al., 2014; Wang et al., 2020). In addition, the presence of CC domains was examined using EMBOSS v6.6.0. The integrity of the identified domains was further validated using the Conserved Domain Database (CDD) search tool provided by NCBI (Wang et al., 2023). Custom scripts were used to extract and calculate molecular weight, isoelectric point, and other protein characteristics (Supplementary Table S1).

The presence–absence variation information of ZmNBSs was obtained from (Hufford et al., 2021). The ggplot2 package in R was used to construct the PAV heatmap.

We performed correlation analysis between the PAV frequency (PAV frequency = number of lines containing the gene/total number of lines) and Ka/Ks ratios of ZmNBS gene subfamilies. We conducted the analysis using the cor.test function in the R software (version 4.3.2), employing Pearson’s correlation coefficient. Prior to correlation analysis, we filtered extreme Ka/Ks values (Ka/Ks > 5) that may represent computational artifacts or recently duplicated genes undergoing strong positive selection, as such outliers can disproportionately influence correlation estimates.

The classification of the ZmNBS subfamily is based on the research by Hufford et al. (2021). The conserved domains of ZmNBS subfamilies were used for phylogenetic tree construction. To generate a robust and interpretable phylogenetic tree from the large set of 1,457 genes, a representative sequence from each of the 129 identified subfamilies was selected based on the integrity and conservation of the NBS domain. Multiple sequence alignment of the conserved domains was performed using MUSCLE, followed by model selection using IQ-TREE v1.6.9 to determine the best-fit evolutionary model (Nguyen et al., 2015). A maximum likelihood (ML) phylogenetic tree was then constructed with 1,000 bootstrap replicates. The resulting tree was visualized and refined using TVBOT (Xie et al., 2023).

To investigate the evolutionary dynamics of ZmNBS genes, two approaches were employed to calculate Ka/Ks values.

The protein sequences of sorghum (Sorghum bicolor) were used as an outgroup reference (https://sorghum.genetics.ac.cn/SGMD/Genome.html) (Supplementary Table S2). Orthologous gene pairs between sorghum and maize were identified through protein sequence alignment, and the corresponding Coding sequence (CDS) pairs were used to calculate Ka and Ks values using the KaKs_Calculator. The replication type of ZmNBS was analyzed based on DupGen-finder (Supplementary Table S2).

Based on the ZmNBS protein and CDS sequences obtained from the 26 maize genomes reported by Hufford et al. (2021), pairwise comparisons among ZmNBS gene copies were performed (Supplementary Table S3). After aligning coding sequences, Ka and Ks values were estimated using the KaKs_Calculator to assess selection pressure among intraspecific paralogs (Wang et al., 2010). The pheatmap and ggplot2 packages were used for subsequent plotting.

Structural variation (SV) data were obtained from the study conducted by Hufford et al. (2021) using the B73 genome as the reference for SV identification. Structural variants were annotated using ANNOVAR (Wang et al., 2010), and those associated with ZmNBS gene family members were subsequently extracted. Correlation analysis was then performed between structural variation and gene expression levels of ZmNBS genes. For genes showing statistically significant expression differences associated with SV presence, bar plots were generated to visualize their expression levels. Subsequently, the maize accession with the highest number of ZmNBS genes overlapping structural variation regions was selected for comparative analysis with the reference genome. MEME Suite v5.5.8 (https://meme-suite.org/meme/tools/meme) was employed to analyze and compare conserved motifs between the ZmNBS proteins of this accession and those from the reference genome (Bailey et al., 2009). Subsequently, the InterProScan 5.76-107.0 program was used to perform functional domain annotation for all identified motifs, with databases including Pfam, SMART, CDD, and PRINTS. The functional annotations obtained through this process were used for the subsequent analysis of how structural variations affect protein function.

RNA-seq data of maize B73 leaves subjected to Spodoptera litura feeding for 0, 6, and 12 h (Project ID: PRJCA003103) were used to analyze the expression patterns of ZmNBS genes. Low-quality reads were filtered using fastp (Chen et al., 2018), and the clean reads were aligned to the B73 reference genome using HISAT2 (Kim et al., 2019). Gene-level read counts were obtained using HTSeq (Putri et al., 2022), and differentially expressed genes (DEGs) were identified with the criteria of |log₂FoldChange| > 1 and False discovery rate (FDR) ≤ 0.05. Heatmaps of ZmNBS gene expression were visualized using the ComplexHeatmap package (v2.6.2) in R (Gu et al., 2016). Finally, the annotation files of the differentially expressed genes were extracted using custom scripts, and gene structure diagrams were generated using GSDS v2.0 (http://gsds.gao-lab.org/) (Hu et al., 2015).

In total, we identified 1,457 ZmNBS genes across the maize pan-genome, including 53 genes from the reference genome. Among the diverse inbred lines, B97, M37W, and Oh43 harbored notably more ZmNBS genes than other genotypes (Figure 1A). Using domain annotation tools such as Pfam, we classified all ZmNBS genes into four structural types (Table 1): CN, CNL, N, and NL. Consistent with previous findings, we identified no TIR or RPW8 domain-containing ZmNBS genes in maize (Lan et al., 2019; Li et al., 2015). Notably, we found that 775 genes contained the CC motif, which is significantly more than the 491 genes carrying LRR motifs, highlighting the overrepresentation of CC-type ZmNBS proteins in maize. A phylogenetic analysis of 129 full-length ZmNBS protein sequences revealed that N and CN types tended to cluster together within the same clade. Interestingly, we also found CNL-type members nested within this group, suggesting potential structural convergence or recent domain fusion events (Figure 1B).

We explored PAV patterns of ZmNBS genes across diverse maize inbred lines and observed that a significant portion of ZmNBS loci were variably present or entirely missing in certain genotypes (Figure 2), reflecting the high plasticity and genomic dynamism of the ZmNBS gene family.

Notably, several subfamilies—including ZmNBS11, ZmNBS13-14, ZmNBS17-19, and ZmNBS21-42—showed strong conservation and were retained across the majority of genotypes, suggesting that these genes may perform fundamental basal immune functions. In contrast, subfamilies such as ZmNBS1-10, ZmNBS12, ZmNBS15-16, and ZmNBS43–60 displayed strong genotype-specific absences, with some genes restricted to only one or two inbred lines. This implies that these variable genes may be involved in lineage-specific adaptive responses.

Additionally, considerable differences in overall ZmNBS gene retention were observed across genotypes. For instance, CML247, Ki3, B97, and HP301 harbored substantially more ZmNBS copies, whereas M162W and Tzi8 retained fewer ZmNBS genes. These patterns may reflect distinct selection histories or genetic bottlenecks associated with specific breeding lineages (Catlin et al., 2025).

To quantitatively assess the relationship between gene conservation and evolutionary constraint, we calculated the PAV frequency for each ZmNBS subfamily and correlated it with the mean Ka/Ks ratio. After filtering out extreme Ka/Ks values (Ka/Ks > 5) that may represent computational artifacts or recent duplications under strong positive selection, we observed a negative correlation between PAV frequency and Ka/Ks (r = −0.151, p = 0.189) (Supplementary Figure S1). Although not statistically significant, this trend suggests that core ZmNBS subfamilies (high PAV frequency) tend to evolve under stronger purifying selection (lower Ka/Ks), while variable subfamilies (low PAV frequency) may experience relaxed selection or positive selection (higher Ka/Ks). This pattern is consistent with the core-adaptive model of resistance gene evolution.

To explore the evolutionary constraints of the ZmNBS gene family, we first identified the duplication modes of ZmNBS genes (1,457) in all the maize lineage using the sorghum genome as an outgroup (Figure 3A). The statistical analysis of duplication modes indicates that DSD was the dominant mechanism for the expansion of the ZmNBS gene family, driving the formation of significantly more genes (755) than any other mode. TD, PD, and WGD accounted for 256, 162, and 134 genes, respectively. The analysis revealed a significant association between the structural classification of ZmNBS subfamilies and their duplication modes. Specifically, CNL- and NL-type genes are mainly derived from DSD; in contrast, the expansion of N-type genes is primarily driven by TD and DSD, indicating that these mechanisms play a more prominent role in the proliferation of non-canonical ZmNBS gene structures.

Subsequently, we evaluated the long-term selective pressures acting on ZmNBS genes in the maize lineage by calculating the Ka/Ks ratios between maize ZmNBS genes and their sorghum orthologs. The results showed significant differences in selective pressures among genes with different structural types (Figure 3B) and different duplication modes (Figure 3C).

At the gene structural level, the orthologous Ka/Ks ratios of CN-type genes were significantly higher than those of CNL, N, and NL types, indicating that CN-type genes may have experienced stronger adaptive selection during evolution. At the duplication mode level, compared with genes generated via WGD, genes derived from TD and PD had significantly higher Ka/Ks ratios with their sorghum orthologs. This suggests that genes of TD/PD origin may have an average faster evolutionary rate throughout their evolutionary history since divergence from the common ancestor and have been subjected to positive selection or relatively relaxed functional constraints. In contrast, genes derived from WGD exhibited the lowest Ka/Ks ratios, strongly implying that they have undergone continuous purifying selection throughout evolution, thereby being more likely to retain important core biological functions.

Taken together, these results collectively support the hypothesis that “duplication mode shapes structure, which in turn influences function”. Our analysis reveals an association between two evolutionary processes operating on different timescales: first, the specific duplication mode that a gene undergoes in the early stages of the maize lineage predetermines the domain composition of its encoded protein (e.g., CNL and N types). Subsequently, this inherent structural feature shapes the selective pressure the gene experiences during subsequent long-term evolution (since the divergence of maize and sorghum)—a pattern reflected in the Ka/Ks ratios between the gene and its sorghum orthologs. Therefore, the duplication origin of a gene, by influencing its structure, ultimately shapes the evolutionary trajectory it exhibits in cross-species comparisons.

We systematically identified the types of SVs in ZmNBS genes and their potential impacts. We detected multiple SVs in ZmNBS genes, primarily deletions and insertions. Most of these SVs are located in intergenic regions and may indirectly affect gene expression through regulatory sequences. A more critical finding, however, is that we identified SVs that can directly disrupt the coding structure of ZmNBS genes. For example, an insertion was found in the exon region of gene Zm00001eb015450; this variation is predicted to directly disrupt protein integrity and function, providing evidence that SVs directly contribute to the functional evolution of ZmNBS genes (Supplementary Table S4).

Furthermore, we analyzed the impact of SVs on gene expression. By comparing the expression profiles of genes with and without nearby SVs, we investigated whether structural variations affect gene expression. Several ZmNBS subfamilies, such as ZmNBS23 and ZmNBS24, showed markedly reduced expression in SV-positive lines (Figure 4A). Meanwhile, Figure 4C displays the differences in motif composition between SV-affected ZmNBS genes and reference ZmNBS genes, further supporting that SVs may impair functional expression by disrupting core domains. We next investigated the potential molecular mechanisms by which SVs alter gene function beyond expression changes. Functional annotation of SV-associated motifs revealed their direct impact on key protein domains. For instance, Motif5 (IPR001024) and Motif7 (IPR058922), present in CML52 but absent in B73, are predicted to alter subcellular localization for membrane-associated pathogen sensing and facilitate nuclear translocation for defense gene regulation, respectively. Conversely, the B73-specific Motif4 (IPR038005), which is absent in CML52, is a hallmark of resistance proteins like Rx and is critical for oligomerization and immune signalosome assembly. Therefore, these are not neutral variations but represent gains and losses of functional modules that directly drive functional divergence among paralogs, underscoring that SVs are a key mechanism in the adaptive evolution of ZmNBS genes.

Among different ZmNBS subtypes (CN, CNL, N, and NL), TD and segmental duplication (SL) were more enriched in N and NL types, while CNL and CN types were predominantly expanded through DSD, indicating a preferential association between duplication mechanisms and ZmNBS structural subtypes (Figure 4B).

To investigate stress-responsive expression, we conducted a time-course analysis following insect herbivory simulation. While many ZmNBS genes, ZmNBS26 and ZmNBS23, were transcriptionally upregulated at 6 or 12 h post-treatment, ZmNBS31 stood out for its high expression even at 0 h (Figure 5A). This pattern was consistent across biological replicates, suggesting that ZmNBS31 maintains a high basal level of expression and may function in early immune surveillance, with its regulation potentially involving a rapid feedback mechanism upon pathogen perception (Figure 5A).

Subsequently, we examined the Ka/Ks values of the ZmNBS31 subfamily across different maize lines. We found that most copies exhibited low Ka/Ks values (average < 0.35), suggesting that this subfamily has been subjected to long-term purifying selection at the population level and retains a high degree of functional conservation (Supplementary Figure S2). This finding is consistent with its high expression under pathogen treatment, indicating that ZmNBS31 members may have been preferentially retained during evolution to maintain basal immune function (Figure 5B). To further investigate the selective pressure acting on ZmNBS31 family members across different maize genomes, we calculated Ka/Ks values for each homolog (Supplementary Figure S2). Notably, the lowest Ka/Ks was observed in CML333 (0.2033), indicating that this allele has undergone particularly strong purifying selection. All ZmNBS31 copies exhibited Ka/Ks < 1, indicating that this gene family has been subject to strong purifying selection and maintains a high degree of functional conservation during evolution (Supplementary Figure S2).

Structural comparison across ZmNBS subfamilies (Figure 5C) revealed that genes derived from WGD and DSD generally possess intact CDS regions and canonical exon–intron structures, whereas PD-derived ZmNBS2 shows signs of CDS truncation and compressed structure, suggesting possible loss of function due to incomplete architecture. In contrast, ZmNBS31, which belongs to the DSD type, exhibits a more complex gene structure with multiple exons and extended introns, which may contribute to its high basal expression level, a configuration that aligns with its role in immediate immune surveillance and subsequent feedback regulation.