Two methods were used to identify maize membrane-trafficking genes. First, the maize genome and protein databases were downloaded from the Ensembl Plants database¹ (B73 RefGen version 3 and version 4). The protein sequences of Arabidopsis membrane-trafficking components were retrieved from the TAIR10 database² and employed as queries to search against local databases using BLASTP and TBLASTN methods (E value ≤ 1.0e-30). The orthology relationships between maize and Arabidopsis genes were further confirmed using the orthologs identification tool in the MCENet server³ (Tian et al., 2018). Second, the hidden Markov model (HMM) profiles of available membrane-trafficking components were retrieved from Pfam⁴ and used to search against the local maize protein database using the BLASTP method. The redundant sequences obtained from both methods were manually removed.

Full-length protein sequences of maize and Arabidopsis membrane-trafficking proteins were aligned using ClustalW with default parameters, and phylogenetic trees were generated by MEGA (version 6.0) software using the neighbor-joining method with 1,000 bootstrap replications. The resulting Newick trees were visualized with the Evolview-v2 online server.⁵ All the conserved domains of membrane-trafficking proteins were identified using the NCBI Batch CD-Search database⁶ (Lu et al., 2020) and Pfam website with default parameters. The conserved domains were visualized with TBtools (Chen et al., 2020a).

An R package, WGCNA, was used to construct gene co-expression networks (Langfelder and Horvath, 2008). For network analysis, two transcriptomic datasets, B73 different tissue RNA-seq data and Illumina HiSeq data (GEO accession: GSE92448), were collected and filtered to remove lower-quality reads and reads with expression variance of <20% (Stelpflug et al., 2016; Kebede et al., 2018). In the weighted gene network analysis, each node represented a gene, and the edge was determined by the similarity of paired genes. Pearson correlation coefficients were calculated based on the similarity of paired gene expression and then transformed into an adjacency matrix of connection strength, which was used to calculate the topographical overlap matrix (TOM). A 1-TOM was calculated using hierarchical clustering, which defined gene modules corresponding to the basis of dissimilarity of gene connectivity. The total connectivity and intramodular connectivity, KME (for eigengene-based connectivity), and the KME P-value were calculated. The eigengene value for each module (ME) with phenotypic traits was calculated, and the most significant module was considered the core module. Each module had a unique color. Gene connectivity in modules balanced the number of edges of a node and was further visualized in Gephi-0.9.2 software.⁷ Gene Ontology (GO) categories of co-expression interactions with membrane-trafficking genes in the module were performed using R 3.6.4 packages clusterProfiler (Yu et al., 2012).

To analyze the spatiotemporal expression pattern of maize membrane-trafficking genes, the expression data were retrieved from the previously described RNA-seq dataset (Stelpflug et al., 2016). To analyze the expression level of ZmSec23a in the seeds of regular and high oil inbred lines, the expression data in 15 DAP kernels of 368 maize inbred lines were downloaded from Yan’s lab website⁸ (Fu et al., 2013). The expression heatmaps were generated using the R.3.6.4 platform.⁹

The single nucleotide polymorphism (SNP) data were downloaded from the HapMapV3 study for 1,210 maize lines (HapMapV3.2.1)¹⁰ (Bukowski et al., 2018). After filtering out SNPs with ≤50% missing data, 28,803,913 SNPs were obtained. Moreover, 3.6 million SNPs were retrieved from the 368 RNA-seq association mapping panel, which include 32 high oil lines¹¹ (Fu et al., 2013). After filtering with a minor allele frequency of ≥5% using the TASSEL software (Bradbury et al., 2007), 557,977 higher-quality SNPs were obtained.

The identification of genomic selective sweep regions and selected vesicle membrane-trafficking genes during maize domestication and improvement was conducted as described in the previous report (Lin et al., 2014). Briefly, the nucleotide diversity (π) ratio and the fixation index (Fst) during domestication (teosinte vs. landrace) and improvement (landrace vs. improved) processes were calculated using VCFtools (Danecek et al., 2011), with a window size of 10 kb and step size of 2 kb. The top 5% of π ratio (1.55 and 1.22 for domestication and improvement, respectively) and the top 5% Fst values (0.25 and 0.089 for domestication and improvement, respectively) were set as the thresholds for putative selection signals. When the values of π or Fst of the genomic region exceeded the threshold, the detected membrane-trafficking genes in the region were defined as “domestication” or “improvement.”

Single nucleotide polymorphisms (MAFs ≥ 0.05) and full-length DNA sequences of candidate genes were extracted from wild teosintes, landraces, improved cultivars, and high oil accessions. These SNPs were used to calculate the square values of the correlation coefficient (r²) between each SNP pair to measure the linkage disequilibrium (LD) using HaploView version 4.2 (Barrett et al., 2005). Haplotypes were generated using DnaSP version 5 (Librado and Rozas, 2009). For haplotype network analysis of ZmSec23a (Zm00001d009150), a sample of 20 teosintes, 22 landraces, 58 maize inbred lines, 23 high oil materials, and 345 regular inbred lines was used to construct the haplotype network using POPART version 1.7 (Leigh and Bryant, 2015).

Putative promoter regions containing 6,000 bp of 5′ genomic sequences immediately upstream of the start code of the ZmVPS37A gene were extracted from the teosinte and maize genome database¹² and were submitted to PlantCARE¹³ for cis-element analysis.

Teosinte (Zea mays ssp. Parviglumis, accessions: W71-2), landrace (303WX), and B73 inbred lines were grown in a greenhouse at 25°C under a light/dark 16/8-h photoperiod. The Fusarium graminearum wild-type strain PH-1 used in this study was cultured at 25°C on a carnation leaf agar for 1 week and then transferred and incubated in liquid carboxymethyl cellulose (CMC) medium with gentle shaking (100 rpm) for three more days. The conidium was subsequently collected by filtering through the two layers of Miracloth (Millipore, 475855) and the number of conidia was counted using a hemocytometer. For the inoculation assay, in the stalk of maize V2 stage seedlings, a hole was created 2 cm above the soil using a sterile toothpick, 40 μl suspension conidia (1 × 10⁶ conidia/ml) or ddH₂O (mock) was injected into the hole with a needleless syringe, and then it was placed in a greenhouse for 7 days (70% humidity, 25°C, a light/dark 16/8-h photoperiod).

Maize stalk tissues were collected and frozen immediately in liquid nitrogen, and total RNA samples were isolated using the Plant RNA extraction kit (TIANDZ, Inc., Beijing, China). One microgram of total RNA was converted to cDNA using HiScript II RT SuperMix for qPCR (Vazyme, Nanjing, China) according to the manufacturer’s instructions. qRT-PCR reactions were performed with the Bio-Rad CFX96™ real-time PCR system using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). Maize Ubc9 was used as an internal control (forward primer, 5′-acgaaggtcttccatccaaacatc-3′; reverse primer, 5′-gtttcatgggcaacaccacaatcg-3′), and the relative expression levels for target genes were calculated using the 2^–ΔΔCt method. qRT-PCRs were all performed with three biological replicates. The primers used for ZmVPS37A are as follows: ZmVPS37A-F 5′-gagctggagatttcctctcttc-3′ and ZmVPS37A-R 5′-gaagactgcggactatctgaagc-3′.

Prior studies by Paul et al. (2014) and Li et al. (2015) defined the maize genomic loci that encode components of secretory and endocytic pathways and autophagy pathways, respectively. Here, with the improved maize genome sequence available [RefGen_v4, Jiao et al. (2012)], we exhaustively search for membrane-trafficking components using BLASTP and TBLASTN¹⁴ with the full-length amino acid sequences of Arabidopsis membrane-trafficking factors as baits (Bassham et al., 2008; Paul et al., 2014; Gao et al., 2017). By sequence alignment and phylogenetic and domain analyses, the membrane-trafficking genes in the maize genome are defined (Figure 1, Supplementary Figures 1–9, and Supplementary Datasets 1–7). Totally, we obtained 479 maize genes involved in membrane trafficking, including 11 genes encoding components for clathrin-coated vesicle (CCV), 35 for AP1-5 complexes, 26 for COP-I complex, 30 for COP-II complex, 10 for ARF GTPases, 52 for autophagy machinery, 63 for ESCRTs machinery, 96 for tethering complexes, 11 for retromer complex, 76 factors for Rab GTPase protein family, and 69 for SNARE protein family.

The factors are highly comparable in their protein length and the number of orthologs between maize and Arabidopsis (Supplementary Figures 1–9 and Supplementary Datasets 1–7). However, we observed a significantly higher copy number for certain factors in maize than in Arabidopsis. For example, there are seven, four, and three gene copies for the B-COP β’ and ε subunits, and F-COP γ subunit, respectively, while there are only three, two, and one genes for the relevant Arabidopsis factors. Another clear example is the genes encoding clathrin heavy chain proteins, with eight and two gene copies in maize and Arabidopsis genomes, respectively. Notably, an interesting exception is that Arabidopsis possesses a slightly expanded set of Rab11 (RabA), with 26 paralogs as compared to 24 paralogs in the maize genome (Supplementary Figure 8 and Supplementary Dataset 6). Moreover, we could not detect any orthologs of PROS (an ESCRT component) and RABA6 protein in maize (Supplementary Datasets 3, 6).

Domain architecture is critical for protein function and is therefore mostly retained by orthologs during evolution for their function performance. Thus, we compared the domain architectures of maize factors with their relevant Arabidopsis components as an additional criterion to assign orthology. Using NCBI Conserved Domain (CDD¹⁵) and Pfam website,¹⁶ we observed that 81.84% (392/479) of maize factors possess the same domain architecture as their Arabidopsis orthologs, which indicates that the membrane-trafficking system is largely conserved in monocotyledonous and dicotyledonous plants. In addition, we found that 6.47% (31/479) of maize orthologous genes are predicted to encode truncated proteins with defective domain structures. For examples, four genes [Zm00001d009915 (ZmCHC3a), Zm00001d026086 (ZmCHC3b), Zm00001d021067 (ZmCHC4a), and Zm00001d012897 (ZmCHC5b)] are predicted to encode truncated clathrin heavy chain protein (Supplementary Figure 1), and two genes [Zm00001d013458 (ZmATG3b) and Zm00001d003156 (ZmATG3d)] are predicted to encode truncated autophagy component ATG3 (Supplementary Figure 4). Further gene expression analysis with MaizeGDB¹⁷ revealed that no cDNAs have been detected for these six genes, indicating that they are likely pseudogenes (data not shown).

In some cases, we observed the presence of additional domains in some maize factors when compared with their Arabidopsis orthologs, e.g., for SNARE protein SYP3 or ESCRT component SNF7. Arabidopsis Qa-SNARE protein AtSYP31 (AT5G05760) and AtSYP32 (AT3G24350) contain a t-SNARE domain at their C-terminal regions, while an additional syntaxin-5 domain is detected at the N-terminal regions of maize orthologs ZmSYP31 (Zm00001d008758) and ZmSYP32 (Zm00001d039614) (Supplementary Figure 9). A similar scenario is also detected for one of the maize ESCRT SNF7 proteins (Zm00001d044480), which contains an additional leucine-rich repeat N-terminal domain besides the C-terminal SNF7 domain (Supplementary Figure 5). These extra domains might provide another layer of regulatory feature for these maize membrane-trafficking factors.

To investigate the spatiotemporal expression patterns of maize membrane-trafficking genes, we performed a gene co-expression (clustering) analysis using a previously described RNA-seq dataset derived from 64 different tissues/stages of the maize inbred line B73 (Stelpflug et al., 2016; Figure 2 and Supplementary Figures 10–16). We observed 3–12 clusters in all protein complexes or families analyzed, indicating that the expression profiles of these genes varied considerably. Notably, with few exceptions (e.g., ZmGOS11), most genes encoding SNARE proteins are grouped into nine clusters, which largely coincide with their subcellular localization (Figure 2). The first and last clusters mostly contained genes encoding TGN/PM-localized SNAREs, including ZmSYP1, ZmVAMP72, ZmNPSN1, and ZmTYN11, but showed distinct expression patterns, indicating that these two sets of SNARE genes might have tissue-specific functions. The genes in the second cluster (TGN-Vac) displayed elevated expression in root and some reproductive tissues (silk and tassel), whereas the genes in the third (TGN-vac or PM), fourth (Golgi), and sixth (ER-Golgi-Vac) clusters showed enhanced expression in endosperm but not in embryo. Interestingly, we also detected increased mRNA abundance in endosperm for multiple membrane-trafficking genes (175 loci) other than SNAREs (Supplementary Figures 10–16). It was similar to the earlier observation for maize autophagy machinery (Li et al., 2015), indicating that the membrane-trafficking system, including the autophagy pathway, might participate in the maturation of maize endosperm during seed development. The fifth cluster contained seven SNARE genes (ZmSNAP33, ZmVAMP728, ZmVAMP729, ZmVAMP712, ZmSYP125, ZmSYP123, and ZmSYP23), which all expressed preferentially in anthers above all other tissues, implying their roles during anther and/or pollen development.

Besides the analysis of spatiotemporal expression patterns, our RNA-seq analyses also detected differential expression at the isoform level for multiple membrane-trafficking genes. For example, VAMP72, which is a subgroup of VAMP7 specific to green plants and is largely involved in trafficking at the PM, had eleven isoforms in the maize genome and could be subdivided into five groups based on their expression patterns (Figure 2). The first group, encompassing ZmVAMP727 and ZmVAMP7210, displayed elevated expression in root and leaf tissues. However, the second group, containing ZmVAMP721, ZmVAMP723, and ZmVAMP725, had enhanced expression patterns in whole seed and endosperm tissues, suggesting a role in regulating seed development. The third group included ZmVAMP728 and ZmVAMP729; they were expressed exclusively in anther tissue, implying their specific role in pollen development. The fourth group contained ZmVAMP722 and ZmVAMP726, which were expressed preferentially at the very early stage of seed development. The last group (ZmVAMP724a and ZmVAMP724b) had ubiquitous expression patterns.

Our spatiotemporal expression analysis revealed that 257 membrane-trafficking genes displayed elevated expression in the endosperm, embryo, or whole seed tissues, suggesting that they might play important roles in seed development (Figure 2 and Supplementary Figures 10–16). To further investigate the contributions of these genes in maize seed development, we performed a weighted gene co-expression network analysis (WGCNA) to analyze the co-expression relationships among 16,888 differentially expressed genes identified from the above-mentioned RNA-seq dataset (Stelpflug et al., 2016). In total, 27 color-coded gene modules were identified according to the similarity of gene expression patterns (Supplementary Figure 17A and Supplementary Table 1). The module-tissue correlation analysis showed that modules turquoise and yellow positively correlated with early and late endosperm development, respectively (r = 0.36, p = 9.45e–4; r = 0.39, p = 3.44e–4, respectively) (Supplementary Figure 17B). The genes of module turquoise peaked at the endosperm of 12 DAP (days after pollination), whereas those in module yellow were expressed highly at the endosperm of 24 DAP (Supplementary Figures 17C,D). Moreover, module black was correlated with embryo development (r = 0.66, p = 2.46e–11), and module dark red was correlated with early seed development (r = 0.56, p = 5.56e–8).

In early endosperm-specific module turquoise, 114 genes encoding membrane-trafficking factors were identified among 3,850 annotated genes, including 18 for CCV and AP complexes, 15 for COPI complex, 15 for COPII complex, three ARF GTPases, six for ESCRT complex, three for retromer complex, 30 for tethering factors, 16 for Rab GTPases, and 8 for SNARE proteins (Supplementary Table 1). Further analysis of the subnetwork harboring components with a weighted value over 0.2 revealed that genes encoding five COPII components (ZmSec13a, ZmSec13b, ZmSec24a, ZmSec23i, and ZmSec23j), four Rab GTPases (ZmRABA1f, ZmRABE1b, ZmRABE1c, and ZmRABE1k), one Golgi/PM localized SNARE protein (ZmSec22), and one ESCRT-II component (ZmVPS25) have higher connectivity with other genes (Figures 3A,B). The COPII complex is responsible for the transport of newly synthesized protein from ER to cis-Golgi, and RabE GTPases are well-known regulators for vesicle transport from the Golgi apparatus to the PM (Speth et al., 2009). The GO enrichment analysis revealed that the genes involved in protein biosynthesis and assembly, ribosome biogenesis, and translation apparatus are overrepresented (Figure 3C), suggesting the role(s) of these membrane-trafficking components in the delivery of nascent biosynthesized proteins during early endosperm development.

To investigate the membrane-trafficking genes responsible for pathogen defense during Gibberella ear rot disease caused by F. graminearum, we also performed a WGCNA to identify gene modules between resistance inbred line CO441 and susceptible line B73 (Kebede et al., 2018). In total, 15 modules were identified among 12,528 differential expression genes across eight experimental conditions (Supplementary Figure 18A and Supplementary Table 2). According to the module-trait correlation analysis (Supplementary Figure 18B), modules red and magenta showed a highly positive correlation with disease defense of resistance line CO441 (r = 0.54, p = 0.04; r = 0.57, p = 0.03, respectively). The genes within module red were highly expressed in resistance line CO441 at 1–2 days after fungi infection (Supplementary Figure 18C). In contrast, genes within module green were highly expressed in susceptible line B73 (r = 0.7, p = 0.04) and peaked at 2 days after fungi infection (Supplementary Figure 18D).

In module green, four membrane-trafficking genes, which encode two ESCRT components (ZmTOL9 and ZmISTL6), one SNARE protein ZmSYP121, and a SNARE-interacting protein ZmSec1b were identified among 603 annotated genes with a weighted value over 0.3 (Supplementary Figures 19A,B). These genes were strongly co-upregulated with several defense genes in B73 2 days after fungi infection (Supplementary Figure 19C). Within the subnetwork containing these four genes, genes involved in defense, SA biosynthesis, jasmonic acid signaling, and amino acid biological process are overrepresented (Supplementary Figure 19D), implying the role of these four factors in response to pathogen infection. The identification of ZmSyp121 here further confirmed its conserved role in pathogen defense in both monocotyledonous and dicotyledonous plants (Collins et al., 2003).

We then analyzed module red and identified eight membrane-trafficking genes among 671 annotated genes (weighted value > 0.15) (Figures 4A,B). These factors all showed a stronger and more rapid gene expression in resistant line CO441 than in B73 upon pathogen infection, along with many defense genes, such as PR1 and PR10 (Figure 4C). Among these genes, those encoding SNARE protein SYP131, Rab GTPase RABH1c, and autophagy component ATG18f showed a higher degree of connectivity (>10). The Arabidopsis orthologs of SYP131 and RabH1b have been reported to be involved in the vesicle trafficking of callose (Gamir et al., 2018; He et al., 2018). The identification of ZmSYP131 and ZmRABH1c here suggested that maize orthologs might have a similar function in callose deposition upon fungi infection. The identification of ZmATG18f also suggested the role of autophagy in the pathogen defense of maize. Further, GO analysis of the subnetwork containing these eight genes revealed an enrichment of defense genes (Figure 4D). Moreover, GO terms related to the biological process of biotic stress such as lipid biosynthetic process, isoprenoid biosynthetic process, and trypotophon metabolic process were also specifically enriched, indicating that membrane-trafficking factors likely coordinated with the defense-related metabolic processes and transported defense proteins and antimicrobial metabolites to the infected regions. In summary, the WGCNA analysis above revealed that the maize genes encoding membrane-trafficking factors were upregulated together with defense-related genes, which conferred basal resistance to F. graminearum in B73, but also with genes involved in various secondary metabolic processes to further improve resistance capability in CO441 line.

The identification of selected genetic loci during domestication and improvement of crops is important for breeding elite varieties and can further facilitate appreciation of the contributions of selected loci to important agronomic traits (Zhou et al., 2015). Molecular evidence for selective sweeps on multiple membrane-trafficking genes during maize domestication and improvement has been reported in several previous studies (Hufford et al., 2012; Jiao et al., 2012; Liu et al., 2015). With more comprehensive and better annotated B73 genome datasets and the complete list of membrane-trafficking genes, we searched the larger and more diverse maize HapMap 3 genotype dataset (Bukowski et al., 2018), which covers 20 wild relatives (Zea mays spp. Parviglumis), 23 landraces, and 59 improved cultivars (Supplementary Table 3), for genomic regions showing positive selection signals during domestication and improvement through statistical analyses of high population differentiation (Fst) and low level of nucleotide diversity (π) (Table 1). The top 5% of Fst and π ratio were set as thresholds to identify potentially selected regions. During maize domestication (teosinte vs. landrace), 15,760 windows were found as candidate regions (Fst ≥ 0.26; π teosinte/π landrace ≥ 1.55). A total of 25 out of 479 membrane-trafficking genes were found to be associated with domestication (Figures 5A,C and Table 1), including five genes encoding components of the autophagy machinery, three for ESCRT complex, four for Coat and AP complexes, one for retromer complex, four for tethering factors, six for ARF and Rab GTPases, and two for SNARE proteins. In the comparison of maize landrace and improved populations, 14,729 slide windows were revealed by two methods (Fst ≥ 0.08; π landrace/π improved ≥ 1.39). Similarly, a total of 14 membrane-trafficking genes were identified to be improvement genes (Figures 5B,D and Table 1), including two for autophagy factors, three for SNARE proteins, two for Rab GTPases, two for tethering factors, three for coat complexes, and two for ESCRT components.

Out of the 39 selected membrane-trafficking genes identified in our study, 22 have been identified in previous evolutionary analyses (Hufford et al., 2012; Jiao et al., 2012; Liu et al., 2015). Moreover, nine out of 39 selected membrane-trafficking genes were detected by QTL mapping or genome-wide association studies (Supplementary Table 4), which further implicates the potential roles of these genes in maize development and biotic and/or abiotic stress responses. For example, four of them are implicated in biotic stress responses [Zm00001d042303 (ZmATG6a), Zm00001d027962 (ZmVPS37A), Zm00001d015102 (ZmGOS13), and Zm00001d011347 (ZmATG7)], and three for ear and seed development [Zm00001d011984 (ZmATG6b), Zm00001d042215 (ZmATG18e), and Zm00001d009150 (ZmSec23a)], and one for leaf senescence [Zm00001d042215 (ZmATG18e)].

Among 114 genes involved in seed development, ZmSec23a on chromosome 8 was also found to have strong selection signals with high-level population differentiation between landrace and improved lines (Fst = 0.211, π landrace/π improved = 1.721) (Figure 6A and Table 1). Given that ZmSec23a is an important COPII component that might regulate the protein biosynthesis during seed development, we first analyzed its mRNA levels using previously reported RNA-sequencing data (Stelpflug et al., 2016). As shown in Figure 6B, ZmSec23a showed elevated expression in the endosperm and seed tissues. Interestingly, we also found that ZmSec23a exhibited increased mRNA abundance in high oil lines compared with regular materials (Figure 6C).

To gain deeper insights into the evolution of ZmSec23a, we then performed haplotype analysis with 468 selected teosintes (20), landraces (22), improved regular lines (403), and high oil lines (23) and identified 24 different haplotypes in its coding region. Among these haplotypes, Hap_20 and Hap_22 are enriched with high oil lines (Figure 6D), with 18 out of 22 being high oil lines. These two haplotypes not only carry the lead SNP (chr8: 38521846; G/T allele) detected by a previous genome association study (Li et al., 2013) but also carry the G/A, C/T, and T/C variants. Interestingly, these variants were also found in Hap_4, which mainly consists of teosintes and landraces, but not modern cultivars (Figure 6D). However, the haplotype network analysis failed to detect the close relationship between Hap_4 and Hap_20/22 (Figure 6E). Moreover, several alleles (Chr8: 38526804, T/C; Chr8: 38526831, G/A; Chr8:38528117, T/C; and Chr8: 38530903, C/T) within the gene region of ZmSec23a have different degrees of linkage disequilibrium with the lead SNP (G/T allele) from domestication to improvement, or from regular inbred lines to high oil materials (Supplementary Figure 20), indicating that these alleles of Zmsec23a could also be associated with high oil content.

Another noticeable selected gene is ZmVPS37A. As a component of ESCRT machinery, Arabidopsis VPS37-1 is required for the endocytosis of plant immune receptor FLS2 after flagellin elicitation and for flg22-triggered stomatal closure to prevent bacterial entrance (Spallek et al., 2013). Moreover, ZmVPS37A is significantly associated with variation in Rp1-D21-induced HR (Olukolu et al., 2014). In this study, we found that the genetic diversity of ZmVPS37A in landrace was significantly lower than that in teosinte (Fst = 0.589, π landrace/π teosinte = 2.154, Figures 5, 7A and Table 1), indicating that it is a domestication gene. To understand why ZmVPS37A is selected, we infected the stem tissues of teosinte (W71-2), landrace (303WX), and B73 seedlings with F. graminearum. As shown in Figure 7B, teosinte exhibited better resistance than 303WX and B73, while B73 displayed obvious lesions and severe stalk rot symptoms. qRT-PCR analysis revealed that the transcript level of ZmVPS37A in teosinte and 303WX was about 1.54-fold compared to B73 control after F. graminearum infection (Figure 7C). Analysis of the cis-regulatory elements in the promoter regions of ZmVPS37A further identified more stress-related elements in teosinte than in B73 (Supplementary Figure 21). These observations suggest that ZmVPS37A, like its Arabidopsis homolog, is critical for pathogen defense and likely lost its defense function during domestication as a trade-off between plant growth and defense.