Genetic basis of Fusarium ear rot resistance and productivity traits in a heterozygous multi-parent recombinant inbred intercross (RIX) maize population

Abstract

Maize (Zea mays L.) is one of the most productive crops worldwide. As a heterotic crop predominantly grown as F₁ hybrid, maize exhibits challenges for genetic studies of complex traits, since homozygous genotypes, which are largely used in these studies, may not accurately reflect what happens in cultivated conditions. To map Fusarium Ear Rot (FER) resistance to Fusarium verticillioides and traits with potential impact on yield, including phenology, we constructed a recombinant intercross (RIX) population. This was achived by crossing pairs of recombinant inbred lines (RILs) derived from a multi-parent maize population. We characterized the RIX population over two growing seasons, employing artificial F. verticillioides inoculation. The heterozygous background of the material enabled the identification of QTL and candidate genes through in silico reconstruction of RIX genotype probabilities. A total of 37 loci were identified using single-year BLUPs while 29 with joint-year BLUPs. These, included several known QTL associated with days to tasseling, kernel row number and a QTL on the chromosome 9 associated with FER resistance. In this region, we could identify candidates based on their predicted functions and potential roles in plant-pathogen interactions and/or resistance mechanisms. These QTL represent potential breeding targets to FER resistance and yield components in commercial maize varieties.

Background

Maize (Zea mays L.) leads the global cereal production with 1.16 billion tons harvested in 2022 [1]. This versatile crop serves multiple purposes: it is a staple food for humans, a primary feed for livestock, and an important raw material for manufacturing and biofuel production. With the world population projected to rise, maize production must increase to meet the growing demands for food, feed, and biofuel, which are essential for achieving net zero emission targets [2, 3]. To keep pace, modern maize breeding techniques must raise their efficiency. Nowadays, selection strategies predominantly rely on genomic information and molecular markers, making the breeding process more accurate and potentially faster [4].

Typically, in maize, mapping of genes responsible for traits that enhance productivity and/or quality is conducted using fully homozygous genotypes, as these make quantitative trait locus (QTL) mapping more tractable [5]. However, these materials may not fully represent the genetic makeup found in field-grown maize varieties, which are typically highly heterozygous plants. The use of these genotypes in cultivated maize reflects the heterotic mechanisms that contribute to the superior performances of F₁ hybrids over their parental genotypes [6]. Indeed, maize populations composed of hybrids should be used to map QTL as they have a genetic background more similar to that of commercial varieties. To date, a few forward genetic studies employed diallel cross designs to develop panels of heterozygous maize materials. These panels have been used to map classical phenotypes such as plant height [7], phenology [6], as well as performance-related traits [8], demonstrating that this approach is doable and efficient.

The ability to pinpoint the genetic basis of complex traits largely depends on the diversity available in the mapping material and by the density of recombination events, which contributes to enhancing the mapping precision. Multi-parent Mapping Populations (MPPs) of recombinant inbred lines (RILs) provide an ideal platform to increase diversity and maximize recombination, using resolved pedigrees that enable haplotype-based mapping [9]. Past studies conducted in animal models used MPP to develop Recombinant Inbred Intercross (RIX) − F₁ hybrids that are created by controlled intermating of pairs of MPP RILs − with the objective of further increasing recombination density in a heterozygous background [10]. Applying this approach to crops like maize can help create panels of hybrids with traceable pedigrees, enhancing the ability to map complex phenotypes that impact maize productivity, and quality.

In maize, the accurate mapping of QTL of agronomic importance is made urgent by a push for sustainable intensification of production. QTL markers and relative candidate genes can be used to tackle complex traits that have seen limited potential for improvement, such as disease resistance. Fusarium verticillioides, the causal agent of Fusarium Ear Rot (FER), is a widespread mycotoxin-producing fungus that limits maize cultivation, and maize resistance represents a trait that should be prioritized. The infection process is complex: its spores can survive extreme conditions and infect maize throughout the cultivation cycle, making agronomic practices largely ineffective for disease control [11]. Effective control requires frequent treatments due to the limited resistance of commercial varieties [12]. Indeed, enhancing host resistance to FER is an optimal control strategy; studies have demonstrated that suitable genetic materials can help identifying genomic regions and molecular markers that can contribute to FER resistance [13]. Perspectives of maize improvement for disease resistance should always be coupled with agronomic productivity, and traits like yield components and phenology should also be studied contextually. There is value in carefully addressing these traits in open field experiments to capture their association and gene by environment interactions.

In this study we produced a new RIX population made of 214 F₁ genotypes derived from the MAGIC maize MPP [14], with the scope of mapping traits of agronomic relevance in a heterozygous background. We derived RIX genotypes in silico starting from the genome reconstruction of the MAGIC RILs. Indeed, we could use the RIX as a platform to measure and map FER resistance following artificial inoculation of F. verticillioides, phenology and ear traits that can contribute productivity. We found that the RIX panel is an ideal, unstructured population of hybrids characterized by wide phenotypic variability. We mapped 29 significant associations, underlying QTL for the studied traits including a major genomic region for FER resistance on chromosome 9. We identify and discuss potential candidate genes, including those associated in well-characterized QTL. The development of RIX demonstrates that creating panels of hybrids from MPPs can yield valuable insights into the genetic control of complex traits.

Results

Genetic Diversity in the RIX is uniformly distributed

The production of the RIX panel followed a chain-crossing scheme where MAGIC RILs were intercrossed in random pairs (Fig. 1A) to establish a panel of 214 RIX hybrids, combining the genomes of 229 MAGIC RILs (Table S1).

Fig. 1

A) Schematic representation of the chain-crossing scheme used to produce the RIX panel: Each RIL was crossed with its adjacent RIL. Exceptions to this scheme occurred due to differences in flowering times. B) Kinship represented as heatmap of pairwise similarities between the 214 RIX genomes reconstructed starting from MAGIC maize RILs: white and red represent low and high similarities, respectively. On the left and top side, a hierarchical tree represents relationship between RIX. Samples were ordered according to R base function hclust(). The histogram in inset represents the count of occurrences (i.e. similarity, y axis) vs the similarity value subdivided by 1,000 bins

The MAGIC RILs were genotyped with Single Primer enrichment Technology (SPET) [15] and the RIL genotyping data of each crossing pair were used to derive the corresponding RIX genome in silico, relying on reconstructed genotype probabilities. The RIX cryptic relatedness was used to measure the diversity in the panel, revealing low structure (Fig. 1B). The distribution of the pairwise dissimilarity values (Fig. 1B, inset) confirms that the panel of hybrids had no significant genetic structure, although groups of individuals, that share a common parent, highlight small clusters, mainly composed of three to six hybrids, when multiple crosses involved the same parents (Fig. 1, darker red squares along the diagonal).

There is large diversity in the RIX for FER and components of production

The phenotyping of the RIX panel was conducted over two consecutive years (2018 and 2019) in Northern Italy (San Bonico, Piacenza, 45.01 N, 9.70 E). We measured Days to Tasseling (DT), number of kernel rows per ear (Krows, n), ear length (EL, cm), ear diameter (ED, mm), cob diameter (CD, cm), and volumetric weight (Vweight, Kg/hL) and FER score based on visual scoring of artificially inoculated cobs (FER_score, 1 to 7, where 1 indicates no infection). The phenotypic values, including FER_score, measured in each experiment were used to extract Best Linear Unbiased Predictors (BLUPs) using single-year and joint-year models (Table S2). Joining molecular based measures of relatedness via kinship and BLUPs, we were also able to derive estimates of narrow-sense heritability (Table S3).

The top heritability values were obtained for CD and DT, with values of 0.78 and 0.77, respectively. The heritability of FER_score across years was 0.52. Krows, CD and EL and ED displayed strong and positive correlations (Fig. 2A). Correlation of FER scores was positive over the two years, but no significant correlation was observed with the other traits, except a significant negative correlation with Vweight (Supplemental Fig. S1).

Fig. 2

A) Pairwise complete observation correlations BLUPs of the MAGIC maize RIX population measured in 2018 and 2019. Non-significant coefficients (P < 0.05) are marked with an ‘x’. B) PCA of scaled joint-model BLUPs of 214 MAGIC maize RIX. Individual RIX are depicted as gray dots while the MAGIC founder lines are higligheted with colors. Abbreviations: DT: Days to Tasseling, Krows: n° of kernel rows; EL: Ear length; ED: Ear diameter; CD: Cob diameter; Vweight: Volumetric weight; FER_score: Fusarium Ear Rot severity scores; C) boxplots of joint model BLUPs of Ear Diameter and D) FER score, with color coding according to panel B

Using joint-model BLUPs, positive correlations were observed among ear traits and Krows (Supplemental Fig. S1). In a principal component analysis (PCA) of phenotypic diversity, the first component separated the RIX by ear dimensions while the second component ordered the RIX by FER score, phenology and kernel volumetric weight (i.e. Vweight). Altogether these components explained about 53% of the total phenotypic variation (Fig. 2B). The PCA and the distribution of the measured traits (Fig. 2C and D) clearly display that the RIX were phenotypically transgressive concerning the measured traits when compared to the founders of the MAGIC population.

The observed differences were especially pronounced in traits like ED (Fig. 2C, Supplemental Fig. S2), which significantly distinguished the founder lines from the RIX in the PCA (Fig. 2B), but not in FER scores (Fig. 2D). A putative heterotic effect of RIX compared to the homozygous MAGIC MPP founder lines was recorded for traits of ear and cob, with most hybrids consistently exhibiting transgressive phenotypes relative to the founders (Fig. S2 C, D and E). The same trend could be observed for Krows, where the majority of founders laid in the first quartile of the distribution of the joint-model BLUPs (Fig. S2 B). Regarding phenology, the earliest genotypes were exclusively hybrids, while the earliest founder, F7, was located in the second quartile of the distribution (Fig. S2 A). For Vweight, the founders displayed values consistently close to the median, with the distribution’s tails occupied exclusively by the RIX (Fig. S2 F). Regarding the distribution of FER scores, the most tolerant genotypes were consistently RIX (Fig. S2 G). No significant correlation was found between DT and FER scores (Fig. S1).

Identification of QTL for FER and components of production

The reconstructed RIX genotype probabilities, along with BLUPs derived from the phenotypic characterization of the studied traits, were used to map QTL. We identified 37 significant signals of associations using single-year BLUPs, some of which overlapping across years (Supplementary Table 1). Using joint-model BLUPs, we identified 29 significant peaks. We focused our attention on a subset of 14 loci by selecting those spanning confidence intervals shorter than 10 Mbp. The definition of confidence intervals is a function of the maximum level of significance (LOD score peaks) and of the founder effects estimated by the model. The selected regions are those where we could have the highest statistical support and are that are likely to contribute for a relatively higher variation of the traits (Table 1). The confidence intervals of the prioritized QTL had an average size of 3.08 Mb, ranging from a minimum of 0.3 Mbp up to 8.55 Mbp (Table 1). The smallest and largest intervals were both associated with the Krows, located on chromosomes 4 and 7, respectively. The selected QTLs were further analyzed to identify loci or suggestive candidate genes based on evidence from literature, predicted annotation and/or functional characterization of Arabidopsis and rice (Oryza sativa L.) orthologs within the mapping intervals (Table 1).

Table 1 Summary of the QTL mapped using the RIX population spanning over regions < 10 Mbp

The top associations, in terms of LOD score, were recorded for DT on chromosome 8 (Fig. S3 A). For these signals, we defined 4 main confidence intervals supported by a LOD drop equal to 3 (Table 1), occupying an overall chromosomic portion spanning about 10 Mbp (Fig. S3B). In this genomic region, the F7 haplotype consistently contributes to earliness (Fig. S3 C), similarly to what observed in MAGIC RILs [14]. The QTLs detected on the RIX correspond to two well-characterized loci: the phosphatidylethanolamine-binding protein (pebp8, Zm00001 d01075), which is associated to the vegetative to generative transition 2 (vgt2) locus [16] laying at about 126 Mbps on Chr 8 (Fig. S2); and the vgt1 locus, where the gene model ZmRap2.7 (Zm00001 d010987) has been previously identified and described for its role in flowering time regulation [17] (Table 1).

Regarding the traits related to the ear, specifically Krows, CD, and ED, the strongest association was found on chromosome 7, with a peak at approximately 137.3 Mb (Table 1). The confidence interval of this QTL spans over 1.5 Mbps, and encompasses 39 gene models. The QTL on chromosome 4 at 208.17 Mb for Krows contains krn4 (Zm00001 d052889—kernel row number4), a gene know to contribute to the variation of this trait [18] (Table S4). Krows was mapped to another QTL on chromosome 1 at 286.3 Mb, approximately 13.6 Mb upstream of the gene Zm00001 d034629, also known as tasselseed6, which has been previously associated with the kernel row number 1 (krn1) locus [19, 20]. Vweight was mapped in five regions across four chromosomes, with the highest LOD score found on chromosome 9. This locus contains 17 gene models and does not overlap with the previously known kw9 locus [21].

For FER we identified a QTL on chromosome 9 (Fig. 3A and B), which contains 70 gene models (Supplementary Table S5). Among these, there are several candidate genes having functions potentially associated with resistance mechanisms (Table S5).

Fig. 3

QTL mapping results of the Fusarium Ear Rot (FER) score phenotyping; A) Genome-wide LOD score plot. B) LOD score plot of chromosome 9; in both figures the red and the green dotted lines represents the threshold of significance of the 99^th and 95^th percentiles, respectively, of the permuted LOD distribution, based on 999 permutations. C) Effect, expressed as coefficient, of each founder at the QTL confidence interval

In this QTL, the pattern of founder effects splits into three groups (Fig. 3C). The haplotype effect of H99 appears to drive the QTL and contributes to FER susceptibility in this region. Observing the phenotypic distribution, the H99 founder, however, does not show any transgressive expression of FER tolerance, being in the second quartile of the distribution of joint model BLUPs (Fig. 2D). Conversely, another group of founders’ haplotypes has positive effects and likely to increase tolerance, with the haplotype of A632 producing the best effect in the population. Despite this, the A632 founder genotype itself did not display a notable phenotypic value underlying possible resistance, remaining close to the population median in terms of FER values (Fig. 2D).

Discussion

MPPs are invaluable for generating genetic variability, enabling researchers to unravel the molecular basis of complex phenotypes. Creating RIX starting from RIL within MPPs can enhance the mapping definition while reducing the number of genotypes to be screened when compared the corresponding RIL panel, due to the combination of different pairs of genomes with traceable pedigrees, into single individuals. Yet in a RIX panel, a complex network of relationships may arise, especially when different hybrids share a common parent, as it happens by design in a population like that used in this study (Fig. 1A). These underlying relationships may cause genetic structure which is not typically observed in RIL populations [10]. The results of our study showed that the hybrids within the MAGIC maize RIX panel − built using a chain-crossing scheme (Fig. 1A) − have relatively limited relatedness and marginal genetic structure (Fig. 1B). From a phenotypic standpoint, the RIX exhibited broad variation: hybrids were transgressive in the expression of several phenotypes, particularly those related to fitness (i.e. ear dimension), relative to the inbred founder lines.

From a breeding perspective, studying these phenotypes using hybrids is relevant, as it may lead to the identification of loci that contribute to the expression of specific traits, such as FER resistance or those contributing to yield in heterozygous genomic backgrounds. RIX panels can also constitute an invaluable resources for investigating the mechanisms underlying heterosis in maize, which remain largely elusive. To enhance our understanding of heterosis, it is crucial to include the RIL parents used to produce the hybrid population in the same field experiments. This would allow to measure the hybrid’s performance versus the average performance of both parents (mid-parent heterosis) or the better-performing parent (best-parent heterosis). The crossing scheme employed in our study (Fig. 1A) can make these types of experiments more cost-effective by minimizing the number of RILs and reducing potential issues related to underlying genetic structure.

To evaluate the effectiveness of our RIX panel for QTL mapping, we started with days to tasseling (DT), a key trait in maize breeding used to improve adaptation to various environmental conditions. We were able to successfully type known regions, including those on chromosome 8 harboring Vgt1 and Vgt2. It was also possible to appreciate the potential contribution and/or interaction of other loci that encompass genes, which play a role in flowering time [22] and for which we observe a potential role in hybrid genomic backgrounds. Indeed, we could map peaks on the chromosomes 1, 3, 5, 6 and 7 (Fig. S3) (Table S5). Among the genes identified within the confidence intervals of these QTL it is noteworthy to mention that we mapped Zm00001 d021748 on chromosome 7, which is ortholog to Arabidopsis’ FLOWERING BHLH 4, which is a known transcriptional activator of the photoperiodic flowering regulator CONSTANS [23].

Concerning the mapping of traits potentially contributing to yield, we could map kernel row number, ear size and kernel specific weight. Interesting candidate genes were identified for these traits. An example is represented by te1—terminal ear1 locus (Zm00001 d042445) [24], which has been extensively studied maize plant development, leaf initiation and its potential contribution to improve yield, which in our case was found associated with Krows, in the chromosome 3 QTL, with its peak located at 169 Mbp.

FER and fumonisin accumulation poses significant challenges to maize cultivation, as still today, current commercial hybrids remain susceptible or, in certain conditions, highly susceptible to ear rots. In recent years, many studies have focused on breeding superior genotypes and searching for QTL and candidate genes that may contribute to improving resistance against these pathologies. Pathogen invasion is among one of the main stresses that plant cells face in the first stage of an infection. This event is known to trigger the expression of stress-responsive genes.

Our QTL on Chromosome 9 encompasses putative candidate genes that were identified based on their functions and their potential role in explaining the observed variations. As these are candidate genes, further characterizations and experiments are necessary to understand and validate their roles. The first putative candidate gene is a Caleosin/Peroxigenase gene (Zm00001 d044732), the corresponding protein of which is known to catalyze the production of several compounds, including oxylipins, which are signal molecules that regulate plant defense [25]. Several authors attempted to better elucidate the oxylipins-mediated plant defense response against F. verticillioides by employing inbred lines and lipoxygenase knock-out mutants as well [26,27,28,29]. In this regard, Battilani and colleagues [30] observed that the disruption of a maize 9-lipoxygenase resulted in increased resistance to F. verticillioides and reduced fumonisin levels. Another candidate gene that may contribute to the observed variation at the chromosome 9 QTL is the ortholog of Arabidopsis NFX-L1 (Zm00001 d044726), which may act as activator of response mechanisms of plants to phytotoxins. In Arabidopsis, the ortholog of Zm00001 d044726 is a negative regulator of the defense response, which is induced by the presence of trichothecenes and T-2 toxins whose artificial application resulted in the occurrence of hypersensitive response in mutant plants [31, 32]. In Triticum aestivum, the ortholog TaNFXL1 is involved in the mechanisms that repress F. graminearum resistance [33]. Although F. verticillioides cannot produce trichothecenes [34, 35], it would be interesting to test whether other toxins, using, including those produced by F. verticillioides, contribute the resistance repression. In plants, barriers to pathogen invasion can rely on the accumulation of secondary metabolites like anthocyanins whose synthesis is promoted by GLKs (Zm00001 d044785), which has been reported to be highly expressed in resistant maize genotypes [36] or via the regulation of abscisic acid, for which we also identify a putative candidate (Zm00001 d044802). Other interesting candidates harbored on the Chromosome 9 QTL are Zm00001 d044775 and Zm00001 d044752 that encode for an enzyme involved in xyloglucan modifications. The xyloglucan is the main hemicellulose component and, in maize, the release of ferulates from hemicellulose has been associated to resistance against F. graminearum and F. verticillioides [37] also with implication on fall armyworm resistance. Physical barriers depend also on callose deposition by the action of DRP2B and DRP1 A proteins (Zm00001 d044757) [38] Previous literature on FER resistance reported several QTL on chromosome 9 associated with resistance measured on kernel; among these meta QTL on chromosome 9, the ZmMQTL9.1 co-localize with our QTL, supporting our finding. FER is a disease significantly influenced by environmental factors, which can mask genetic effects [12], potentially leading to varying disease incidence across consecutive years. The RIX panel revealed a moderate but significant correlation among the two experimental seasons in agreement with previous findings [39].

Interestingly, when screening MAGIC maize RILs for FER resistance, but at the seedling stage, Septiani and colleagues [40] mapped different resistance loci that were not recovered in the RIX. This result suggests that the same fungal strain may have different effects at various developmental stages, with other loci potentially playing a more significant role, as observed by Stagnati and colleagues [41], or as result from the interaction between additive and non-additive (i.e. epistatic) effects. Even though we have no FER resistance data produced on MAGIC RILs, we could observe that hybrids outperform parental lines in terms of FER resistance. This trend would likely be expected when testing RILs as well.

To understand potential sources of positive alleles we estimated the founder haplotype effect at the QTL region on chromosome 9. We identified the H99 parental haplotype as a negative contributor to the expression of FER resistance, likely driving the identification of the QTL. This haplotype originates from the Illinois Synthetic 60 C population. In Italy, where our experiment was conducted, this population has been used to develop lines, like Lo924 and Lo977, aimed at improving yield and tolerance to Maize Dwarf Mosaic Virus [42]. These materials were released in 1983 by the Maize Station in Bergamo (now CREA-CI). The haplotype that showed the best contribution to FER resistance was derived from the founder line A632. Developed in the ‘50 s and first released in 1963, A632 was originally selected for its resistance to stalk rot and corn borer [43]; the accession from which the line originates is also part of the Goodman diversity panel [44]. Stalk rot is caused by fungi of the genus Fusarium, while corn borers act as vectors, promoting the infection and spread of Fusarium species within stalks and ears. In our phenotyping approach, we mimic a mechanical infection with the pin bar-method to simulate a realistic insect attack. Interestingly, the A632 allele in this region improves resistance to FER triggered by this infection mechanism. We may speculate that this trait, for which A632 has been selected in the past, may have contributed to the observed FER resistance within the context of our RIX population.

Indeed, within the same QTL on chromosome 9, we also identified three genes (Zm00001 d044762, Zm00001 d044765 and Zm00001 d044763) encoding Anthranilate O-methyltransferase (aomt) isoforms (1, 2 and 3, respectively). Kumari and colleagues [45] − investigating maize response to mechanical wounding − reported the induction of the aomt1 gene in response to tissue damage. AOMT1 converts anthranilic acid into methyl anthranilate, which acts as a volatile compound induced during herbivory. Insect injuries play a significant role in spreading the disease, and using maize genotypes resistant to ear damage could help reduce disease incidence and mycotoxin contamination [46].

Mycotoxigenic fungi are a major threat for maize cultivation as they not only damage the plant and deteriorate the grains, but also produce toxins that are highly dangerous to humans and animals alike [47, 48]. Harvests having mycotoxin contamination above policy threshold are banned from food and feed use and FAO estimates that 25% of global maize production is contaminated by mycotoxins [49], with high economic impacts [50]. Mycotoxins-producing fungi also include Fusarium verticillioides, the causal agent of Fusarium Ear Rot (FER), a disease widespread in maize-growing regions worldwide. F. verticillioides spores can survive extreme environmental conditions and are able to infect maize plants throughout their entire cultivation cycle, rendering agronomic practices ineffective [11]. Indeed, effective disease control requires frequent treatments as commercial hybrids have limited resistance [12]. Improving host resistance would be an optimal means of controlling FER, and appropriate genetic materials may allow to identify promising genomic regions and molecular markers significantly associated with resistance variation [13, 51, 52].

Conclusions

Understanding the genetic basis of complex traits is crucial for advancing agriculture and ensuring food security. Traditional bi-parental mapping approaches are often challenged by limited genetic variability, few recombination events, and by the need of screening large panels to obtain high statistical power and resolution. Advanced genetic materials like Nested Association Mapping (NAM) and MAGIC populations, which incorporate greater genetic diversity, can mitigate these issues. Our study shows that RIX panels derived from MPPs can further improve these critical aspects, provide consistent results and allow for the measurement of loci contributions to the phenotypic variance directly into heterozygous genetic backgrounds. We believe that similar approaches to that used in this study, enabling the identification of valuable genomic regions that will support the ongoing maize breeding efforts.

Materials and methods

Development of RIX

The MAGIC population was initially developed from eight inbred lines: A632, B73, B96, F7, HP99, HP301, Mo17, and W153R. To address the failure of some crosses, a ninth founder line, CML91, was introduced [14]. The resulting population comprised over 1600 recombinant inbred lines (RIL). In this research, we derived RIX from an initial set of 229 randomly chosen MAGIC RIL. The production of RIX took place in Caleppio di Settala (45.44 N, 9.39 E, Milano, Italy) in 2015, at the experimental stations of CREI-CERZOO experimental station in San Bonico (45.01 N, 9.70 E, Piacenza, Italy), and in San Piero a Grado, (43.67 N, 10.35 E, Pisa, Italy) in 2017. In each location, RILs were sown in plots consisting of 25 plants in single rows and field management was conducted following standard agricultural practices. During pollination, a chain-crossing scheme was used, where each RIL was crossed with its adjacent RIL in the field. This method ensured that each RIL both received pollen from the preceding RIL and donated pollen to the following RIL. Deviations from the scheme occurred in instances when differences in flowering time between adjacent RILs were too large to enable a successful cross. In such cases, pollen from the respective RILs was used to pollinate other randomly selected plots. This implies that some RILs were used as pollen donor/receiver multiple times in the construction of the RIX panel. At full maturity, the products of the crosses were manually harvested and stored. Only RIX with enough seeds to enable the planned phenotyping experimental design were retained.

Phenotyping and artificial inoculation with F. verticillioides

Phenotyping on the successfully developed RIX was conducted over the growing seasons 2018 and 2019 at the CREI-CERZOO experimental station. The founders of the MAGIC maize population, along with the 214 RIX, were sown using a Randomized Complete Block Design (RCBD), with two replications. Each plot consisted of a single row of 25 plants, spanning 5 m. Rows were spaced 0.8 m apart with 1 m aisle on the hedges. A standard commercial maize hybrid variety surrounded the experimental field. The trials were sown on April 28 th in 2018 and April 17 th in 2019. Agronomic management was carried out with standard agricultural practices consisting of chemical weed management with ADENGO® XTRA (0.4 L/ha), nitrogen fertilization, 100 kg/ha of urea at the V4-V5 stage and irrigation as required.

A spore suspension of F. verticillioides ITEM10027 (MPVP 294) was produced to a final concentration of 1 × 10⁶ spores/ml and used for inoculation as previously described [27, 53]. Five plants were randomly chosen within each plot and subjected to artificial inoculation of F. verticillioides using the pin-bar method, i.e. perforating the ear with a three-pronged skewer carrying the spores [39, 54]. The inoculation took place 15 days after pollination, targeting the primary ear. The fungal strain is derived and maintained in the culture collection of the Department of Sustainable Crop Production, Catholic University of the Sacred Heart (Piacenza, Italy).

In both seasons, RIX were measured for phenology and production traits: DT (days), recorded when 50% of the plant in the plot showed visible silks, counting from the sowing date; after hand harvesting of the cobs, the following traits were recorded: Krows (number), EL (cm), ED (mm), CD (cm), and Vweight (Kg/hL) with a GAC 2100 cereal analyzer (DICKEY-john). Krows, EL, ED, CD were measured according to the Community Plant Variety Office (CPVO)TP2/3 protocol. Plants carrying the inoculated ears were hand-harvested separately from each plot. The inoculated ears underwent a drying process in a greenhouse for 2 weeks. Once dried and cleaned of the husks, the inoculated ears were evaluated for susceptibility to FER using a visual scoring approach. The severity of disease in inoculated ears was assessed by estimating the percentage of each ear covered by symptoms and pinkish or white mycelia. A seven-point evaluation scale was then used to categorize the severity where 1 indicated no infection, and the subsequent points corresponded to the following levels of infection: 2 (1–3%), 3 (4–10%), 4 (11–25%), 5 (26–50%), 6 (51–75%), and 7 (76–100%) [39, 55].

Phenotypic data analysis

Unless otherwise specified, all statistical analyses of this study were conducted using the R programming language. Best linear unbiased predictions (BLUPs) of the traits characterized in this study, including FER disease scores, were estimated using R/ASReml [56] with the model reported in Eq. 1:

$${y}_{ikm}= \mu + {g}_{i}+ {r}_{k}+{y}_{m}+g{y}_{im}+ e$$

(1)

Here, the observed phenotypic value is ${y}_{ikm}$, $\mu$ is the overall mean of the population, ${g}_{i}$ is the random effect for the i^th genotype, ${r}_{k}$ is the random effect for the k^th replication, ${y}_{m}$ is the fixed effect for the m^th year and $e$ is the error. For the calculation of BLUPs within a single year, data were subset and analyzed with a reduced model in accordance with Eq. 1. Correlation analysis and t-tests were conducted on pairs of BLUP distributions to determine trait similarity across different years. Pearson's pairwise correlations were calculated among BLUPs obtained from the two experiments, comparing each pair of phenotypes using all complete pairs of observations. A Principal Components Analysis (PCA) was conducted on BLUPs resulting from the joint model. To avoid unequal contribution of different phenotypes to the PCA, the joint-model BLUPs were scaled.

In silico genotyping of the RIX

DNA extraction and genotyping data were generated on MAGIC RILs as described in Ferguson, Caproni and collegues [15]. Briefly, RILs and founders of the MAGIC were germinated in rolled towels and genomic DNA was extracted from bulked leaf samples of five-day-old seedlings. The samples were genotyped at IGA Technology Services (Udine, Italy) with Single Primer Enrichment Technology (SPET), with a set of custom designed probes [57]. The resulting trimmed reads were aligned against the Z. mays reference genome version 4 [58]. Variant calling was done with the HaplotypeCaller algorithm [59], and SNPs were filtered based on a Phred-scaled score > 30 and SNP missingness per individual and locus below 20%. SNP loci that were homozygous on the MAGIC founders were used to estimate the genotype probabilities of the RILs that were used to produce RIX. A Hidden Markov Model (HMM) implemented in R/qtl2 [60] was used to produce a probabilistic reconstruction of each RIL genome as in [15]. RIX genotypes were then derived in silico: if the pair of RILs originating a RIX respectively had “B73” and “A632” founder haplotypes in a specific locus, being the RILs fully homozygous then the resulting RIX must be heterozygous for the “B73/A632” founder haplotype in that locus. This was done as follows: haplotype probabilities computed in each locus in each of the two RILs used to produce the RIX were halved and then summed together to derive a RIX-specific probability. After obtaining all the RIX genotype probabilities, we calculated the overall RIX kinship using the function calc_kinship() implemented R/qtl2. We represented this diversity as a heat map and dendrogram where the samples were ordered according to R base function hclust(). The measure of the relatedness among RIX also allowed to estimate the additive component of the genetic variance. We used this estimate to calculate narrow-sense heritability via a linear mixed model, using the hest_herit() function on R/qtl2 for all the phenotypes.

QTL mapping and candidate gene identification

QTL were mapped using a Linear Mixed Model (LMM) that incorporates kinship with the function scan1(), including random polygenic effects as in the Eq. 2:

$$y=\mu + {\beta }_{k}{X}_{k}+\sum_{i\ne k}{\beta }_{i}{X}_{i}+e$$

(2)

where y is the phenotype of all individuals, µ of the mean phenotype, β_k is the effect of the genotype X_k, and Σ_i≠kβ_iX_i combines the effects of all other genotypes except genotype k. The error (or residual) is represented by e. We estimated the effects of all other genotypes using the kinship matrix.

We defined significant peaks and QTL intervals, using the function find_peaks() specifying a LOD drop equal to 3 and significance level based on 999 trait-specific permutations, where the 99^th percentiles of the LOD permuted scores for each phenotype were set as phenotype-specific thresholds. Founder’s haplotype effects at each QTL were estimated by fitting a generalized linear model using the function scan1blup(), which calculates BLUPs of QTL effects, using a single-QTL model, treating the QTL effects as random. After identifying the confidence intervals, we prioritized those covering relatively short physical distances by setting a suggestive threshold of 10 Mbp. This approach ensured that candidate gene discovery was focused on genomic regions that did not exceed this threshold. In these regions, we derived all the annotated gene models from the reference genome V4 and used to match with the most similar Arabidopsis and rice orthologs, based on putative protein sequence identity Phytozome V13 [61]. This procedure, along with a thorough literature review, helped us identify suggestive candidate genes based on their putative functions and/or the functional characterization of their orthologs.