In the autumn of 2020, seeds of three soft red winter wheat varieties were planted into maize debris at two field sites in Illinois (Ewing: 38.096440, −88.845311; Urbana: 40.080776, −88.224160; Figure 1a) which were located over 225 km apart. At each site, a randomized complete block design was used where each of the four blocks contained one replicate plot each of all three varieties (i.e., 2 sites × 4 blocks ×x 3 varieties). The varieties planted were all AgriMAXX Wheat Co. varieties (492, 495, and 505; Mascoutah, IL), which have FHB resistance ratings of 5, 8, and 9, respectively (AgriMAXX internal scale: 1 = susceptible, 10 = resistant). Following conventional practices for the region, all seeds were treated with broad-spectrum pesticides: sedaxane, difenoconazole, mefenoxam, and thiamethoxam.

Given that our study aimed to evaluate fungal communities under realistic field conditions, the agronomic practices and field plot sizes were not standardized across sites. Instead, agronomic practices for each region were followed and local planting choices were selected based on site differences in weather, climate, and soil type (Supplementary Table S1). Full management details are provided in Supplementary Table S1. Briefly, neither Ewing nor Urbana was tilled, and maize debris was visibly present on the soil layer prior to planting. The seeding rate was 145.7 kg ha^–1 in Ewing and 160.3 kg ha^–1 in Urbana. Plot dimensions were 6.1 m by 30.5 m in Ewing and 9.1 m by 35.4 m in Urbana, with filler wheat along site edges. To allow for natural FHB establishment, fungicides were not applied during either the vegetative or reproductive stages. Conventional fertilization practices for each region were followed.

The Ewing field site was harvested on June 24, 2021 after dry-down using a Gleaner K2. Due to poor stand density across some plots caused by high soil moisture, 4.0 m was removed from the north end of all plots (0.0161 ha plot −¹). Whole plot yields were taken with a weigh wagon, while moisture content and test weight were calculated using a Dickey John mini-GAC (Auburn, IL). The Urbana field site was harvested on July 5, 2021 using an Almaco SPC 20 (Nevada, IA) with a HarvestMaster H2 Grain gauge for continuous weighing and moisture content measurements. In Urbana, two 1.5 m x 15.2 m strips from uniform sections of each plot were used to calculate yield, moisture, and test weight (0.00465 ha plot^–1).

Sample collection for mycobiome analysis targeted three tissue types (maize debris, wheat heads, wheat leaves) and five key developmental and seasonal changes in winter wheat growth, as well as key points for FHB disease pathogenesis (designated T1-T5; Figure 1b). Specifically, we sampled during: (T1) pre-winter tillering (Feekes 2); (T2) post-winter pseudo-stem erection (Feekes 4–5); (T3) heading (Feekes 10.1–10.5); (T4) flowering, when F. graminearum infection occurs (Feekes 10.5.1–10.5.3); and (T5) dough fill, when FHB symptoms present (Feekes 11.2–11.3).

For each tissue type and time point combination, either eight (leaf and debris) or four (heads) samples were purposefully collected to span the entire length and width of a field plot and then combined into one sample in sterile plastic bags (see Supplementary methods for details). Whole wheat leaves were harvested across the first four time points (T1-T4), typically targeting the second-most fully elongated leaf—to balance the need for leaf maturity with sufficient time since emergence for fungal colonization—except at T4, where the flag leaf was sampled as the ultimate leaf. Entire wheat heads were collected after emergence at the final three time points (T3–T5). For the T5 timepoint, two sets of heads were collected: heads displaying either low (<30%) or high ( > 30%) visual FHB symptoms. Lastly, maize debris was sampled at the first two time points (T1-T2), where the focus was a piece of decaying maize stem inclusive of a node. Scissors and gloves used for collection were sterilized with 70% ethanol between plots. Bagged samples were placed on ice for transport back to the lab and then frozen within 1 day of collection.

Leaf samples from each plot were randomly subsampled into 2 mm² fragments. Forty-eight fragments per plot per time point were then surface sterilized by successive washes in 95% ethanol (30 s), 0.5% NaOCl (2 min), 70% ethanol (2 min), and sterile water (30 s), placed into a sterile tube, and flash frozen. DNA was extracted from leaf fragments using the Synergy 2.0 Plant DNA kit (OPS Diagnostics LLC, Lebanon, NJ). Head samples were processed as is by lyophilizing, grinding to a powder using a tissue homogenizer (1,650 rpm for 8 min, Geno/Grinder 2010; Cole Parmer, Vernon Hills, IL), and extracting DNA from ∼20 mg of tissue using the Synergy 2.0 Plant DNA kit. Debris samples were lyophilized, ground in three passes through a Stein Mill (Steinlite Co., Atchison, KS) for 10 s each, and DNA extracted in duplicate from 20 to 40 mg of ground tissue using the DNeasy PowerSoil kit (Qiagen, Germantown, MD). Debris DNA was pooled by plot prior to PCR amplification.

DNA concentration was fluorometrically quantified (Quant-It dsDNA kit, Thermo Fisher Scientific, Waltham, MA) on a Synergy H1 microplate reader (Agilent, Santa Clara, CA) and normalized to 10 ng μL^–1. Gene amplification of fungal DNA was performed using two-stage PCR with the fITS7 (5′-GTGARTCATCGAATCTTTG-3′; Ihrmark et al., 2012) and ITS4-ngs (5′-TCCTSCGCTTATTGATATGC-3′; Tedersoo et al., 2014) primers modified with Illumina adapters (Smith and Peay, 2014). Primers were frameshifted using 3–6 random nucleotides to increase sequence complexity (Lundberg et al., 2013). Negative PCR controls were included at each stage of amplification. The first stage of PCR consisted of 12.5 μL 1 × Phusion High Fidelity PCR Master Mix (Thermo Fisher Scientific), 0.2 μM of each primer, 0.5 mg mL^–1 BSA, and 25 ng DNA template. In addition, a custom-designed peptide nucleic acid (PNA) was added at 0.75 μM final concentration to the PCR reaction to block host genomic ITS amplification (5′-CAAGTTGCGCCCGA-3′; Whitaker, 2025) with PCR-grade water up to 25 μL total reaction volume. Thermocycling conditions were 30 s at 98°C, then 25 cycles of denaturing for 10 s at 98°C, PNA annealing for 10 s at 74°C, primer annealing for 30 s at 58°C, and extension for 30 s for 72°C, and a final extension for 5 min at 72°C. PCR reactions were cleaned using AMPure XP beads (Beckman-Coulter Inc., Indianapolis, IN). Each index PCR reaction consisted of 25 μL 1 × Phusion High Fidelity PCR Master Mix, 5 μL each of two indices (Nextera XT Index Kit v2, Illumina Inc., San Diego, CA), 5 μL DNA template, and 10 μL water. All products were cleaned using AMPure XP beads, quantified in triplicate using Quant-It on a Synergy H1 Microplate reader, and PCR fragment size determined using Tapestation (D1000 ScreenTape, Agilent Technologies, Santa Clara, CA). Lastly, all samples were pooled in equimolar ratios and sequenced on an Illumina MiSeq platform (v3 chemistry, 2 × 300 bp). Raw sequencing data are provided in the NCBI SRA BioProject PRJNA1297200, accession numbers SRR44969627–SRR44969870.

At the final timepoint, an additional set of wheat heads was collected to determine disease at the plot level. Specifically, 32 bulk wheat heads were randomly collected from eight points spanning the length and width of the field plot (with 4 heads collected per point) in order to avoid bias in FHB disease state and best represent microclimate differences across the plot. For these bulk wheat heads, we sampled heads irrespective of disease state or tiller height. The bulk wheat heads were lyophilized, hand threshed to remove chaff and stems, and the grains ground to a flour-like consistency using a Stein Mill. In between samples, the mill cups were cleaned with 70% ethanol and soap, and the blade/shaft cleaned with 70% ethanol. Disease was assessed on the flour via quantification of F. graminearum biomass and DON contamination.

DNA extractions were performed in triplicate, to account for sample variation across the large volume of collected grain, using the Qiagen DNeasy Plant Pro kit (∼20 mg each). F. graminearum biomass (hereafter Fusarium biomass) was determined by quantitative polymerase chain reaction (qPCR) on each replicate DNA extraction using previously reported protocols (Kemp et al., 2020). The qPCR was performed using the Juno and BiomarkHD instruments (Fluidigm, San Francisco, CA). The LinRegPCR program (v.1.0; Ruijter et al., 2009) was used to correct raw fluorescence values for variable amplification efficiencies by first determining a baseline, then estimating amplification efficiency for each sample individually, and computing a starting quantity. An arithmetic mean was calculated from four technical replicates per assay and then the geometric mean was computed from the three primer assays per organism (i.e., F. graminearum and wheat, primers and probes listed in Supplementary Table S2). To minimize error associated with variation in DNA extraction efficiencies, Fusarium biomass was expressed as the ratio of Fusarium DNA to wheat DNA. A final arithmetic mean was calculated on the three DNA extraction replicates for each field plot as an average. Lastly, the coefficient of variation was calculated and assessed as a measure of accuracy (Supplementary Figure S1).

Toxin analysis was performed in triplicate. Ground grain material (∼0.50 g) was extracted with 10 mL of acetonitrile-water (86:14 vol/vol) in a 50 mL Falcon tube with shaking for 15 min. Following centrifugation, a 5 mL portion of each extract was purified through a small column packed with 1.5 g of a mixture of silica gel C18, alumina, and activated charcoal (150:450:2.5). Then, a 2 mL aliquot of each purified extract was dried under a stream of air. Trimethylsilyl (TMS) derivatives were prepared as described previously (Whitaker et al., 2023). TMS derivatives of purified DON (0.3125–80 μg) were used to construct a standard curve for quantitation.

GC-MS analyses were performed on an Agilent 7890 gas chromatograph fitted with a HP-5MS column (30 m, 0.25 mm, 0.25 μm) and a 5977 mass detector as described previously (Whitaker et al., 2023). Similarly to the Fusarium biomass results, the coefficient of variation was assessed as a measure of accuracy (Supplementary Figure S1).

Sequence data were processed using DADA2 (v.1.26.0) (Callahan et al., 2016) in R (v.4.2.2) (R Core Team, 2024) to determine amplicon sequence variants (ASVs). Primers were removed using cutadapt (v.1.18) (Martin, 2011) and default parameters were used for all other steps. After filtering, denoising, merging, and chimera check, the raw paired read count of experimental samples was reduced from 12,657,525 to 7,998,719. Non-fungal ASVs, including plant- and algal-origin and chimeras not screened by DADA2, represented < 0.5% of the remaining data and were removed in two steps. First, by using a BLAST search (-word_size 50, -perc_identity 60) to a manually curated database of full length ITS sequences for wheat, maize, and barley and the UNITE database (Kõljalg et al., 2013). Second, for sequences with neither high confidence matches to plants nor fungi, using a BLAST search to the NCBI nr database (Clark et al., 2016). Analysis of the negative controls revealed that two ASVs were present in two plates of the first round of PCR, one of which was unique across the dataset and the other which was not. No contaminants were present in the second round of PCR. Thus, we simply removed both ASVs from further analyses, resulting in 2580 ASVs. Putative taxonomy was assigned using the QIIME2 implementation (Bolyen et al., 2019) of a naïve Bayes classifier (classify-sklearn) (Bokulich et al., 2018) with the UNITE fungal ITS database (v10, Euk, RefS, dynamic) (Abarenkov et al., 2025). All taxonomic assignments were manually curated to replace instances of “incertae sedis” with “NA” that were otherwise unidentified at the genus level. Fungal richness (vegan v.2.6–4) (Oksanen et al., 2022) was calculated after quality filtering.

Statistical analyses were performed in R (v.4.4.1) and made substantial use of the DADA2 (v.1.28.0), phyloseq (v.1.42.0) (McMurdie and Holmes, 2014), and DESeq2 (v.1.38.3) (Love et al., 2014) packages. All data figures were created using ggplot2 (v.3.4.1) (Wickham, 2016). The bioinformatics and statistical scripts, as well as experimental data, have been archived on Zenodo (Whitaker and Gdanetz, 2025).

Yield was 17.0 kg ha^–1 higher on average in Urbana compared to Ewing (P < 0.001), but there were no significant varietal differences (P = 0.29, Figure 1c; Supplementary Table S3). By contrast, DON and Fusarium biomass varied by variety (P = 0.0401, P = 0.0269, respectively) at the Ewing site, but not at the Urbana site (Figure 1c; Supplementary Table S3). The highest DON contamination ( > 1ppm) was detected at the Ewing site in the two varieties with resistance ratings of 5 and 8 (moderately susceptible and moderately resistant, respectively). Neither DON nor Fusarium biomass predicted yield (P = 0.40 and P = 0.33, respectively).

In the wheat phyllosphere and maize debris mycobiome, we identified 2,580 ASVs across the growing season. Most ASVs belonged to the Ascomycota and Basidiomycota phyla (66.8 and 29.9% respectively). The genera that accounted for the most sequencing reads included Cladosporium (18.5%), Fusarium (12.3%), Vishniacozyma (8.5%), Alternaria (8.1%), Hannaella (5.0%), and Sporobolomyces (4.9%). The top ten most abundant ASVs varied in their relative abundance across treatments (Supplementary Figure S2) and accounted for 53.0% of all sequencing reads. They included two Cladosporium, as well as a Fusarium, Alternaria, Vishniacozyma, Colletotrichum, Sporobolomyces, Bullera, Papilotrema, and one taxon unnamed at the genus level. The majority of ASVs were unique to a single tissue type (84.5%; 2,179 ASVs), though some were still shared across all three tissue types (3.8%; 99 ASVs; Supplementary Figure S3).

Tissue type was a significant predictor of fungal richness (P = 0.001; Figure 2a), with leaf communities displaying the least richness relative to debris and head tissues. Similarly, mycobiomes were strongly structured by tissue type (P = 0.001; Figure 2b; Supplementary Table S3), with a clear visual separation in community structure along PCoA axes 1 and 2 by tissue type superseding differences between the two sites (Figure 2b).

Nevertheless, additional analyses by tissue type revealed that site, variety, and timepoint significantly affected community structure for debris and head tissues (significant 3-way interactions, both P = 0.001; Figures 2c,e; Supplementary Table S4). For leaf tissues, community structure was affected by significant two-way interactions between site and variety (P = 0.006), as well as between site and timepoint (P = 0.001), but there was no significant three-way interaction between site, variety, and timepoint (P = 0.274; Figure 2d; Supplementary Table S4). Visual assessment revealed differences in fungal community structure between the two sites in all tissue types (Figures 2c–e), as well as a clear temporal shift in the community structure results for debris and heads (Figures 2c,e), but a less clear temporal shift for leaf communities (Figure 2d). For the heads collected at the final timepoint, site, variety, and disease significantly impacted community structure (significant 3-way interaction P = 0.001; Figure 2f; Supplementary Table S4), with moderate shifts in the fungal communities visible in the PCoA (Figure 2f).

Variance partitioning analysis revealed that the importance of treatments in explaining variation in the mycobiome differed across tissue types. Specifically, site and timepoint explained nearly equal conditional variation in fungal community structure in leaf and head tissues (7.4–10.8%; Figure 3a). However, site explained 2.5 times as much variation in community structure as timepoint for debris (26.4% vs. 9.8%, respectively; Figure 3a). Host variety only significantly explained variation in community structure in wheat heads (1.2%; Figures 3a–c), but not in debris or leaves (0.3–0.6%; Figure 3a). At both sites, varietal differences between wheat head communities tended to be stronger earlier in development (e.g., T3), than later (e.g., T5; Figure 3b,c). For the heads collected at the final timepoint (T5), disease and site significantly explained fungal community variation, with site explaining 5x as much (2.7% vs. 13.3%, respectively; Figure 3a). Across all tissue types, unexplained variation in fungal community structure predominated (66–85%; Figure 3a).

The source-sink analyses of tissue-to-tissue fungal transfer revealed that debris was a major source of fungi to wheat leaves during early development (i.e., T1), but not later (i.e., T2; Supplementary Tables S5, S6). Specifically, the estimated source proportion from debris to leaves was 0.70 on average before the over-wintering period compared to 0.20 after (Figure 4a). On the other hand, leaves were consistently a low source of fungi to wheat heads, both during heading and at flowering (T3 and T4; Supplementary Tables S5, S6), with average estimated source proportions of 0.10 and 0.20, respectively (Figure 4b).

Next, we used source-sink analyses to identify how much the preceding timepoint was a source to subsequent timepoints, for each tissue type. The preceding timepoint was a high source of fungi to the subsequent timepoint for both debris and head tissues, regardless of the preceding timepoint tested as a source or disease status of the heads tested as a sink (average source proportion = 0.78–0.91; Supplementary Tables S5, S6; Figure 4c). However, estimated marginal means analysis revealed that high disease heads during dough filling (T5) were a significantly greater recipient of fungi from flowering heads (T4), than were low disease heads (0.91 vs. 0.85 on average; Supplementary Table S6; Figure 4c). In other words, low-disease heads were colonized by a greater proportion of fungi from an unknown source and slightly lower proportions of fungi that had colonized during flowering. By contrast, early timepoint leaves were not a major source of fungi to later timepoint leaves (Supplementary Tables S5, S6), with average source proportion varying from 0.29 to 0.45 across the three timepoints tested (Figure 4c).

Given the importance of the FHB-pathosystem in wheat, we used DESeq2 to identify fungal taxa that were differentially enriched under high or low disease conditions (T5 timepoint). Only 14 ASVs were significantly differentially abundant according to disease state, after accounting for site and variety differences (Supplementary Tables S7; Figure 5a). Of these, most (10 ASVs) were enriched under low disease, with the remaining (4 ASVs) enriched under high disease. The taxa enriched in high disease heads all belonged to the Ascomycota, including two Fusarium (ASV2, ASV278), a Pseudopithomyces (ASV122), and a Cladosporium (ASV151). The identification of two enriched fusaria ASVs corroborated our visual selection for high disease symptoms in the field. However, some Fusarium infection was still identified in low disease heads (mean relative abundance of all Fusarium ASVs in low vs. high disease: 6.5% vs. 30.2%). Of the taxa enriched for low disease, most were basidiomycetous yeasts (8 ASVs) and included two Filobasidium (ASV11, ASV155), two Sporobolomyces (ASV27, ASV96), two Tilletiopsis (ASV31, ASV99), one Entyloma (ASV331), and one unknown Tremellales (ASV201). Two Ascomycota (ASV1:Cladosporium and ASV29:Zymoseptoria) were also enriched in low disease heads. Overall, more ASVs were unique to the low disease heads (42.6%) (Figure 5b) than were unique to the high disease heads (22.8%) (Figure 5b); though many ASVs were still shared between the two disease states (34.6%) (Figure 5b).

Statistical comparison of co-occurrence networks indicated expected complexity based on sample type, with maize debris displaying highest complexity (Table 1, Supplementary Table S8 and Supplementary Figure S4), followed by the wheat head networks (Table 1, Supplementary Table S9, and Figure 6). As observed in other pathosystems, high-disease wheat head networks had reduced complexity relative to low disease head networks, as indicated by fewer nodes, edges, and hub taxa, as well as smaller average edge weight—suggesting reduced community stability (Figures 6a,b and Tables 1, 2).

Four hub taxa (ASV69:Bullera, ASV70:Symmetrospora, ASV371:Udeniomyces, ASV494:Entyloma) were shared between networks of the two disease states (Table 2). Three hub taxa each were shared between the maize debris and high-disease (ASV4:Alternaria, ASV12:Hannaella, ASV100:Cercospora) or between the maize debris and low-disease (ASV28:Papiliotrema, ASV37:Didymellaceae, ASV141:Rhodotorula) networks, respectively. No hub taxa were shared across all three networks (Table 2). The DESeq2 analysis identified 14 ASVs as significantly different between disease status of heads; four of these were also identified as hub taxa in the high disease network (including ASV2:Fusarium) and two in the debris network (Table 2).

The greatest number of Fusarium spp. nodes (26 ASVs) were present in the maize debris network (Supplementary Table S8), with six Fusarium nodes shared with at least one of the head networks (Supplementary Table S10). Four Fusarium spp. nodes were present in both wheat head networks and two Fusarium sp. nodes (ASV2, ASV278) were shared across all three networks (Supplementary Tables S9, S10). There were similar numbers of unique ASVs negatively associated with Fusarium spp. nodes in both head networks, with taxonomic assignments to Basidiomycota (ASV331:Entyloma, ASV94Sporobolomyces, ASV144: Sporobolomyces, ASV93:Symmetrospora, ASV576:Symmetrospora) and Ascomycota (ASV118:Parastagonospora). Taxa that positively co-occurred with Fusarium nodes in both disease states included: Alternaria (ASV431, ASV1154), Ophiosphaerella (ASV166, ASV294), Papiliotrema (ASV88, ASV106, ASV145, ASV266, ASV362), Phaeosphaeriaceae (ASV167), Sporobolomyces (ASV96, ASV144), and Vishniacozyma (ASV74, ASV87, ASV289, ASV410). Phylloplane yeasts; Vishniacozyma, Aureobasidium, Filobasidium, Hannaella, and Occultifur were frequently associated with Fusarium spp. in maize debris networks and rarely found in head networks (Supplementary Table S11). When these genera were detected in heads, it was in low-disease status samples, mirroring the DESeq2 analysis.