We sequenced the genome of AF_NPK13tox to understand its relationship to AO_RIB40 and AF_AflaGuard. A. flavus 13tox was grown on potato dextrose agar (PDA) at 30 °C for 48 hours. DNA was extracted directly from the spores using the method outlined by (Zhao et al., 2019). DNA concentrations were measured for each extraction using a Qubit fluorometer. PCR-free 150-bp paired-end libraries were prepared and sequenced by Novogene (https://en.novogene.com/) on an Illumina NovaSeq 6000 platform. The raw whole-genome sequencing data for A. flavus 13tox can be accessed through the BioProject accession number PRJNA1277813. Raw paired-end fastq files were adapter and quality trimmed using trim_galore (Martin, 2011) with the following parameters: “-q 30”, “--stringency 1”, and “--length 80”. Adapter and quality trimmed reads were assembled using SPAdes (Bankevich et al., 2012) using the “--careful” option and K-mer sizes of 55, 75, 95, 105, and 115.

Eleven genome assemblies spanning the diversity of A. oryzae and A. flavus populations were obtained from NCBI and analyzed with the AF_NPK13tox assembly. These assemblies were A. flavus AflaGuard (GCA_012896875.1), OK-A-8-1-S (GCA_024676875.1), NRRL 3357 (GCF_009017415), A. oryzae TK14 (GCA_009685055.1), TK31 (GCA_009686725.1), TK49 (GCA_009687085.1), TK56 (GCA_009687225.1), TK57 (GCA_009687245.1), TK60 (GCA_009687305.1), RIB40 (GCF_000184455.2), 14160 (GCA_019097755.1) (Machida et al., 2005; Weaver et al., 2017; Watarai et al., 2019; Drott et al., 2020; Chacon-Vargas et al., 2021), and Aspergillus parasiticus CBS-117618 (GCA_009176385.1). For this analysis, gene models for all 11 genomes were predicted using Augustus v3.4.0 with the following parameters: “--species=aspergillus_oryzae,” “--strand=both,” “--gff3=on,” “--uniqueGeneId=true,” and “--protein=on” (Stanke and Waack, 2003). We used Augustus for gene prediction even for genomes with existing annotations to avoid bias from differences in gene model predictions. Single-copy orthologs across the 13 proteomes were identified using OrthoFinder v2.5.5 (Emms and Kelly, 2015; Emms and Kelly, 2019). The identified single copy ortholog proteins were concatenated and aligned using MAFFT v7.481 (Katoh and Standley, 2013). A phylogenetic tree was then constructed from the concatenated protein alignment with IQ-TREE v2.1.3 using 1,000 bootstrap replicates (-B 1000) and the LG amino acid replacement model (Minh et al., 2020).

AO_RIB40, AF_NPK13tox, and AF_AflaGuard were selected for rice fermentation and subsequent analysis because these strains represent distinct populations within the A. oryzae/A. flavus species complex (Figure 1). AO_RIB40 was obtained from ATCC, AF_AflaGuard was obtained from the USDA SRRC culture collection, and AF_NPK13tox was obtained from Nancy Keller’s lab at the University of Wisconsin Madison.

We used Lundberg (Richvale, CA) Organic California Sushi Rice, which has a starch content ~75%, for all experiments involving rice fermentation. Freezer stocks of conidia for each strain were rapidly thawed and cultured on steamed rice 7 days at 30°C in petri dishes, at which time, the plates were flooded with 50% water + 49.9% glycerol + 0.1% tween to harvest conidia. Conidia were then quantified on a haemocytometer and normalized to 1x10⁶ conidia per mL. Rice was prepared for fermentation by autoclaving short grain white rice at 170% hydration for 15 minutes at 121°C. 12 g of sterilized short grain rice was inoculated with 1 mL of normalized conidia stocks (~1 million conidia) in 60 mm petri dishes sealed with parafilm and incubated for 48 hours at 30°C. Fermentations for amylase-activity and e-tongue were performed in triplicate, while fermentations performed for DH-GC-MS were performed in duplicate. However, we only analyzed a single unfermented rice (negative control) sample for e-tongue and DH-GC-MS analyses. Thus, control values serve as a baseline reference rather than a replicated measure of variability. As a negative control we also performed the same experiment but added 1 mL of the conidia harvesting solution (50% water + 49.9% glycerol + 0.1% tween) without conidia (rice only).

At the end of the 48 hours period, samples for e-tongue and DH-GC-MS were transferred to 50 mL conical tubes, 25 mL of water was added, and samples were vortexed for 60 s. A 48-hour incubation was selected because it reflects the typical duration used in standard koji production, during which A. oryzae reaches peak amylolytic activity and produce characteristic early-stage flavor compounds. Water was incorporated at the end of fermentation solely to enable e-tongue measurements, which require liquid-phase samples. Samples were flash-frozen in liquid nitrogen immediately after collection and shipped on dry ice via overnight express to Medallion Labs (Minneapolis, MN). Transit time was less than 12 hours, and samples were stored in a freezer immediately upon receipt. An unfractionated aliquot of the fermented rice/water slurry was used for E-tongue and DH-GC-MS analysis, minimizing the potential for volatile loss. Our overall experimental design is depicted in Figure 1.

We measured alpha amylase activity during rice fermentation across the three fermentation starters, because alpha amylase activity is essential for rice fermentation and is one of the defining features of A. oryzae domestication (Hunter et al., 2011; Gibbons et al., 2012). To measure alpha amylase activity, 12.25 g of fermented rice samples were combined with 10 mL of molecular water, vortexed for 1 minute, and a 2 mL aliquot was withdrawn with analysis. The 2 mL aliquot was fractionated by centrifugation at 1,000 g for 10 min. A 1 mL aliquot of the resulting supernatant was withdrawn and used for alpha amylase quantification. Alpha amylase activity of each strain was evaluated with Ceralpha α-Amylase Assay Kit (cat. No. K-CERA) in biological triplicate. This method utilizes non-reducing end-blocked p-nitrophenyl maltoheptaoside in the presence of an excess of thermostable α-glucosidase. The supernatant was diluted 1:24 with the assay buffer before proceeding with the assay protocol. Absorbance values (405 nm) were recorded with a SpectraMax i3 microplate reader (Molecular Devices, USA). The negative control for this assay was developed by combining all reagents in the same proportions as test samples, but without incubation such that the active enzyme was denatured and thus no longer capable of hydrolyzing starch. This was to confirm that any changes in absorbance (405 nm) were due to enzymatic activity rather than non-enzymatic activity or baseline fluctuations.

Aflatoxin quantification was carried out using a modified protocol adapted from Alshannaq et al. (2018) and Acevedo et al. (2025) (Alshannaq et al., 2018; Acevedo et al., 2025). Isolates were grown in slant cultures prepared with 2 mL of potato dextrose (PD) medium in 10 mL glass test tubes. Approximately 10⁵ conidia were inoculated into each tube using a sterile loop. The tubes were placed at a 45° angle on a rack and incubated at 30 °C for 7 days to allow for mycelial growth and toxin production. After incubation, aflatoxins were extracted using an organic solvent method. One mL of chloroform was added directly to each culture tube, thoroughly mixed, and then centrifuged at 5,000 × g to facilitate phase separation. From the chloroform phase, 0.5 mL was carefully collected and left to evaporate at room temperature. The resulting residue was reconstituted in 0.5 mL of High-performance liquid chromatography (HPLC) mobile phase consisting of water, methanol, and acetonitrile (50:40:10, v/v/v). Samples were filtered through a 0.45 μm membrane filter prior to HPLC analysis.

HPLC was used to quantify aflatoxins AFB1, AFB2, AFG1, and AFG2. Analyses were performed on an Agilent 1100 series system equipped with a degasser, autosampler, quaternary pump, and a 1260 Infinity diode array detector (Agilent Technologies, CA, USA). Separation was achieved using a Zorbax Eclipse XDB-C18 column (4.6 mm × 150 mm, 3.5 μm particle size) with a flow rate of 0.8 mL/min. Aflatoxins were detected at a wavelength of 365 nm.

Traditional descriptive sensory testing with human subjects was not possible because AF_NPK13tox is an aflatoxin producer. Thus, we analyzed our samples via e-tongue, which measures flavor attributes by using an array of chemical sensors combined with pattern recognition algorithms to detect and differentiate flavor-related compounds in a liquid sample. E-tongue analysis assessed these compounds and provided an interpretation of how these might impact eight different flavor profiles (i.e., umami, aftertaste richness of umami, saltiness, sourness, bitterness, aftertaste of bitterness, astringency, and aftertaste of astringency). These attributes are based on a relative comparison to internally developed references. Data were reported as relative comparisons to the control, using a scale where 1 unit represents a level of differential sensitivity equivalent to an estimated 20% difference in taste. Three biological replicates were analyzed for each fermentation, and one sample was analyzed for the unfermented rice negative control to establish a baseline. For each flavor attribute, we conducted an ANOVA between the fermentation starter strains assuming a p-value < 0.05 cutoff. For comparisons in which p-values < 0.05, we conducted post hoc Tukey Kramer HSD tests to compare pairwise differences between samples, imposing a p-value < 0.05 cutoff. Statistical analysis was conducted in JMP v18.2.0.

All sample processing, compound identification and quantification were performed at Medallion Labs (Minneapolis, MN). Samples were analyzed as received by dynamic headspace on an Agilent 5977 MSD coupled to an Agilent 7890A gas chromatography unit with a Gerstel multipurpose sampler (Agilent, USA). Compounds were identified and matched by a proprietary spectrum library, with no other standards or calibrations. Quantitation of compounds present was determined by comparison of their response to that of an internal standard added to each sample. Prior to testing, Parts Per Million (PPM) values were calculated relative to the total sample weight as received (wet basis) (as samples consisted of frozen fermented rice/water slurry). Because each fermentation condition was represented by only two biological replicates (and the unfermented control by one), we did not perform inferential statistical tests on the DH-GC-MS volatile compound data. Instead, for each detected compound we summarized relative abundance across replicates using descriptive statistics (mean and range) to characterize strain-dependent patterns.

We analyzed the phylogenetic relationship of the three rice fermentation starters (AO_RIB40, AF_AflaGuard and AF_NPK13tox) in relation to 7 additional A. oryzae genomes, 2 additional A. flavus genomes that are representative of the major populations based on previously published work (Gibbons et al., 2012; Watarai et al., 2019; Drott et al., 2020; Chacon-Vargas et al., 2021; Kim et al., 2025), and a strain of A. parasiticus which was used as an outgroup. The genome of AF_NPK13tox was sequenced and assembled for the current study (BioProject ID PRJNA1277813), while all other genome assemblies were publicly available (see methods for accession numbers). We used OrthoFinder to identify single copy orthologs (Emms and Kelly, 2015; Emms and Kelly, 2019), and then built a maximum likelihood phylogenetic tree based on a protein alignment of 8,092 single copy orthologs (which included 4,445,407 amino-acid sites and 48,990 parsimony informative sites). The phylogenetic tree showed three phylogenetically distinct clades of A. oryzae, and two clades of A. flavus (Figure 1), consistent with previous studies. Importantly, the three fermentation strains resided in distinct clades.

We measured alpha-amylase activity at the end of our fermentation because variability in alpha-amylase production could impact volatile compound formation (e.g. higher sugar concentrations could lead to more ethanol production by yeasts which in turn could increase the pool of substrates for ester synthesis), and taste perception (e.g. higher alpha-amylase activity would increase sugar content and sweetness of the fermentation product) (Yoshizawa, 1999). Unsurprisingly, alpha-amylase activity was significantly different between isolates (Figure 2) (ANOVA, p-value = 5.2e-6) and significantly higher in AO_RIB40 compared to the A. flavus isolates (Tukey-Kramer: AO_RIB40 vs. AF_NPK13tox: p-value = 1.0e-4 and AO_RIB40 vs. AF_AflaGuard: p-value = 8.6e-6).

We quantified the production of aflatoxins (AFB₁, AFB₂, AFG₁, and AFG₂) by AO_RIB40, AF_NPK13tox, and AF_AflaGuard under culture conditions known to promote aflatoxin biosynthesis (Alshannaq et al., 2018; Acevedo et al., 2025). Aflatoxins were not detected in cultures of AO_RIB40 or AF_AflaGuard. In contrast, AF_NPK13tox produced approximately 97 ppm AFB₁ and 34 ppm AFG₁. These results are consistent with prior knowledge: AF_NPK13tox is closely related to the aflatoxigenic strain A. flavus NRRL 3357 (Figure 1), whereas AO_RIB40 and AF_AflaGuard harbor well-documented loss-of-function mutations in the aflatoxin biosynthetic gene cluster (Gibbons et al., 2012; Alshannaq et al., 2018).

Fermented rice water slurries underwent e-tongue analysis to measure and compare flavor attributes between starter cultures. Compounds known to associate with perceived saltiness, sourness, umami, umami aftertaste, astringency, aftertaste-A (astringency aftertaste), bitterness and aftertaste-B (bitter aftertaste) were measured for the unfermented rice negative control, and the biological replicates of each fermentation starter culture. Four of the 8 comparisons were statistically significant (ANOVA, all p-values ≤ 0.002, Figure 3). Astringency is a tactile sensation perceived as dryness, roughness, and a tightening feeling in the mouth, often accompanied by a slight contraction of the tongue and oral tissues, and it is often accompanied with acidity or bitterness sensations (Pires et al., 2020). Interestingly, astringency and aftertaste-A was significantly lower in the AO_RIB40 fermentation than in the A. flavus fermentations (Figure 3). Additionally, bitterness was substantially lower in all samples relative to the rice control and significantly lower in the AO_RIB40 fermentation compared to the A. flavus fermentations which were themselves significantly different (Figure 3). Lastly, umami aftertaste, described as the richness, complexity, and full bodied, was significantly higher in the AO_RIB40 fermentation compared to the A. flavus fermentations which, were also significantly different from each other (Figure 3).

We performed DH-GC-MS of the fermented rice water slurries to identify unique compounds and compare the relative abundances of compounds between starter cultures. We detected 34, 35, 48, and 61 unique compounds in the rice slurry without a starter culture, the AF_NPK13tox fermentation, the AF_AflaGuard fermentation, and the AO_RIB40 fermentation, respectively (Figure 4, Table 1). Between the starter culture samples, 30 compounds were detected in all samples, while 1, 4 and 19 compounds were uniquely detected in the AF_NPK13tox, AF_AflaGuard, and the AO_RIB40 fermentations, respectively (Figure 4).

To understand how samples volatile profiles were related, we performed a PCA on each sample and replicate using the compound abundance (in parts per million) matrix as input. PC1 explained 62.5% of variance, and separated rice, AF_NPK13tox, and AF_AflaGuard, by AO_RIB40, while PC2 explained 14.9% of variance and separated the AF strains and the rice sample (Figure 4).

To compare compound abundance patterns between fermentations, we summarized the two biological replicates for each starter culture using descriptive statistics, rather than inferential tests, because replication was limited (n = 2 per fermentation; n = 1 control). Several clear trends emerged. Across the panel of detected compounds, AO_RIB40 consistently produced substantially higher levels of many alcohols, aldehydes, ethers, esters, and ketones compared with the A. flavus starter cultures (Table 1). For example, some compounds of note include 2,3-butanediol and R-2,3-butanediol, which was uniquely detected in AO_RIB40 and contributes to the floral and fruity aroma of soy sauce, 3-methyl-Butanal, which is an aldehyde detected in moromi fermentation (Li et al., 2023), and 2-methyl-butanal, which has a cocoa or coffee-like aroma and has been detected in fermented sausages (Olivares et al., 2009) (Table 1, Figure 4). Conversely, a small number of compounds were detected in both A. flavus fermentations but not in AO_RIB40. Interestingly, one such compound, ethyl decanoate, has been linked to both sweet and vinegar aromas (Xiao et al., 2016), and was not detected in the AO_RIB40 fermentations.