Large-scale multi-omics profiling reveals environmental and evolutionary drivers of fungal phylogeographic and metabolic diversity

Domain: biology
Published: 2026-03-18
Editorial outcome: major_revision

Abstract
Chemical innovation is essential for fungi to adapt to ever-changing ecological environments. However, the environmental and evolutionary drivers of fungal metabolic differentiation remain ambiguous. Here, we show the phylogeographic diversity of 1,052 Aspergillus flavus strains across four continents, as conducted through phylogenetic and biogeographical analysis, including 544 newly sequenced strains from China. These strains exhibit varying levels of population-specific mycotoxin production, as determined by population metabolomics analysis. We report a toxigenic subpopulation from China, identified through comparative population genomics analysis. Pan-metabolome analysis reveals strong phylogeographic metabolic patterns associated with specific ecological niches. Low-mycotoxin production clades harbour distinct uncharacterized biosynthetic gene clusters and produce different specialized metabolites instead. This discrepancy is only partially explained by variation in biosynthetic pathway genes, and changes in regulation and primary metabolism appear to mainly drive differentiation of specialized metabolite profiles across fungal populations, as indicated by pangenome profiling, metabolites-genome-wide association study, genotype-environment association study, pan-transcriptome analysis, and gene knockout experiments. Altogether, our results reveal how environmental shifts drive the fungal metabolic evolution, and provide insights for predicting the risk of harmful fungal outbreaks and for biogeographically-informed, precise control measures.

Introduction
Fungal specialized metabolism has a profound impact on the Earth's ecosystem1, especially, since the world is currently experiencing rapid anthropomorphic climate change, which is increasing the risk of harmful fungal outbreaks that infect human, destroy crops or poison foods with their toxic specialized metabolites2. Climate shifts promote fungal evolution, altering host-pathogen interactions, and contributing to the emergence of new pathogenic strains3. Fungi are well-known for their capacity to produce a wide range of bioactive compounds. These molecules possess important ecological functions in inter-microbial competition, defense, signalling, nutrient acquisition, and development4. Many fungi secrete mycotoxins to facilitate pathogenesis or to defend a specific ecological niche5. Unfortunately, the mycotoxins, such as aflatoxin B1, mainly produced by the A. flavus, is classified as a Class I carcinogen by the World Health Organization, may leads to poisoning events in humans due to their carcinogenicity. According to FAO data, ~25% of the world's food is contaminated yearly by mycotoxins, causing large socioeconomic losses globally and posing a major challenge to food security and public health6. Moreover, A. flavus is also the second most common opportunistic pathogen causing aspergillosis. More than 2.113 million people suffer from aspergillosis every year worldwide, with a mortality rate of 85.2%7.

An improved understanding of phylogeographic diversity of harmful fungal specialized metabolism would inform the design of precision control strategies to mitigate the growing threat of fungal infection and mycotoxins. The extraordinary diversity of fungal specialized metabolites is maintained by evolutionary adaptations to challenges posed by complex environments8. Specialized metabolite repertoires have been shown to evolve rapidly through the functional divergence, horizontal transfer, and de novo assembly of biosynthetic gene clusters (BGCs) that encode for their production5. Recent computational genome mining efforts, including global analysis of 1,000 fungal genomes9, investigation of inter- and intraspecies variation in Aspergillus section Nigri10, and a comparative genomic study of 23 Aspergillus species from section Flavi11 have shown that fungal biosynthetic diversity is vastly untapped. Even within individual fungal species, considerable variation in BGC repertoires exists12, and, sometimes, closely related pathogenic and non-pathogenic strains show distinct differences in their biosynthetic capabilities13. However, little is understood about the processes that shape the biogeographical diversity of fungal specialized metabolism within species under climatic shifts. The hypothesis is that the environment selects for diversification; yet, to what extent the environment really does select for specific specialized metabolic repertoires and their evolution direction, and how this is encoded genetically, has remained unclear. While recent key studies showed the existence of considerable pan-genomic diversity of accessory BGCs12, 14 and population-specific differences in metabolite production15. Traditionally, BGC variation is considered to be the main driving force; however, other drivers, such as regulatory, primary metabolism and environmental selection of gene variation, have not been systematically studied. Therefore, a transcontinental phylogeographic analysis is required to systematically identify the major genetic and environmental drivers that lead to different metabolic profiles across populations around the globe.

Here, we select the ubiquitous species of broad societal impact, A. flavus, as a model system to analyze the evolutionary drivers of phylogeographic diversity of fungal specialized metabolism. We perform de novo genome sequencing of >550 representative isolates from China and also investigate these strains using untargeted eco-metabolomics; the genomics data were then analyzed together with 508 previously published environmental (n=412) and clinical (n=96) A. flavus genomes from nine other countries. We report an aflatoxigenic subpopulations that mainly reside in southern and central China. Strains from different populations exhibit conserved BGC repertoires, yet show phylogeographically distinct specialized metabolite production. Notably, low-aflatoxin clades tend to harbour clade-specific unknown BGCs and were found to produce different metabolites instead, including other mycotoxins, which may prompt a rethinking of food safety management strategies. To our surprise, these differences in metabolic output are only partially explained by polymorphisms of BGC genes, but mainly driven by evolutionary rewiring of key transcriptional regulation and primary metabolism. Genotype-environment association (GEA) analysis further indicates that environmental local adaptation promotes these processes. Altogether, our multi-omics approach reveals evolutionary drivers of the phylogenetic diversity of fungal specialized metabolites in the context of a changing environment.

Results
Multi-omics dataset across 1,052 A. flavus strains
To facilitate dissecting the nature and origins of fungal specialized metabolic differentiation from a comprehensive overview of variation at the species level, 544 representative A. flavus and five Aspergillus section Flavi strains were selected from an in-house strain library from China. The chosen sampling sites represent six climatic ecological zones to maximize their ecological breadth and geographical origins (Figure 1a, Supplementary Figure S1a, and see for detailed strain sampling site coordinates in Supplementary Table 1). The analysis workflow of this study is shown in Supplementary Figure S1b. Subsequently, these isolates were deep-sequenced using Illumina paired-end sequencing with a 97-fold mean sequencing depth and 99% mean mapping rate (Supplementary Figure S2a). The reads of all isolates were de novo assembled for obtaining the whole genome sequence. The average genome size was 37.9 ± 2.2 Mb, with 47.5 ± 0.3% GC content and 2.6 ± 0.7% repetitive sequence content. In addition, 508 previously sequenced A. flavus and two Aspergillus section Flavi strains genomes were downloaded from NCBI GenBank (Supplementary Figure S1a and Supplementary Table 1,2), mainly from three studies: 95 of them originate from a study of USA populations15, 225 isolates from infected plant parts of corn or soil samples16, and 96 representive clinical samples mainly originating from a study by Hatmaker et al.17. The scaffold number of newly sequenced assemblies (median value = 109.5 ± 91.2) is significantly lower than the number (median value = 317 ± 653.9) of draft A.flavus genomes that have been published so far (Supplementary Figure S2b and Supplementary Table 3), indicating that our genome assemblies are of high quality. We further selected 27 representative strains from a species phylogeny for PacBio long-read sequencing to obtain chromosome-level assemblies for the construction of high-quality pangenomes. The scaffold number (median value = 20.5 ± 4.9) of the 27 PacBio-sequenced strains was significantly lower than the Illumina assembly results (median value =109.5 ± 91.2), near chromosome-level (Supplementary Figure S2c & 2d). The HiC sequencing technology was further used to sequence 8 strains, which assisted in the assembly to the chromosome level. Most single nucleotide polymorphisms (SNPs) are present at low frequencies, 92% with a minor allele frequency (MAF) < 0.3 (Supplementary Figure S3a). Extended results on detailed genetic differentiation and diversity characterization, such as InDels (Insertion/Deletion, InDel), structural variation (SVs), and transposon variation in the population are available in the Supplementary Information. Untargeted metabolomics profiles for >550 Aspergillus sp. isolates were simultaneously acquired by ultra-performance liquid chromatography coupled with Orbitrap Fusion high-resolution mass spectrometry using data-dependent acquisition (DDA) of mass fragmentation spectra (UPLC-HRMS/MS) (see Methods). All metabolomes were annotated in detail using a range of computational metabolomics procedures that utilize spectral alignment, spectral similarity networking, and de novo structure prediction(see Methods). The raw metabolomic data of 95 strains from the USA generated by Drott et al.14 were downloaded from GNPS. In addition, 28 strains (with three biological replicates each and a total of 84 samples) from different clades with different phenotypes were selected to obtain species-level pan-transcriptome variation profiles. Finally, to conduct genotype-environment association (GEA) analysis, we collect the data of 21 climate variables and two geographical factors from the China weather data website (http://data.stats.gov.cn) and soil survey data of 24 soil physiological metrics from our lab.

Discovery of an Aflatoxigenic Subpopulations from China
To provide a strong foundation for transcontinental-scale phylogeographic analysis, we reconstruct a high-resolution phylogenetic tree from all 1,059 genomes (including seven outgroups) using a genome-wide SNP data matrix (including 387,282 SNPs). This facilitate the study of phylogenetic relationships of strains isolated from different geographical origins and niches (Figure 1a and Supplementary Table 1). The phylogenetic analysis of A. flavus on a transcontinental scale is genetically differentiated into eight subpopulations. W-clade 6 was discovered and contain high-frequency isolates with aflatoxin-producing capacity(APC) sampled mainly from China (Figure 1a). Previously, Geiser et al.18 constructed a phylogenetic tree of USA A. flavus isolates using three marker genes and divided the population into two groups (I, II). Drott et al.15 recently reported that the USA A. flavus population is subdivided into three genetically differentiated subpopulations (Pops) (Pop A, B, and C) at genome-wide resolution. The study of Hatmaker et al.17, recently inferred that the A. flavus population could be grouped into five populations: Pop A, B, C, D, and S-type, which were also identified in our phylogeny. However, our study at the largest scale to date revealed that of which Pop A was further differentiated into four subpopulations (Pop A1, A2, A3, A4) (Figure 1a and Supplementary Figure S4). Specifically, Pop B is synonymous with the W-clade 2 found in this study, and Pop C (containing a mix of clinical and environmental samples) was the same as our W-clade 3. Pop D, equivalent to W-Clade 8, contains the vast majority of clinical samples, as shown by Hatmaker et al. Pop A mainly overlapped with the W-clades 4,5,6,7 of this study17. Notably, 95% of the strains in W-clade 6 were from southern and central China, and only five strains were from the USA, which may result from long-distance migration or geographical spread by atmosphere, with 70% of the strains in W-clade 6 having an aflatoxigenic chemotype. We further construct the neighbor-net network to cross-validate the reliability and accuracy of population structure inference (Supplementary Figure S4a). The network topology is highly similar to the phylogenetic tree, and the W-clades 4 and 5 strains from China also differentiated into the branches. Principal component analysis (Supplementary Figure S4b) and population structure inference (Supplementary Figure S4c) confirm the rationale of dividing the A. flavus population into eight subpopulations. Notably, 94% of the geographical locations of the strains of these clades were from southern and central China in the subtropical or mid-subtropical regions at low latitudes (Figure 1b). We further characterize the genetic diversity and degree of genetic differentiation among different subpopulations. W-clade 6 has high genetic diversity(π = 0.173, p=0.01). In contrast, W-clade 2 (also known as Pop B) has the lowest genetic diversity (Figure 1c). There is a large genetic differentiation (Fst > 0.25) between W-clade 2 and other clades, including W-clade 6 (Figure 1c). The clades show moderate to large genetic divergence from the other clades (Figure 1c). In conclusion, this combined transcontinentally sampled whole-genome phylogeny reveals fine-scale population structure and genetic differentiation across the A.flavus species.

Phylogeographic Diversity Patterns of A. flavus
The genetic diversity of A. flavus shows strong associations with the environmental origins of the strains at a continental scale. Phylogenetic analysis reveals that clinical and environmental samples were clearly clustered, showing genetic differentiation of ecological niches. For example, 51% of clinical samples are distributed in W-clade 8, 20% in W-clade 3, 22% in W-clade 5, and the remaining 7% are scattered in other subpopulations. This pattern is consistent with a recent study on the pathogenicity of A. flavus samples from clinical settings17. On the other hand, it is worth noting that 65% of environmental strains from China in W-clade 8 were phylogenetically interspersed with clinical strains, suggesting that environmental non-aflatoxigenic strains of this subpopulation may have the ability to infect humans19. For example, we observe that the genomes of these clinical A.flavus strains (W-clades 3,8) harboured higher copy numbers of homologue virulence gene 1 (Afu5g03790, Iron transport multicopper oxidase fetC) after large-scale homology clustering comparison of 1003 A. flavus genomes and 149 virulence gene protein sequences from A. fumigatus strains that have been reported to infect humans19 (Figure 1a,d). In addition, the homologue virulence gene 2 (Afu8g00230, Chain A, Verruculogen synthase) is mainly distributed in W-Clade 2 environmental non-aflatoxigenic clade strains. However, overall, the vast majority of the virulence genes were found in both clinical and environmental strains, and the phenomenon of homologues of virulence genes over-represented in clinical strains versus environmental strains was not widespread (Supplementary Figure S5a,b). In environmental samples, sclerotial size S (W-clade 1) and L morphotypes (W-clades 2-8) with different aflatoxin profiles and niche adaptability characteristics were also clearly classified (Figure 1a). The phylogenetic relationships show that phylogeographically distinct clades of A. flavus show clear differences in aflatoxin profiles (Figure 1a). Analyzing from a transcontinental scale, the proportion of non-aflatoxigenic subpopulation strains in W-clades 2 and 8 exceeds 50%. In contrast, W-clades 4, 5, 6, and 7 are dominated by aflatoxigenic strains (Figure 1a,b). The results show a clear latitudinal gradient pattern, that is, the aflatoxin-producing strains were frequently found in low-latitude areas, whereas the non-aflatoxigenic strains are enriched at higher latitudes both within and between clades (Figure 1b). Analyzing strains from China alone, 84% aflatoxigenic strains originate from the south subtropical climate zones (Figure 2a). 62% of the strains isolated from middle subtropical regions had medium or high aflatoxin-producing capacity.

In addition to aflatoxins, pan-metabolome analysis further reveals that phylogeographically distinct clades show extensive metabolic differentiation. Since it is difficult to simultaneously obtain consistent metabolome data from other locations using the same culture conditions and comparable metabolome methods, we used 551 samples from China as a case dataset to discover patterns. A total of 16,432 metabolic features were extracted from 551 metabolome profiles. We find that 36.5% (5,996/16,432) core metabolome features were shared between all isolates. We also observe the presence of many clade-specific metabolites. For example, W-clade 4 and W-clade 5 possess 11.3% (1,861/16,432) and 4.6% (762/16,432) clade-specific metabolite features, respectively (Supplementary Figure S6a). Five largely discrete metabolic patterns could be distinguished following differential abundance analysis between different sets of clades based on annotated metabolites (Figure 2b). We use feature-based molecular networking (FBMN) to organize our LC-HRMS/MS data by grouping mass features with similar mass fragmentation spectra20, accompanied by extensive metabolite annotation efforts using various (GNPS) mass spectral libraries, to boost annotation rates (see Methods). We found that the distribution of a molecular family including cyclopiazonic acid, followed a largely opposite trend compared to that of the aflatoxins molecular family (Figure 2c,d). The highest abundance of the aflatoxin molecular family associated with Chinese clades 2, 4, 5, and 6 (W-clades 4, 5, 6 and 7), while it was found at much lower abundance in Chinese clades 1 and 3 (W-clades 2 and 8) (Figure 2d). This corroborates recent results by Drott et al.14, who showed for the US population that aflatoxin B1 was produced more in Pop A (W-clades 4, 5, 6 and 7) than in Pop B and C (W-clades 2 and 3). Our data therefore demonstrates that most A. flavus isolates have a wide mycotoxin-producing capacity and that low-aflatoxin producing capacity clades often produce other metabolites, including the mycotoxins of Cyclopiazonic acid (CPA).

Evolution of Accessory Genes Drives Metabolic Differentiation
Next, we conduct an exhaustive analysis of the possible drivers that affect the formation of phylogeographic diversity patterns of fungal specialized metabolism. A pangenome representation was built based on 977 representative genomes selected from the phylogeny to characterize commonality and uniqueness attributes. These 977 strain sources cover 9 countries, including 75 clinical strains17 and 902 environmental strains isolated from soil, peanuts, corn, cotton, and other sources14, 16. Orthology inference revealed that the A. flavus pangenome is composed of 15,628 protein-coding genes, which is 3,181 genes more than the reference genome of A.flavus NRRL3357 (GCA_009017415.1)15. These non-redundant orthogroups could be subdivided into core (in all isolates, n = 9,584), accessory (in 5-95% of the isolates, n = 5,085) and ‘unique’ (in <5% of the isolates, n = 959) genomes (Figure 3a) after rigorous removal of bacterial and other prokaryotic sequences by aligning NR database. Among them, the conserved core genome accounts for 61.3%, the accessory genome accounts for 32.5%, and the unique genes account for the remaining 6.2% of genes. Our results and recently reported 59%16 and 42.5%17 proportion of accessory genomes illustrate that the A. flavus population exhibits significant plasticity. However, the ratio of core metabolome (39%) (Supplementary Figure S6a) is lower than the core genome (61.3%). The pangenome exhibits a closed trend as the population size increases (Figure 3b). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis confirm that genes responsible for genetic information processing are dominant in the core genome, alongside genes related to primary metabolism and cellular processes (Figure 3c). 81% of accessory genes and 72% of unique genes were annotated in biosynthesis of primary and specialized metabolites(Figure 3d,e). Clustering analysis of the gene gain/loss of accessory genes (Figure 3f) and pan-transcriptome data indicates corresponding clade-specific expression profiles (Figure 3g), which may be linked to these accessory gene variations directly or indirectly. In conclusion, pangenome and pan-transcriptome evidence indicate that the evolution of accessory genes could drive metabolic divergence in direct or indirect ways.

Variations in BGC Genes Only Partially Explain Metabolic Differentiation
A logical explanation for the differentiation in specialized metabolism could be variation in BGC repertoires between the clades. A total of 52,511 BGCs were identified across all genomes, which were grouped into Gene Cluster Families (1,707 GCFs) using BiG-SCAPE (Supplementary Figure S7a). However, the number of unique BGCs after the first round of automatic deduplication was obviously overestimated, due to e.g. contig breaks inside gene clusters in subsets of the genomes. Therefore, we extracted 2,268 core gene protein sequences of BGCs from these GCFs for detailed curated deduplication by pairwise sequence comparison. Finally, 103 unique BGCs were confirmed in the entire A. flavus population by large-scale sequence alignment and manual validation based on BGC class rules, which amounts to 11 uncharacterized BGCs compared with previous A.flavus pan-BGCs studies14. We matched these 103 BGCs with the 92 BGCs reported in previous literature14 and numbered the BGCs’ ID in the same order. The extrapolation results suggest that there is still a relatively large biosynthetic potential within A. flavus, considering the differences in BGC microevolution sites between strains (Supplementary Figure S7b) and a more stringent sequence similarity threshold to perform similarity clustering and de-duplication on the BGC core gene. The rarefaction curve showed that the entire species could contain ±120 distinct families of BGCs (Supplementary Figure S7c).

Previous research by Drott reported that genetic differences in BGCs core gene can result in differences in SMs production14, highlighting how population-specific SNPs and InDels result in clade-specific differences in BGC core gene content and presence/absence. Lind's analysis of the A. fumigatus population emphasized that the high-frequency variation mainly comes from SNPs and gene gain/loss polymorphisms12. However, partially because most BGCs in A.flavus have not been characterized (16 BGCs products are confirmed by experiments21), previous studies did not statistically or experimentally link core and flanking genes variation to differentiation in the production of known SMs to validate the consequences of genetic variation. To this end, we used a knowledge-based machine learning and mGWAS approach to couple larger paired datasets to comprehensively analyze the variation of BGC core genes and flanking genes and their association with metabolic differentiation patterns. We found that 77% of BGCs core genes showed a conserved pattern, 15% showed presence-absence variation (PAV), and 9% showed a dispersed distribution pattern (Supplementary Figure S8). For example, there is an overall trend of conservation of aflatoxin BGC in A.flavus (Supplementary Figure S8, S9) and Aspergillus species (Supplementary Figure S10) based on large-scale collinearity analysis. This corroborates previous work by Inge et al. in Aspergillus section Flavi, based on comparative genomics11. In addtion, by leveraging the power of paired omics integrated analysis using NPLinker22 and based on structural annotations using reference MS/MS spectra of several known metabolites, we were able to link detected metabolites to candidate corresponding BGCs (Table 1). As a positive control, the aflatoxin B1 and B2 MS/MS spectra are associated with this method to the experimentally characterized aflatoxin BGC (MIBiG:BGC0000006) with high metcalf scores of 4.78 and 6.58, respectively (Table 1). Furthermore, vioxanthin and viopurpurin, two dimeric naphthopyrones produced by Aspergillus species23, were annotated in our metabolomics dataset. These specialized metabolites protect filamentous fungi from a wide range of predators23. The biosynthetic pathway of these two structurally similar compounds is still unknown in A. flavus, and the mass spectra of these molecules were putatively linked to the putative naphthopyrone BGC homologous to MIBiG gene cluster BGC0000107 from Aspergillus nidulans (Table 1). Additional plausible BGC-mass spectral links for seven different A. flavus BGCs with metabolites of three different chemical compound classes pave the way to further dissect their detailed biosynthetic pathways and ecological roles using genetic and biochemical studies (Table 1).

To assess the effect of BGC core and flanking gene variation on differences in metabolite production, we compare in pairs 2 experimentally characterized BGCs with corresponding metabolites detected in the metabolome dataset (Supplementary Figure S11). The results reveal that variations of BGC genes can only partially explain metabolic differentiation. Specifically for the aflatoxin BGC, 34% of W-clade 2 and 14% of W-clade 3,8 strains showed core or multiple peripheral gene loss events (Figure 4a). However, 44% strains with intact clusters still showed no/low aflatoxin production. This is consistent with Drott et al.’s finding that the loss of aflatoxin BGC core gene occurred only in the PopB population14. For Cyclopiazonic acid (CPA) BGC, 29.8% (34/114) of clade 1 and 14.3% (13/91) of clade 3 strains showed whole BGC loss. These strains therefore did not produce CPA. However, 73.3% (11/15), 84.7% (89/105), 70.9% (61/86), and 58.1% (75/129) contains all BGC genes in clades 2, 4, 5, and 6, respectively, while producing little or no CPA. We then analyzed the transcriptional profiles of 58 BGCs (Expressed BGCs genes numbers = 440) of 28 strains in different clades (Supplementary Table 5) and found that 64% (282/440) of BGCs genes were not significantly differentially expressed in different clades, and only 36% (158/440) were significantly differentially expressed (log2|FC| ≥ 1, p.adjust < 0.05) (Supplementary Table 6). The aflatoxigenic Chinese clade 4 (W-clades 5,6) differentially expressed aflatoxin and cyclopiazonic acid pathway genes, while the low-aflatoxins Chinese clade 1 (W-clade 2) highly expressed other BGCs instead, especially kojic acid, aflatrem, and aspirochlorine BGCs (Figure 4b). KEGG enrichment analysis also demonstrated significant differential expression among different subpopulations in aflatoxin biosynthesis pathways (Figure 4c,d,e,f,g). In addition, a linear mixed model (mGWAS) method was used to demonstrate the association and contribution rate of genetic variation and metabolic differentiation in BGC. The mGWAS results were intersected with all predicted BGCs genes in A.flavus, and it was further found that SNP variations in 28% (124/440) of BGCs genes were significantly associated with metabolic differentiation (Figure 5b). Among them, 44 genes showed significant differential expression across clade strains (Figure 4b) . Thus, both analyses can link a subset of variation in BGCs to metabolic differentiation. This suggests that BGC absence/presence variation only partially explains the differential expression at the transcriptional level and the metabolic differentiation outcomes.

Environmental Change Promotes Metabolic Pattern Evolution
The above transcriptome analysis of strains on different clades also revealed that 12 regulators associated with environmental factors such as light (veA and laeA,B), pH (pacC), ion starvation (mscA, mscB, pmr1, hapB, hapE, and hapX), temperature (AFLA_037820), and sensors of carbon/nitrogen (C/N) source (creA and areA) show significant differential expression (Supplementary Table 7 and Figure 5a). We speculate that geographic environment selects certain functional regulators to accelerate the formation of phylogeographic metabolic diversity. The long-term environmental adaptation of fungi does select for specific specialized metabolic repertoires (Figure 2b)14, 24, for example, the latitudinal gradient pattern of aflatoxins (Figure 1a,b)24, 25. However, which environmental factors drive fungal metabolic differentiation remains poorly understood. Jointly with genetic factors, these environmental factors could shape phylogeographical metabolic differentiation. To explore this possibility, the data of 21 climate variables and two geographical factors were collected from the China weather data website (http://data.stats.gov.cn) and soil survey data of 24 soil physiological metrics from our lab. These environmental data correspond to the sampling points one by one. Two widely used genotype-environment association (GEA) pipelines in landscape genomics were employed to identify and interpret the signatures of local adaptation across loci. A total of 2,651 local adaptive genes were identified (G‐Scores > 60, a score threshold) by a logistic regression model using SamBada (Figure 5b). Bayenv2 identified 3,652 SNPs were located in 2,336 protein-coding gene regions (Figure 5b); among these, we identified a total of 803 common adaptively selected genes via two models (Figure 5b). There are 191 overlaps with 2,522 mGWAS genes, 34 intersections with known regulatory genes (n=306), and 33 overlaps with BGCs genes (n=441) (Figure 5b). The important environmental influencers prioritized by a random forest model included soil bulk density, temperature, precipitation, average relative humidity, and soil pH value4 et al. (Figure 5c), all of which showed a notable correlation with latitude across different years (Supplementary Figure S12a~d). We further performed differential expression gene (DEGs) analysis from different origins at transcription level. The genes with the most intersections with DEGs of environmental selection genes were mGWAS genes (n=67), followed by regulatory genes (n=10) (Figure 5d). This demonstrates that environmental changes can promote the evolution of genes associated with different metabolic differentiation patterns (mGWAS) or regulatory genes. It is possible that these environmental factors that vary with the latitudinal gradient jointly shape the phylogeographic metabolic diversity by selecting the above-mentioned environmental selected genetic loci. These environmentally mediated genes have broad functions, mainly related to metabolic regulation (Figure 5e) and physiological responses (Figure 5ef). Among them, 304 environmental selection genes were significantly differentially expressed between strains of different clades. For example, AFLA_057410, a NAD-binding Rossmann fold oxidoreductase, and a NADH-cytochrome b5 reductase of AFLA_029440, and an ATP citrate lyase subunit (Acl) encoded by AFLA_106350 participate in energy metabolism in the tricarboxylic acid cycle, providing energy and precursor materials such as acetyl-CoA for the biosynthesis of specialized metabolites. Fatty acid synthase alpha subunit, encoded by fasA (AFLA_117420) (High Gscore= 90) and fatty acid desaturase of AFLA_089170, known to be linked to the biosynthesis of polyketide specialized metabolites such as aflatoxins, was also detected as a local adaptation locus and in mGWAS. Environmental selection-related signals appear in the tryptophan synthase alpha subunit encoded by AFLA_112810 and a serine/threonine protein phosphatase PPT1 (AFLA_056690); these genes are also significantly differentially expressed between strains from different clades. Previous transcriptome and metabolome enrichment analyses also revealed significant differences in the tryptophan pathway in different clade strains. In addition, a regulatory gene (sfaD, AFLA_093240) for local adaptation (Gscore=92) encodes a guanine nucleotide-binding protein subunit. Environmental selection signals in their gene intervals and significant differential expression were simultaneously detected. This suggests that environmental changes, such as climate change, can accelerate the evolution of genes involved in primary metabolism and regulators. In terms of genes involved in physiological responses, we detected strong environmental interaction signals between bulk soil density environmental factors and two genes encoding a pH signal transduction protein, palA (AFLA_113560, Gscore=112), and a pH-response transcription factor, pacC (AFLA_030580, Gscore=65), using SamBada. In particular, the expression levels of the pacC gene across strains from different clades are significantly different. In addition, four heat shock protein genes, HSPs (AFLA_037820, AFLA_035620, AFLA_084590, AFLA_084460), also showed local adaptation (Gscore=62~63). These genes may be related to environmental temperature adaptability, especially the AFLA_037820 gene (heat shock protein Hsp30-like) (Figure 5a), which has significant differences in expression levels among different geographical origins and clade strains.

Regulatory Variations Drive Metabolic Rearrangement
Deleterious mutations in regulatory genes may trigger larger phenotypic effects through a cascade effect. For example, a partial loss-of-function of the regulatory gene veA has been shown to greatly impact SMs production in fungal co-cultures26. However, some nonsense mutations do not cause metabolic phenotype changes, the mutations in regulatory genes generally cannot be directly associated with metabolic outcomes. We analyzed the differential expression of 306 experimentally validated regulatory genes extracted from peer-reviewed literature. Significant differences were found in the expression of different clades of multiple regulatory genes (Figure 5a). We demonstrated statistical association of mutation sites by intersecting the 2,522 genes previously associated with 60 representative metabolites, for example, aflatoxin B1. We found that mutations in ~28% (86/306) of reported regulatory genes were significantly associated with the five metabolic differentiation patterns (Figure 5d). Among these 86 associated regulatory genes, ~30% (26/86) showed significant differential expression (Figure 5b,d). This illuminated that regulatory variation were significantly associated with metabolic differentiation, significant differences in transcription level was also observed (Supplementary Figure S13). We further examined the deleterious variation (dividing them into low-, moderate-, and high-impact variants) in 36 major regulatory genes that govern specialized metabolism in filamentous fungi. This analysis revealed that clade-specific deleterious variants appeared in regulatory genes, such as the pathway-specific regulator genes aflS, aflR, and the global regulatory gene veA (Supplementary Figure S14). Other variations in specific regulatory genes and the distribution patterns of these deleterious variants in the population are described in the supplementary information. Between the clades, 47% (143/306) known regulatory genes were differentially expressed (log2|FC| ≥ 1, p.adjust < 0.05). For example, the expression of aflR aflS, acon, laeB, msnA, nsdD, veA, and creA regulatory genes was significantly different between the aflatoxigenic (Chinese clade4, PopA3) and non-aflatoxigenic (Chinese clade1, PopB) clades (Supplementary Table 7 and Figure 5a). Excitingly, the functions of a large number of other uncharacterized genes associated with mGWAS have not yet been experimentally verified, providing us valuable avenues for future research (Figure 5b,d). This also reveals that ~86% (C2160/2522) of uncharacterized metabolic regulatory genes or other functional genes associated with metabolic differentiation (Figure 5b) play an important role in driving metabolic differentiation. To verify the association results at the molecular level, we selected four known (sakA, nsdD, apsA, gprJ) and 12 unknown function putative regulatory genes associated with aflatoxin synthesis by mGWAS and conducted gene knockout experiments in LNZW-1 strains (from Chinese clade 5) (Figure 6a,b). Nine of these genes had a significant impact on metabolic differentiation. Compared with the wild-type strain, the AF210 (nsdD), AF420 (apsA), and AF89 (uncharacterized gene) gene knockout strains (Figure 6c) have smaller growth diameters, especially the AF210 and AF420 gene knockout strains (Figure 6c). Gel electrophoresis images of wild type (top) and transformant (bottom) of the gene of AF210 (nsdD) (Figure 6d), AF420 (apsA) (Figure 6e), and AF89 (uncharacterized gene) (Figure 6f), demonstrating successful gene knockout. The untargeted metabolomic analysis of these mutant strains validated that a single gene mutation can underpin metabolic switches (Figure 6g). The concentrations of aflatoxin B1, G1, cyclopiazonic acid and kojic acid were significantly downregulated in 9 out of the 16 knockout strains compared to the wild-type strain (Figure 6h-k). The differential and pathway analysis based on the metabolomes demonstrated that single gene mutations resulted in significant metabolic perturbations (Supplementary Figure S15a,c,e), mainly enriched in primary metabolism pathways, such as phenylalanine and tryptophan biosynthesis (Supplementary Figure S15b,d,f). These regulatory mutations changed the chemical profile of deletion strains (Chinese clade 4, PopA) to look more like the chemical profile of strains in Chinese clade 1, PopB.

mGWAS results had linked primary metabolic genes in glycolysis pathways (AFLA_133830 and AFLA_053990) to aflatoxin B1 production (Figure 6a). The differentially expressed genes of Chinese clades 1 and 4 were significantly enriched in the glycolysis pathway, TCA cycle, and aflatoxin pathway (Figure 4c), thus supporting this mGWAS association result and the general hypothesis that variations in primary metabolism genes are significantly associated with specialized metabolism differentiation. To further address this hypothesis, we examined evidence from three levels: (1) Genome-wide selective sweep analysis (XP-CLR method27), a well-established method in population genetics natural selection pressure analysis, was used to screen for genes showing signs of recent natural selection in each clade. KEGG enrichment analysis revealed differences in the type or number of metabolic genes selected across subpopulations. 58% to 66% genes with signatures of selection were part of pathways belonging to primary metabolism (Figure 6l,m and Supplementary Figure S16 a1-14). Among them, amino acid metabolism (such as tryptophan metabolism), lipid metabolism, and carbohydrate metabolism are over-represented (Supplementary Figure S16b,c1-14). In addition, 17% to 29% genes with signatures of selection were related to specialized metabolism, and again signatures were often specific to one or a few clades, for example, aflatoxin, aflatrem, and aspirochlorine BGCs. (2) We performed KEGG enrichment analysis on the genes associated with mGWAS of different 60 metabolites in the previous five metabolic patterns and differentially expressed genes from different clades, looking for those primary metabolic genes that are significantly associated with metabolic differentiation. The results show that the top five metabolic pathways with the most genes associated with mGWAS are carbohydrate metabolism, amino acid metabolism, lipid metabolism. (3) Lastly, we looked at the correspondence (co-clustering) of transcriptome profiles with metabolome profiles of precursor biosynthesis pathways (e.g., amino acids, etc) in primary metabolism across the 28 strains. Differentially expressed genes from different clade strains were also significantly enriched in glycolysis, citrate cycle (TCA cycle), pentose phosphate pathway, tryptophan biosynthesis, fatty acid biosynthesis, glycerolipid metabolism, sphingolipid metabolism, glycosphingolipid biosynthesis, pyruvate metabolism, and aflatoxin biosynthesis(Figure 4c,d,e,f,g). Correspondingly, the metabolomic data also showed that primary metabolite abundances show consistent clade-specific differences consistent with those of SMs. For example, we found that tryptophan is present at higher levels in Chinese clades 1 and 3 than in clades 2, 4, 5, and 6 (Figure 7a). The distribution trend of the abundance of cyclopiazonic acid, which incorporates tryptophan as a precursor, is consistent with that of tryptophan (Figure 7b). Similar patterns exist between metabolic abundances between other primary and specialized metabolites in the corresponding pathway, such as oleic acid ethyl ester, averantin (AVN), versicolorin A and aflatoxin B1. It is also worth noting that many important primary metabolites, including nicotinamide adenine dinucleotide (NAD) and UDP, have higher abundance in Chinese clades 1 and 3 (Figure 7c,d). The mirror plot demonstrates the high reliability of these four metabolites' annotations (Figure 7e,f,g,h). Varying abundance of cofactors and associated energy metabolism thus seems to be associated with differences in activity of specialized metabolism28. Based on the above data, we summarised and described a model of the main genetic and environmental drivers of fungal secondary metabolic diversity (Figure 7i). The evidence emphasizes that the transcription levels of primary metabolic pathways differ between strains. These differences likely ultimately mediate the differentiation of specialized metabolism. Altogether, the above evidence indicates that variation in primary metabolism flux and regulatory gene sequence significant impact on diversity in specialized metabolite production.

Methods
A. flavus strains library construction and mycelium cultivation
We collected soil, peanut, or corn samples across China from 2014 to 2019, and isolated 3,567 Aspergillus flavus strains by single spore isolation, morphology combined with ITS sequencing identification. The latitude of sampling points in the China region spans from 24.3 north (Zhanjiang, Guangdong province) to 46.4 (Qiqihar, Heilongjiang province). Longitude from 87.3 E (Changji, Xinjiang Uygur Autonomous Region) to 123.42 (Qiqihar, Heilongjiang Province) in the west. Altitudes range from 5.3 m (Lianyungang, Jiangsu Province) near sea level to 2330m (Zayü County, Tibet Autonomous Region)24. All strains were isolated from independent samples. Detailed strain sampling site coordinates were described in Supplementary Table 1. We selected ~600 representative A. flavus strains from the library as research objects by integrating information such as geographic origin, phylogenetic relationships constructed using ITS sequences, climate zone, and APCs.

A. flavus was incubated on a Potato Dextrose Agar (PDA) medium and then put on Petri dish plates at a constant temperature and humidity incubator in dark conditions at 29±1 °C for 10 days. We used a 0.1% Tween-80 solution to wash the conidia and obtain a suspended spore mixture. Conidia suspension mixture was used, which was quantified using a hemocytometer and a light microscope38. 2.5×105 conidia/mL mixture was inoculated into 50 mL of autoclaved liquid medium, containing 0.25% yeast extract, 0.1% K2HPO4, 0.05% MgSO4 · 7H2O, and 10% glucose (pH= 6.0). The liquid culture flasks containing conidia of A. flavus were placed in a shaker (180 rpm) at 29±1 °C in dark conditions. After 5 days, the mycelia were acquired by quickly filtering the cultures from the flasks with multiple layers of cheesecloth and washed with 10 mL of 4 °C saline solution (0.9% NaCl), and the mycelia were squeezed with sterile absorbent paper. Then, A.flavus mycelium was quickly transferred into a 50 mL centrifuge tube and quenched with liquid nitrogen. We used a split strategy to sequence one of the collected mycelium samples for acquiring genome data and the other for metabolome research. The mycelial samples were directly weighed for the DNA extraction experiment. However, for metabolomics experiments, the mycelium samples were freeze-dried and then weighed.

DNA extraction, library preparation, and genome sequencing
We used the following DNA extraction protocol: 200 mg of fungal tissue was weighed in an EP tube with a one-thousandth balance, quenched with liquid nitrogen, and placed in a freezer grinder to homogenize it into powder. 1.2 mL CTABC lysis buffer was added, and the mixture was vortexed for 60s. The mixture was incubated at 65°C for 60 minutes. The mixture was cooled and centrifuged at 13,201 × g for 6 minutes. Then, the supernatant was transferred to a 2.0 mL EP tube, an equal volume of phenol/chloroform/isoamyl alcohol (25:24:1) reagent was added, and the mixture was shaken. The mixtures were centrifuged at 13,201 × g for 15 min, and their supernatant was transferred into a 1.5 mL EP tube. Then, 2/3 volume of -20°C pre-chilled isopropanol was added into the tube, shaken up and down to mix evenly, and then put in a -20°C refrigerator for more than 2 hours. The mixtures were centrifuged at 13,201 × g for 15 min, and then their supernatant was aspirated into a new 1.5 mL EP tube, supplemented with 750 μL of 75% ethanol, and the pellet was rinsed with a pipette. Finally, the tube was centrifuged at 13,201 × g for 5 min at ambient temperature, and the supernatant was aspirated. After centrifugation, the residual liquid was removed, air-dried for 3-5 min, and solubilized in appropriate amounts of TE solution (EB solution for Pacbio sequencing samples).

For next-generation sequencing (NGS), 1 μg of genomic DNA was taken and disrupted by sonication using a Covaris instrument. Fragment selection was performed on the fragmented samples so that the sample bands were concentrated around 200-400bp. Then, a reaction system was prepared, the end of the double-stranded cDNA was repaired, and an A base was added to the 3' end. Furthermore, a linker ligation reaction system was constructed to connect the linker to the DNA. A PCR reaction system was then prepared and a reaction program was set to amplify the ligated products. Amplification products were purified and recovered by magnetic beads. After the PCR product was denatured into single strands, a circularization reaction system was prepared, and the reaction system was fully mixed and incubated to obtain a single-stranded circular product. After digesting the linear DNA molecules that had not been circularized, the final library was obtained. Qualified NGS libraries were performed on the Illumina (HiSeq X-Ten) platform at BGI (Shenzhen, China) to generate the 150 bp paired-end raw reads.

For whole-genome long-read sequencing, 1 μg of genomic DNA was fragmented using a Covaris instrument. The fragments longer than 20 kb were selected for the PacBio library construction according to the manufacturer’s instructions, and qualified libraries were sequenced on the PacBio Sequel II system at BGI (Shenzhen, China).

Quality filtering and genome-wide variation detection
We used SOAPnuke (v1.5.6)39 to filter out adapter contamination reads, PCR duplication reads, and low-quality and ambiguous reads with the parameters “-q 0.2 -l 0.2 -n 0.05 -d” for the NGS data. High-quality sequencing reads from each individual were aligned to the reference genome of A. flavus NRRL 3357 (Accession No. GCA_000006275.2) using Burrows-Wheeler Aligner with the BWA MEM algorithm. Subsequently, duplicate reads and realigned indels were processed using Sentieon (https://www.sentieon.com/), and then the gvcf file of each individual by using the Haplotyper method of Sentieon. GVCFtyper was used to integrate the variations detected from all individuals. Thereafter, SNPs were filtered by GATK(v4.0.2.1), which match the condition “QD < 2.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || FS>60.0 || SOR>3.0”. High-quality InDels were filtered under the parameters of “QD < 2.0 || ReadPosRankSum<-20.0 || FS>200.0 || SOR>10.0”. Single-nucleotide polymorphism (SNP) variants used for population analysis were further processed by vcftools with parameters “--maf 0.05 --minDP 4 --minGQ 10 --minQ 30 --max-missing 0.99”. The impact of both SNPs and InDels was annotated by SnpEff(v4.3) with default parameters.

De novo assembly, decontamination, and assessment of genome assemblies
For the Illumina sequenced data, to obtain the ideal assembly results, we use SPAdes (v.3.15.1)40 to assemble the high-quality clean reads into genomic sequences with k-mer ranges from 33 to 83 by a step size of 10. Genome assembly for Pacbio RSII sequencing data was performed using wtdbg241 with the parameters “-L 5000 -x preset2”. Clean reads from Illumina sequencing data were used to correct INDEL and SNP errors for the assembly sequences, and Pilon (v 1.24) was used to perform three rounds of polish error correction for possible assembly INDEL and SNP errors with parameters “--diploid --changes –verbose”.

To remove the bacterial contamination sequences, the Non-Redundant Protein Sequence (NR) and Nucleotide Sequence (NT) database (v20180315) was used to align against the assembled genomes using BLAST. The alignment results with E-value < 1e-5, align length > 100bp, and identity > 70 for NT, >40 for NR were deemed significant. Sequences not matching to the Aspergillus species were considered as contaminations and then removed to obtain the final assembled genome. The completeness and redundancy of genomes were assessed based on single-copy ortholog analysis with BUSCO (3.0.2) using the Eurotiomycetes_odb database. The Virulence genes were collected from the J. Rhodes et al. research42, and OrthoFinder(v2.4.0)43 with default parameters was used to find homologous genes in China and America A. flavus samples, respectively.

Genome annotation, functional classification, and enrichment analysis
There are two main types of repeats in the genome: one is tandem repeats (Tandem Repeat, TR); the other is interspersed repeats (Transposable element, TE). For repeat region annotation, we used Tandem Repeats Finder (TRF, v4.0C) to annotate tandem repeats. TEs' homolog annotation was carried out by compared with a known repeat sequence library (Repbase, v16.02) by using Repeatmasker and Proteinmask.

Gene structure annotation was performed by two methods: homology annotation (Homolog) and model-based de novo prediction (De novo), respectively, and EVidenceModeler integration software was used to combine different annotation results: (1) homology-based annotation: Five closely related species sequences (GCA_000002655.1, GCA_000184455.3, GCA_009193485.1, GCA_009176385.1, GCA_000006275.2) were selected and downloaded from NCBI, and the genome was annotated using GeneWise; and (2) de novo annotation: gene prediction of assembled genomes using Augustus and GeneMark. The weight of the homologous annotation result was set to 10, and the weight of the two de novo software prediction results was set to 1.

The A. flavus genes were functionally annotated as detailed as possible using the six mainstream functional databases, including NR, Interpro, GO, KOG, KEGG, and SwissProt. We used Blastn (v2.2.31) for NT database (v20180315) annotation with E-value < 1e-5 and identity > 70%, Blastp (v2.2.31) for NR (v20180315), KOG (v20090331), KEGG (v84) and SwissProt (release-2017-09) annotation with E-value < 1e-5 and identity > 40%, and Blast2GO (v2.5.0) combined with nr-based annotation results for GO (v20171220) annotation. The phyper function in R was used for enrichment analysis of GO and KEGG annotations, based on Fisher’s exact test.

Population structure and phylogenetic analysis
By using admixture (v1.3.0)44, we assessed population structure from k=2 to k=10 with default parameters. The maximum likelihood estimation was applied without using the prior population information, and the most likely number of ancestors was determined using cross-validation (CV) error. Principal component analysis (PCA) was performed with Plink v1.945, and the first two eigenvectors were plotted. The θπ ratio and population genetic differentiation index (Fst) between each clade were calculated using VCFtools v0.1.13, respectively, with a 5 Kb non-overlapping slide window. By screening the genome in windows, we detected candidate regions with selection signatures with the largest differences in θπ ratio (θπ-clade-i /θπ-clade-j，bottom and top 5%), and the top 5% FST(population differentiation coefficient values) regions between two clades. The intersections of the θπ ratio and FST candidate regions were defined as the regions that were under selective sweep. Then, the genes within or overlapping the regions showing a selection sweep were defined as candidate genes. Linkage disequilibrium (LD) decay of each clade was calculated by using PopLDdecay v3.40 with default parameters.

The genome-wide SNPs were used to construct the phylogenetic tree for 1,064 samples using the Neighbor-Joining method using TreeBeST (v1.9.2) with parameters “-b 1000”. For other group phylogenetic analysis, we employed the maximum likelihood method in IQ-TREE2 (v2.0.4)46 with parameters “-alrt 1000 -bb 1000”, then, the tree output was displayed by the iTOL (http://itol.embl.de). A.flavus has three mating types (MAT1-1, MAT1-2, and MAT1-1-MAT1-2 ), which are encoded MAT alpha 1 (MAT1-1) gene (Genebank Accession Nos: EU357934) and the HMG-box protein MAT1-2 (MAT1-2) gene (Genebank Accession Nos: EU357936)47. We downloaded these two genes from GeneBank and then used Blast to scan the genome for mating typing.

Pangenome construction and analysis
An A.flavus pan-genome included 25 PacBio assembly genomes and 977 NGS assembly genomes was built by EUPAN (v0.44)48. Genomic quality assessment was first performed using QUAST using the longest contiguity genome HNSY-2B (Has the highest N50 value) as a reference. The unaligned contigs from individuals were merged, and then the redundant contigs were removed to generate the non-redundant contigs using CD-HIT (v4.6.3) with sequence identity < 0.9. Subsequently, a BLAST-based method was used to cluster the reference sequences and non-redundant unaligned contigs using the “eupan rmRedundant blastCluster” command with parameter “-c 0.5” to generate the comprehensive A.flavus pan-genome dataset, and the pan-genome sequences were obtained from the previous sequences annotation. We mapped raw reads to the comprehensive sequences to determine gene presence-absence. A “map-to-pan” strategy subsequently was utilized to determine gene presence-absence and gene family presence-absence by using Bowtie2 (v2.2.5). The gene body coverage and CDS coverage of each gene were calculated after mapping the reads to the pan-genome sequence. Finally, we determined gene presence-absence patterns considering gene body coverage value (>0.8) and CDS coverage values (> 0.95). If one member of a gene family satisfying the above criteria was present in a given A.flavus genome, the gene family was considered as present. The output of the presence/absence variation (PAV) profile file was further used to plot the rarefaction curves and perform statistical analysis.

Biosynthetic gene cluster deduplication and validation
We used antiSMASH49 to predict BGCs across all 651 de novo assembly and annotated A. flavus genomes from China and the USA. BGCs in all genome sequences were clustered to deduplicate similar BGCs sequences into Gene Cluster Families (GCFs) and compared to experimentally characterized reference BGCs from the MIBiG (v2.0) repository using BiG-SCAPE (v1.0.0)50. The biosynthetic gene clusters (BGCs) in all genomes were identified with antiSMASH (v6), using the genome assembly results and gene annotation gff3 files as input files. These samples have corresponding metabolomic data, which facilitates subsequent paired comparisons between BGCs and secondary metabolites. BiG-SCAPE analysis was performed using default settings, separating the analysis according to the BGC product type and creating network directories for each class. In addition, mixing all classes and retaining singletons was also performed to deduplicate similar BGCs sequences into Gene Cluster Families (GCFs) using a threshold of 0.3 to 0.6, with an interval of 0.05, and then identify which BGCs are shared between which strains. This resulted in a total of 52,511 BGCs that were grouped into 1707 GCFs. The core gene protein sequences of these BGCs were further extracted from these 1707 GCFs, and then 2268 core genes were obtained (some BGCs have more than 2 core genes). Pairwise comparisons were performed using diamond (v0.8.23.85) with the following parameters: --evalue 1e-5 --sensitive. Then, according to the comparison results, the core genes were redundancy-filtered according to different similarities (amino acid identity value = 40, 50, 60, 70, 80, 90) . After taking their intersection, we further manually checked and removed the duplicates. The network files were visualized using Cytoscape (v3.8.2). The BGC types encoding production of miyakamides and desferri-ferrichrysin were manually annotated, since MIBiG (v2.0) had not yet included these two BGCs. Based on the results, a gene cluster family (GCF) presence/absence matrix was constructed. This was used to screen for conserved BGCs and clade-specific BGCs. To identify clade-specific GCFs, we performed a dimensionality reduction analysis based on the presence/absence matrices of GCFs to obtain the linear combination matrices and compared the combination matrices by principal component analysis (PCA) function in R.4.1.0. Rarefaction curves analysis was constructed by “iNEXT” package of R to extrapolate biosynthetic diversity potential of the whole A. flavus populations. The t-Distributed Stochastic Neighbor Embedding (t-SNE) analyses were conducted with the “Rtsne” package of R.

Biosynthetic gene cluster variation and comparative analysis
Two recently developed pipelines were used for finding clusters of co-located homologous sequences (cblaster v.1.2.951) and visualizing the BGCs comparison (clinker v0.0.2052) for a large-scale study population BGC variation. The specific analysis steps are as follows: 1. Building a cblaster database using the assembled genome and annotation files of each sample with the command “cblaster makedb input.gff myDB”. A total of 570 sample databases were constructed, including an old reference genome version of A.flavus NRLL 3357 and three new versions of A.flavus NRLL 3357 genomes, including GCA_009017415.1, GCA_014117465.1, and GCA_014117485.1. 2. Detecting the presence/absence of the gene cluster in each sample in turn according to the protein sequence of each BGC, to facilitate comparison with BGCs that have been experimentally characterized. BGC protein sequences were predicted by SMURF to determine the BGC start nucleotide of the first biosynthetic gene (5′ end) and the stop nucleotide of the last biosynthetic gene (3′ end) with the AFLA_x gene ID, and these coordinates were manually curated. Then, these BGCs were fed into cblaster as input files to locate and align the BGC in a specific genome with the command “cblaster search -m local -db myDB.dmnd -jdb myDB.json -qf input_pep.fa -p output.html -o output.alignment”. 3. Creating a json file for the next step to extract the best alignment sequence gbk file. 4. Using the cblaster extract_clusters module to extract the gbk files aligned for each BGC in each sample for clinker input. 5. Leveraging cblaster to extract the gbk files of each sample (each BGC) as clinker input; call clinker to get the visualization results of each group of gene clusters. In order to facilitate the subsequent comparative synteny analysis, we selected every 20 genomes as a group for clinker analysis. The BGC core (backbone) gene was identified by antiSMASH 549. The protein sequences of these core genes were then extracted and fed into cblaster to scan all isolate genomes for analyzing the identity values of the BGC core gene. A threshold of 80% was then employed to determine the presence/absence of BGC core genes.

Sample preparation and quality control (QC) for metabolomics analysis
The 566 fungal samples (biological replicate n=1) were separated to two aliquots, with one for genome sequencing and one for metabolomics data acquisition. In our subsequent subpopulation comparative analysis, we treated different strains within the same subpopulation as biological replicates. The liquid culture flasks containing conidia of A. flavus were placed in a shaker (180 rpm) at 29±1 °C in dark conditions. After 5 days, the mycelia were acquired by quickly filtering the cultures from the flasks with multiple layers of cheesecloth and washed with 10 mL of 4 °C saline solution (0.9% NaCl), and the mycelia were squeezed with sterile absorbent paper. Then, A.flavus mycelium was quickly transferred into a 50 mL centrifuge tube and quenched with liquid nitrogen. The mycelium samples were freeze-dried and then weighed. To obtain metabolite extracts, an optimized two-step extraction method24 is compatible with polar and non-polar metabolites from mycelial samples that maximizes the scope of extraction. The specific steps are as follows: 1) the mycelium samples were freeze-dried, and then 50 mg per lyophilized mycelial sample was extracted by vortexing vigorously for 1 min with 1 mL of extraction solutions 1 (methanol: acetonitrile: water 2:2:1 v/v/v), which contained the internal standard of 2-chlorophenylalanine (20.0 µg/mL) and camphanic acid (25.0 µg/mL). 2) We added 6~8 beads in each tube and homogenized in a ball mill for 4 min at 45Hz, then it was treated for 5 min with ultrasound of 700W (incubated in ice water). Hereafter, all tubes were sonicated for 5 min with 700W (incubated in ice water). This process was repeated one more time using solution 2 (methanol: dichloromethane: ethyl acetate 1:1:1 v/v/v) to further extract the nonpolar metabolites. 3) The proteins were precipitated by incubating the homogenate for 1h at -20 ℃ followed by centrifugation at 20000 g/min for 10 min at − 20 °C. 1 mL of the resulting supernatant was transferred into a fresh 2 mL LC-MS/MS glass vial for the UPLC-HRMS (Ultra performance liquid chromatography-tandem high-resolution mass spectrometry) analysis. Based on the minimum quality control requirements of the metabolomics community, three quality control methods were used for data quality control throughout the metabolomics analysis process in this study. Quality control measures such as the introduction of quality control samples, the addition of internal standards, and blank samples ensured that we could collect high quality metabolomics raw dataset during this period53. Firstly, two internal standards were added to the extraction solvent to monitor experimental deviations in the sample pre-treatment process to ensure the quality of the data. Meanwhile, a widely used assessment criterion was employed to assess the quality of metabolomic data. Namely, If the data from quality control (QC) samples were all closely clustered to the origin in the PCA scores plot, indicating that the quality of data was suitable for subsequent analysis. Furthermore, we inserted a sample of blank solvent per 10 samples to monitor the cross-contamination between samples. Quality control (QC) samples were prepared by pipetting 10 µL from each sample extract into a QC injection vial using a 100 µL pipette. We set up a tube containing two extraction solvents as an extraction blank. This tube did not contain any A. flavus mycelium sample, but instead contained 4~6 pieces of bread and two internal standard solvents. Solvent blank samples were prepared by pipetting 1 mL of pure methanol solvent (99.99%) into the injection vial named “blank”. The 551 samples were divided into 6 batches for analysis using UPLC-HRMS, because the sample tray of the UPLC-HRMS instrument can only accommodate 120 sample vials, and the mobile phase is only sufficient for one batch. A diagram of the sample run sequence shown in our previous research articles of supplementary figure S124. Each batch contains 120 samples, one solvent blank sample, and one extraction blank sample. However, in the LC-MS run sequence, the solvent blank sample and quality control sample will be injected repeatedly every 10 samples to assess cross-contamination between samples and instrument stability. Before collecting mass spectrometry data from A.flavus mycelial samples, we first collected 10 quality control samples as technical replicates to allow the UPLC-HRMS instrument to reach mobile phase equilibrium and a stable state. The mycelial samples and extraction blank samples were collected only once. For the actual mycelial samples, we did not set up technical replicates; we used QC samples to assess the reproducibility of the samples directly. Therefore, each batch generated 12 solvent blank samples and 22 quality control samples, resulting in a total of 72 solvent blank samples and 132 quality control samples across six batches, for which metabolomic data files were obtained. Ultimately, we collected metabolomic data from 673 A.flavus strains. However, due to failures in the genome sequencing process for some strains, we had to discard the metabolomic data for those strains. There were 551 strains that possessed both genomic and metabolomic profiles. The quality control analysis results shown in our previous research articles of supplementary figure S124. The high-quality data were proved by QC replicates, which are tightly clustered in PCA space. We also used the peak alignment results of the metabolite Versicolorin B to assess its distribution in different A.flavus strain samples, quality control samples, and blank samples. We found that cross-contamination was minimal and reproducibility was excellent, indicating that we obtained high-quality metabolomic data24.

Metabolomics data acquisition by UPLC-HRMS
All the isolates’ non-targeted metabolome raw data were collected via our optimized and standardized metabolome platform24, 38. A C18 column (Hypersil Gold, 100 mm × 2.1 mm (i.d.)) with 3 µm pore size (Thermo Fisher Scientific, USA) coupled with the Ultimate 3000 system (Dionex, Sunnyvale, CA, USA) was employed for liquid chromatography separation. Mobile phase A is water: methanol (95:5, v/v) containing 0.1% formic acid, and 10 mM ammonium formate. Mobile phase B is composed of methanol: water (95:5, v/v) containing 0.1% formic acid and 10 mM ammonium formate. Both mobile phases were mixed following the program below in the UPLC system at a flow rate of 300 μL/min as a gradient mobile phase for effective separation of different components. The gradient elution program was set to: 0~1 min: 15% B (mobile phase B), 1~3 min: 15–50% B, 3~5 min: 50-70% B, 5~10 min: 70-100% B, 10~13 min: 100% B, 15min: 100~15%B, and from 15~20 min: 15% B. Subsequently, 2 μL of the extracted metabolites was injected into Orbitrap Fusion mass spectrometer (Thermo, USA) for mass spectrometry raw data acquisition.

In this study, full scan mass spectrometry (MS1) data were used for relative quantitative analysis, and mass fragmental (MS2) data for metabolites identification and annotation. We acquired the mass spectrometry data containing MS1 and MS2 mass spectra by a full-scan MS coupled with data-dependent fragmentation (DDA) mode under Ultra Performance Liquid Chromatography Tandem Orbitrap Fusion High-resolution Mass Spectrometry (UPLC-HRMS). For parameter optimization, we refer to the optimal method parameter optimization results of François’ lab using LC-Orbitrap fusion mass spectrometry for metabolite annotation54. The positive and negative mode raw data were separately acquired to enhance data quality. The HESI (Heated Electron Spray Ionization) ion source was used to efficiently ionize the UPLC effluent components. The -1.9 kV capillary voltage and 320 °C ion transfer tube temperature parameter settings were selected. The specific parameters of the ion source are as follows: 320 °C heater temperature for mobile phase volatilization and ionization; sheath gas (Arb): 40; aux gas (Arb): 5; and sweep gas (Arb): 0; spray voltage: static, positive ion: 3.5kV, negative ion: 3.0kV. The full-scan MS1 stage was identical to that for the positive and negative models except for polarity. The specific parameters were as follows: Orbitrap detector resolution: 12000, scan range (m/z): 100~1200, RF lens (%): 60, AGC target: 5.0e5, maximum injection time (ms): 100, and source fragmentation: disabled. The acquisition of mass fragmentation (MS2) spectra requires a smaller number of ions; otherwise, it will affect the number of primary mass spectra (MS1) and finally affect the quantification. Filter parameters included: intensity threshold at 1.0e4, charge state: 1~2, and dynamic exclusion after 1 time. The exclusion duration was set to 60 seconds. Top speed mode was selected in the data-dependent acquisition mode, and the number of scan event types: 1s. The DDA section specific parameter settings were identical as for positive and negative mode, except for the high-energy collision-induced dissociation (HCD) setting. A type of stepped collision energy in HCD activation in positive mode was set at 40  5 eV and 33% for fragmentation of the isolated precursor ions, and the HCD cell pressure was set to 8 mTor. In negative ionization mode, the value of HCD was set at 30  5 eV in negative mode. Other parameters for both ionization modes included: Orbitrap detector mass resolution of 30000 FWHM, AGC target: 5.0e4, and the maximum injection time was set to 100. The quadrupole with a narrow isolation window of ±1 Th to isolate the MS2 ions. Submission of data acquisition sequences, instrument control, and partial data processing task via Xcalibur 4.0 software (Thermo Fisher Scientific).

Metabolome features extraction and data quality control
The metabolome feature extraction, feature alignment, and gap filling steps were performed in MZmine 4.2755 (Linux version). First, LC-HRMS positive and negative raw mass spectrometry data files were converted to mzXML format by the ProteoWizard using default parameters. In MZmine 4, the batch processing mode (mzwizard tab) was selected to achieve streamlined mass spectrometry data processing. The MS1 and MS2 noise levels were set at 1E5 and 6E3, respectively. The final MS2 noise threshold (6E3) selection was determined by comparing the number of core and unique metabolites identified on the GNPS platform by Feature-Based Molecular Networking (FBMN) analysis across eight preset MS2 thresholds (1E0; 1E2; 1E3; 2E3; 4E3; 6E3; 8E3; 1E4) in MS2 peak detection step (Supplementary Figure S25). A minimum group size of scans set 5, a minimum highest intensity set 2E5, and an m/z tolerance set 0.01 Da in the second step of the ADAP chromatogram builder. We selected the baseline cut-off deconvolution algorithm to perform the peak deconvolution with the following settings: min peak height: 1.0E5; peak duration range(min): 0.05-2.0; Baseline level: 1E5; m/z range for MS2 scan pairing (Da): 0.01; RT range for MS2 scan pairing(min): 0.3. We used the isotopic peaks grouper algorithm with an m/z tolerance of 0.01 Da (or 10 ppm) and an RT tolerance of 0.2 min for grouping isotopes. The join aligner module (m/z tolerance = 0.01 Da, weight for m/z = 2, weight for RT = 2, absolute retention time tolerance = 0.2 min) was used to do feature alignment. Feature filtering was performed by feature list rows filter selection, with a minimum peaks in a row = 2, and a minimum peaks in an isotope pattern = 2. The peak gap was filled with the RT = 0.2 min and m/z tolerance of 0.01 Da. Eventually, a .csv quantitative file was exported and then was submitted to GNPS20 and MetDNA56 after modifying the data matrix format for metabolites identification or molecular networking analysis and structural annotation. We followed the metabolomics data quality control methods and data quality assessments as previously reported24. The raw metabolomic data of 95 strains from the United States generated by Drott et al.14 were downloaded from GNPS MassIVE repository (ID no. MSV000087134, https://doi.org/doi:10.25345/C54226). They used a target method of PRM mode to acquire the data on the Thermo Scientific Q Exactive Orbitrap mass spectrometer. For the 94 USA strain samples, they only acquired MS1 information for quantification and then collected 24 qualitative files to acquire the MS/MS spectral information containing target compounds (We finally contacted Drott and Tomás to upload the 24 qualitative files containing the MS/MS mass spectra to the GNPS repository). Drot's PNAS article only studied 13 experimentally verified pathways for secondary metabolites, but did not study those pathways that have been isolated and identified in A. flavus but whose corresponding unknown BGCs have not yet been resolved. The identification of 13 metabolites was completed by Tomás and Drott, and they provided MS/MS spectra information of 13 secondary metabolites in the PNAS article14 (Table S4 and SI Appendix, Figs. S5–S29 in this article). Four of the secondary metabolites belong to annotation level 2, and the other nine metabolites belong to annotation level 3. We merged the metabolite quantification table sent to me by Drott with our data matrix based on the MS1 accurate mass (±5 ppm, m/z tolerance window) and MS/MS spectral information from the PNAS literature (annotation level 3). However, this data fusion method does have risks. Therefore, we performed MS1 and MS2 detection using the same Mzmine 4 workflow and FBMN analysis on the updated files containing MS/MS data with the same parameters, aiming to merge the two datasets based on MS1 accurate mass (±5 ppm, m/z tolerance window) and MS/MS spectrum comparison results (annotation level 2). We only retained the MS1 metabolic features that had MS/MS spectra. From the Mzmine analysis results, we can see that Drott et.al’s dataset contains only hundreds of metabolite MS/MS spectra information, which is consistent with their goal of targeting the analysis of a few secondary metabolites (Figure S17). However, compared with the number of MS1 reported in their article14 (Fig 2A in the PNAS article), the number of metabolic features (MS1) containing MS/MS spectra is small. Therefore, it is difficult to perform FBMN joint analysis with our dataset, rich in MS/MS spectra information. We therefore performed a separate FBMN analysis and only matched 3 metabolites with the GNPS reference library, including aflatoxin B1, cyclopiazonic acid, and aflatrem (GNPS job results link: https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=3532f01f1f8a47c9b2e36eca093b89bc). This is consistent with the results of annotation level 2 shown in Table S4 in their PNAS14 article (annotation level 2), except for the metabolite Ditryptophenaline. We couldn't find it in the GNPS library. In order to comply with the minimum standards proposed, we therefore selected two A. flavus secondary metabolites, including Aflatoxin B1 and Cyclopiazonic acid (±5ppm, m/z tolerance window, and with MS/MS spectrum alignment mirror plot ) for data fusion to verify that the presence/absence variation of the BGC gene in the strains from China and the United States in our manuscript can only partially explain the changes in secondary metabolites.

Metabolite annotation, identification, and molecular networking analysis
Metabolite identification and annotation were performed using a four-step approach. First, we used our in-house mycotoxin library, which contains 262 commercial or isolated chemical standards, and a manually curated compound list based on accurate mass (m/z, ± 5 ppm), MS/MS spectra, retention time, and isotope patterns. If we found a match, we considered this MSI Metabolite Identification Level 1. Second, metabolites were further identified based on accurate mass, isotope pattern, and MS/MS spectra against public databases by spectrum alignment of tandem mass spectrometry (MS/MS), including mzCloud, GNPS20 (all public mass spectral libraries included there), MoNA, MassBank, and NIST. Specifically, tandem mass spectrometry (MS/MS) data of samples were acquired using an Orbitrap Fusion mass spectrometer (Thermo Fisher, USA). To annotate more likely metabolites, we used three widely available commercial or open-source metabolome software, including Compound Discovery 2.0, GNPS, MetDNA3. For Compound Discovery 2.0, the detailed parameter settings of Compound Discovery 2 (Thermo Fisher Scientific, USA) were the same as in our previously reported method24. The detailed steps are as follows: 1) metabolomics raw data were imported into the software of Compound Discovery (CD) 2.1(Thermo Fisher Scientific, USA) to generate a data matrix that consists of the retention time (RT), mass-to-charge ratio (m/z) values, peak intensity, and annotation of metabolites. 2) Metabolite annotation consists of the following steps: ① Unknown peaks were aligned and detected with the parameters of RT tolerance of 0.8 min and 5 ppm mass deviation. Minimum peak intensity set to 10000 and S/N threshold=3. ② The compounds were annotated with different types of databases. Among them, the mzCloud database was used to identify compounds on MS/MS level with a mass tolerance of 10 ppm. Compound class set to All, Match ion activation type = False. Chemspider, BioCyc, and KEGG database were used to annotate features based on exact mass (MS1) with a mass tolerance of 5 ppm. An endogenous metabolites database of 4400 compounds embedded in the CD internal database also was used to annotate metabolites. 3) In addition, an in-house mass spectral library of >3,384 microbial natural products and mycotoxin reference standards was used to search by the search mass lists module in CD. ③mass tolerance set to 5ppm, RT tolerance(min)=0.05, S/N threshold=1.5 infill gaps (missing values) module. We used Xcaliber 4.0 (Thermo Fisher Scientific, USA) to manually check and identify the metabolites identified above. Genesis peak detection method was selected to detect the peak area. Minimum peak height(S/N) set to 3. Metabolic features(m/z) with peak intensity greater than 1E4 in QC raw data files were retained. The metabolites that did not meet the requirements were discarded. Finally, compounds were identified and selected by integrating positive and negative ion mode to build a peak list for subsequent biomarker analysis. For GNPS, feature-based molecular networking (FBMN)20 analysis was performed using MZmine 4 and the GNPS platform to build the molecular network. The .csv quantitative and .mgf mass spectrum profiles file (from MZmine 4 output above) and metadata file were submitted to the GNPS20 platform (http://gnps.ucsd.edu) to generate the molecular network. We used a parent mass tolerance of 0.02 Da and an MS/MS fragment ion tolerance of 0.02 Da to cluster the spectra and to subsequently create consensus spectra from which the network is constructed. All networks where edges were filtered to have a cosine score above 0.6 and more than five matched peaks. The mass spectra in the molecular network were then matched to the mass spectral libraries of GNPS with a threshold score above 0.6 and at least five matched peaks recorded. Other networking parameters were left at default levels. The obtained molecular network was then visualized using Cytoscape (version 3.9.1). For the MetDNA3 software, the metabolites features table (from MZmine 4 output above) and MS/MS data file (.mgf format) were upload and the default parameters were used. Furthermore, those metabolic features not identified by the above three software were searched in the online public databases, METLIN (https://metlin.scripps.edu) and NIST (https://www.nist.gov/). Then, the MS/MS spectra of all the above matched metabolite features was exported from Xcalibur, and these matches were manually checked to confirm the annotation, which was considered a level 2 identification. In addition, if metabolite features with MS/MS spectra were not matched in public databases as analyzed by the above tools and we considered them a MSI level 4 unknown identification. Third, to further annotate these unknown metabolic features, we construct a self-built MS1 mass spectra Aspergillus metabolites library from NPAtlas (https://www.npatlas.org/) and Reaxys (https://www.reaxys.com/). We used accurate masses (±5 ppm) to deduplicate the metabolic features using Compound Discovery 2 (Thermo Fisher, USA), and gave an MSI level 3 identification. Fourth, we used the computational metabolomics annotation tools MolDiscovery and Dereplicator+57 to annotate microbial metabolites, mainly peptides natural products. They are both implemented in the GNPS ecosystem and we used default parameters. The mzCloud, MetDNA2, DERREPLICATOR+, and MOLDISCOVERY annotation results were summarized in Supplementary Figure S6b. Mirror plots of query and GNPS mass spectral library MS/MS spectra of representative metabolites that show high to decent matching qualities (Supplementary Figure S6e-j). We uploaded 355 metabolites MS/MS spectral alignment mirror plots to demonstrate the reliability of the metabolites annotation results, which are stored in https://zenodo.org/records/18432127. MetDNA3 only outputs MS/MS similarity score results, and their database is not open source; therefore, we cannot obtain mirror plots. In addition, we provide chromatographic retention times (Supplementary Figure S27) and MS/MS chromatograms of standards and A.flavus mycelial extracts of six precursor metabolites in the aflatoxin biosynthetic pathway, including 5-methoxysterigmatocystin, Norsolorinic acid, Averantin, Averufin, Versicolorin B, and Versicolorin A (Supplementary Figure S28-35).

Fungi RNA extraction, sequencing, and transcriptome analysis
We selected 28 representative strains from different clades on the phylogenetic tree for cultivation, with the same culture conditions with our previous research24 and three biological replicates. Transcriptome sample preparation, sequencing, and data analysis were consistent with our previous study methods58. After 3 days of cultivation, we harvested A. flavus hyphae for RNA isolation. The TRIzol method was used to extract the fungi RNA. Firstly, the filamentous fungus was ground to a powder using liquid nitrogen, and the powder was transferred into the 2 mL tube contains 1.5 mL Trizol reagent. The mixture was shaken for 3 min, and placed for 5 min at normal temperature; then it was centrifuged at 10000g for 5 min. 200 μL of chloroform/isoamylalcohol (24:1) was added to the supernatant with 1 mL of lysis reagents. After centrifugation at 10000 g for 10 min, the supernatant was transferred into another new tube with an equal volume of isopropanol and then refrigerated at -20 °C for 1 hour. After centrifugation at 13600g for 20 min at 4 degrees celsius, the supernatant was precipitated by 1 mL of 75% ethanol and let dry for 5 min. The RNA pellet was dissolved in 30~100 microliters of DEPC water or RNase-free water. The concentration of the extracted RNA samples was determined using a Nanodrop system(NanoDrop, Madison, USA), and the integrity of the RNA was examined by the RNA integrity number (RIN) using an Agilent 2100 bioanalyzer (Agilent, Santa Clara, USA). The mRNA was purified with Oligo(dT)-attached magnetic beads and subsequently fragmented into small pieces with a fragment buffer at appropriate temperatures. Then, the first-strand cDNA was generated by PCR, and the second-strand cDNA was synthesized as well. The final product of the reaction was purified through end-repair and subsequent incubation. The cDNA fragments with adapters were amplified by PCR before they were purified with Ampure XP beads. The quality and quantity of the library were assessed using an Agilent 2100 bioanalyzer before it was subjected to DSN treatment. The DSN-treated library's quality was evaluated two methods to ensure high-quality sequence data. First, we checked the distribution of the fragment size using the Agilent 2100 bioanalyzer, and secondly, we quantified the library using real-time quantitative PCR (qPCR). The qualified library was amplified on cBot to generate the cluster on the flow cell, and the amplified flow cell was sequenced single-end on the Hiseq 4000. Illumina sequencing was commissioned from the Beijing Genomics Institute(BGI) for RNAseq sequencing.

After obtaining the raw sequencing data, low-quality, adaptor-polluted sequences and reads with a high content of unknown base (N) were filtered. The clean reads were mapped to an A. flavus reference genome (GCA_000006275.3) using HISAT and Bowtie 2 with default parameters. The gene expression level for each sample was determined with RSEM59. The FPKM (Fragments Per Kb of exon per Million reads) and TPM method were used to calculate and normalize the expression level of genes. The Dr.Tom (BGI) cloud platform was used to generate a transcriptome expression quantitative data matrix. Gene annotation, GO and KEGG enrichment analysis, GSEA enrichment analysis, and WGCNA analysis were also completed in this platform. Filtering, normalization, and statistical analysis of transcriptome quantification data matrix were performed in Expressanalyst60 with default parameters. Log2-counts per million (logCPM) transformation method was used to normalize. Limma61 R package was selected to do the differential expression analysis against a common control of Clade1, pairwise comparisions, and nested comparisons method in simple metadata penel. To perform multi-factor comparison analysis for complex metadata, we leveraged linear models with covariate adjustments of limma for its high-performance implementation in the Expressanalyst60 cloud platform. We executed the plotting in Windows 10 R.4.1.0 environment with the ggplot2 package (https://r-graph-gallery.com/ggplot2-package.html).

Statistics and Reproducibility
To compare the differences between the two groups, we chose the two-sided Wilcoxon test or Student's t-test, with a p-value of P < 0.05. The correction method chosen is Benjamini-Hochberg (BH). The gene pathway enrichment analysis in Figures 3c, d, e and Figures 4c, d, e, f, g was performed using the hypergeometric test. T-tests, principal component analysis (PCA), and Projections to Latent Structures Discriminant Analysis (PLS-DA) model were used to pre-screen differential metabolites between the two clades. Univariate analysis used the Student's t-test, with a p-value of P < 0.05. A variable influence on projection (VIP) value ≥1.0 was used as a selection threshold to select differential metabolites in the PLS-DA model. A statistical method specialized for multi-testing, SAM (Significance Analysis of Microarrays) to validate the credibility of differential metabolites between the two clades in metabolome-wide analysis using MetaboAnalystR62 package. The sample-wise permutation (default by the samr package) and distribution-independent ranking tests (based on the Wilcoxon test) were used to ascertain significance (false discovery rate, FDR <0.05) for the SAM analysis. The default recommended parameters were used in SAM analysis. The differential metabolites identified and screened by the steps mentioned above with MSI levels 1, 2, and 3 were pooled together. Enrichment analysis and the pathway analysis module in MetaboAnalyst 5.0 were used to perform the metabolite set enrichment analysis (MSEA) as well as metabolic pathway analysis (MetPA), respectively, using all MSI-level-1-and-2-annotated metabolites as input. All the default parameters were used in MetaboAnalyst 5.062. All statistical analyses were performed on a Windows 10 R.4.1.0 environment and using Excel 2019 (Microsoft, USA).

Gene knockout experiments
We referred to the method of Perng-Kuang Chang's laboratory63 to construct an A. flavus fluorescent engineering strain with the PyrG gene using an aflatoxigenic A.flavus strain LNZW-1A isolated from the peanut field in Zhangwu, Liaoning province. The detailed gene knockout experimental procedure is as follows: 1) Primers were designed based on the target gene sequence using Primer3 (https://www.primer3plus.com/). Primer synthesis was performed by Wuhan Tianyi Huayu Gene Technology Co., Ltd. Primer information used in A. flavus regulatory gene knockout experiments was provided in the supplementary table 19. 2) First, using the A. flavus LNZW-1 genome as a template, the 5′ upstream homologous arm and 3′ downstream homologous arm sequences (each 1000-1200 bp) of the target gene, as well as the full-length pyrG gene, were amplified (reaction system: 1 µL template DNA, 1 µL Primer-F, 1 µL Primer-R, 12.5 µL 2×Hieff Robust PCR Master Mix, and ddH2O added to a total volume of 25 µL); 3) Amplification procedure: Pre-denaturation at 94℃ for 3 min, denaturation at 94℃ for 10 sec, annealing at 58℃ for 20 sec, extension at 72℃ for 15 sec/kb; denaturation, annealing, and extension were repeated for 35 cycles, followed by a final extension at 72℃ for 5 min; PCR instrument: (Bio-RAD, T100 Thermal Cycler, USA). 4) The upstream homologous arm, the full-length pyrG gene, and the 3' downstream homologous arm fragment were fused together using the Double-Joint PCR method (reaction system: 1 µL of 5' upstream homologous arm, 3 µL of full-length pyrG, 1 µL of 3' downstream homologous arm, 5 µL of 2×Hieff Robust PCR Master Mix, and ddH2O added to a final volume of 20 µL). 5) Amplification procedure: pre-denaturation at 94°C for 2 min, denaturation at 94°C for 30 sec, annealing at 60°C for 1 min, extension at 72°C for 3 min, denaturation, annealing and extension repeated for 10 cycles, final extension at 72°C for 10 min. The genetically transformed fragment was obtained by nested PCR amplification, and the transformed fragment was recovered by gel extraction (E.Z.N.A Gel Extraction Kit, Omega Bio-tek, D2500-02). 6) Finally, the transformed fragments were introduced into A. flavus protoplasts using a polyethylene glycol (PEG)-mediated method to obtain the knockout strain. The mutant strain was then verified for successful gene knockout by amplifying the target gene and its upstream and downstream homologous arm fragments. The preparation methods and composition of the key reagents used in the experiment:
(1) Enzymatic lysis solution: 200 mg Lysing Enzymes (Sigma, USA, 37340-57-1), 50 mg Driselase (Sigma, USA, 85186-71-6), 20 µL β-glucuronidase (85000 µg/mL, Sigma, USA, 9001-45-0), 1.4 g NaCl, 8 µL of 1 mol/L CaCl2, and 10 mmol/L NaH2PO4 (pH 5.8) to a final volume of 20 mL. The lysis solution was filtered through a 0.22 µm filter after preparation. The lysis solution was prepared on ice and used immediately after preparation. (2) 1.2 mol/L STC: 1.2 M Sorbitol, 10 mM Tris-HCl (pH 7.5), 50 mM CaCl2; Weigh 21.86 g of sorbitol, add 1 ml of 1 M Tris-HCl (pH 7.5) solution and 5 ml of 1 M CaCl2 solution, and dilute to 100 ml with deionized water. Sterilize at 121℃ for 20 min. (3) 50% PEG: 50 g of polyethylene glycol (PEG4000), 1 mL of 1 M Tris-HCl (pH 7.5) solution, and 1 mL of 1 M CaCl2 solution were heated and stirred until dissolved. Deionized water was then added to a final volume of 100 mL, and the solution was sterilized at 121°C for 20 minutes. (4) Regeneration medium: 34.2 g sucrose, 0.1 g casein hydrolysate, 0.1 g yeast extract, 1.5 g agar, deionized water to a final volume of 100 mL, sterilized at 121℃ for 20 min. (5) PDA and PDB culture media: Boil 200 g of peeled potatoes, filter the liquid through four layers of gauze, add 20 g of glucose, and adjust the volume to 1 L. Sterilize in aliquots at 121°C for 20 minutes. Adding agar (15 g/L) results in PDA medium.
Protoplast preparation: A. flavus spores were inoculated into PDB medium and cultured overnight at 37°C and 200 r/min in a constant temperature shaker (Shanghai Zhicheng Analytical Instruments Manufacturing Co., Ltd., ZWY-2102C, Shanghai). The mycelia were filtered through four layers of lens paper and washed three times with sterile water. The mycelia were then transferred to lysis buffer and incubated at 30°C and 100 r/min for 4 hours. The number and morphology of protoplasts were checked every half hour using a microscope (Nikon, ECLIPSE E100, Japan) after 2 hours. The lysis buffer was filtered through four layers of lens paper into a 50 mL centrifuge tube. The tube was centrifuged at 5000 r/min for 10 minutes at 4°C in a low-temperature centrifuge (Hunan Xiangyi Laboratory Instruments Development Co., Ltd., H1850R, Hunan). The supernatant was carefully decanted, leaving a white precipitate at the bottom of the tube. 2 mL of STC was added to the opposite side of the precipitate, and the centrifuge tube was gently rotated to wash the precipitate. The washing solution was decanted. An appropriate amount of STC was added, and the precipitate was resuspended by pipetting. The suspension was centrifuged at 5000 r/min at 4°C for 10 minutes, and the supernatant was carefully decanted. An appropriate amount of STC was added, and the precipitate was resuspended by pipetting. The protoplast concentration was adjusted to 1 × 10⁶ cells/mL for further use. Add the prepared DNA fragments (total concentration > 5 ng) to a 10 mL centrifuge tube, then add pre-cooled STC buffer to a final volume of 200 µL, and mix thoroughly with a pipette. Add 100 µL of freshly prepared A. flavus protoplasts to the mixture, mix thoroughly with a pipette, and incubate on ice for 30 min. Then, sequentially add 200 µL, 400 µL, and 800 µL of PEG Buffer, gently invert to mix, and incubate on ice for 20 min. Finally, add 1 mL of pre-cooled STC buffer and incubate on ice for 10 min. Add all the transformation mixture to regeneration medium at approximately 40°C, mix thoroughly by shaking, and pour onto plates. Invert the plates and incubate in a 30°C constant temperature incubator (Shanghai Jinghong Experimental Equipment Co., Ltd., DNP-9022, Shanghai) in the dark for 36 hours. Pick the hyphal tips of single colonies and transfer them to PDA medium. After 3 days of culture, extract DNA and amplify the fragments to verify the transformants. Supplementary Figure 37-39 show the original gel images for the three gene knockout (ΔAF210, ΔAF420, and ΔAF890).

Paired omics data analysis
NPLinker, a metabologenomic analysis software framework22, was employed to link genomic and metabolomic data and connect the identified microbial specialized metabolites to their candidate biosynthetic genomic regions. First, we downloaded the molecular networking results by feature-based molecular networking (FBMN) analysis (metabolomics) of 544 strains (7 strains with low data quality were filtered) from the GNPS platform20. Metabolome data, antiSMASH-predicted BGCs, and BiG-SCAPE-derived GCFs were compiled as input for NPLinker. In the configuration file, we set the Metcalf threshold to 4 (the default value is 2) in order to speed up NPLinker and filter out links with a low Metcalf association score. Other parameters were default. After running NPLinker, the results were examined for promising BGC-mass spectral links in the NPLinker web app running from a Docker instance.

Genotype-Environmental Analysis
Sambada v0.5.364, based on logistic regression and Bayenv 265 by bayes factor analysis models were recruited to identify putatively adaptive loci associated with environmental factors. A single SNP locus was identified as a candidate locus when the log-likelihood ratios (G scores > 60) were significant. G‐Scores were corrected for multiple testing with the Bonferroni method at a 95% confidence level. In addition, Bayenv2, a covariance matrix based on putatively neutral markers was used to assess the correlations between SNP and environmental variables at the Markov chain Monte Carlo (MCMC) model. Bayesian factor(BFs) was used as an indicator to evaluate the degree of genetic and environmental association. An averaged log10(BF) value >1.2 is considered a threshold support for the model where parameters significantly affect allele frequencies.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper describes genome sequencing of a large number of Aspergillus flavus isolates from across China and a detailed comparative genome analysis of this species through comparing the new sequences with earlier sequenced genomes. The authors also perform a detailed metabolome analysis of samples obtained from the isolates and the study represents the most comprehensive analysis of a single fungal species known today! Considering the importance of A. flavus as a pathogen and a known producer of mycotoxins the study is extremely valuable, and a particular interesting finding is the detailed analysis of variations in secondary metabolite production among the individual strains due to both genetic variations as well as mutations in genes of the biosynthetic gene clusters. The paper reads very well and presents the results in a clear fashion and I do not have any comments for improvements. I would just like to congratulate the authors with such a comprehensive paper.", "reviewer_2": "This study presents a large-scale, multi-omics analysis of Aspergillus flavus, integrating genomics, transcriptomics, and untargeted mass spectrometry-based metabolomics to investigate phylogeographic patterns and regulatory drivers of specialized metabolism. The manuscript provides valuable insights into fungal secondary metabolite diversity and proposes connections between biosynthetic gene clusters (BGCs), environmental origin, and metabolic output. Given the growing interest in microbial metabolomics in ecological and evolutionary research, the work is timely and potentially impactful. However, several critical aspects of the mass spectrometry methodology (particularly related to signal-to-noise thresholds, gap filling practices, metabolite identification confidence, and dataset merging) require clarification and improvement to enhance the reproducibility and interpretability of the results.\n\nMajor Revisions\nL910. The authors report on using an MS2 noise threshold of 100 intensity units (1E2) in the MZmine feature processing workflow. This threshold is unusually permissive for high-resolution Orbitrap data, where baseline noise in complex biological matrices (e.g., fungal extracts) typically exceeds this value. The inclusion of such low-intensity fragment ions increases the likelihood of incorporating noise-derived or non-informative peaks, which can degrade the quality of MS/MS spectra.\nThis setting undermines the reliability of spectral similarity scores and increases the risk of false-positive annotations, spurious molecular networking connections, and diminished reproducibility. The authors are strongly encouraged to reprocess the data using a more stringent MS2 noise threshold (≥1E3–1E4), as recommended for Orbitrap-based untargeted metabolomics (Mahieu and Patti, 2017).\nL920. The use of gap filling to address missing features further compounds the concern raised above. Although gap filling is widely accepted to mitigate missing values caused by stochastic MS2 sampling in DDA workflows, it can introduce false-positive intensities, particularly when applied to features near the noise threshold. In this study, the permissive MS2 threshold amplifies this risk, increasing the probability that baseline fluctuations or background ions are retrospectively integrated as genuine peaks. The authors should clarify how gap-filled features were treated in downstream analyses and whether any filtering based on reproducibility or detection frequency was applied.\nL924. The manuscript describes the fusion of MS1-only data from Drott et al. with the authors’ own MS1/MS2-rich dataset. This reviewer finds this practice debatable. Without MS/MS data, no structural confirmation is possible, and without retention time alignment or adduct correction, accurate mass alone is insufficient for robust feature matching.\nThis approach violates standard and well-published practices in cross-dataset untargeted metabolomics integration and raises significant concerns about false feature merging and annotation bias (Dunn et al., 2011; Broadhurst et al., 2018). The authors should, at minimum, report the m/z tolerance window used for matching, whether adduct species were deconvoluted, and how annotation uncertainty was handled during matrix fusion.\nL932. The authors state that metabolite identification was performed using an in-house mycotoxin library of over 260 compounds, with confirmation based on accurate mass (±5 ppm), MS/MS spectra, retention time, and isotope pattern. However, no evidence is provided (such as GNPS molecular networking results, mirror plots, or cosine similarity scores) to support the claim of Level 1 or Level 2 identifications.\nThis omission limits the reviewer’s ability to verify structural assignments. The authors should provide: GNPS job links (not currently listed in L945–953), representative mirror plots showing MS/MS spectral matches between standards and samples, and cosine similarity metrics for key annotated compounds. Without these supporting materials, the claim of confirmed identification using authentic standards remains unsubstantiated.\n\nMinor comments:\nL853. Change polarized metabolites to polar metabolites.\nL855. Define “proper amount”, is there any range of beads used to sample homogenization?\nL863. Did the author use extraction blanks? The authors state that only methanol is used as a solvent blank.\nL870-873. Use names instead of molecular formulas to define water, methanol, formic acid and ammonium formate.\nL873-877. Gradient must be addressed based on organic mobile phase (B) rather than water (A). This section needs to be corrected.\nL1077. Change porsitive for positive.\nL1079. Massive data is not available.\nInclude any potential limitations related to mass spectrometry and metabolomics data analysis in the discussion section.\nAlso, the authors should ensure that mass spectrometry-related terminology is consistently and correctly formatted throughout the manuscript. For example, the mass-to-charge ratio should be written as m/z (in italics and with a lowercase \"m\") following standard conventions. Please review the text to correct any inconsistencies in the formatting of this and related terms.\n\nConclusion\nThe integration of omics layers in this study is conceptually strong, and the manuscript addresses a highly relevant biological question. However, multiple issues in the mass spectrometry data processing pipeline and metabolite annotation framework require careful revision. Addressing these concerns will significantly strengthen the manuscript’s technical foundation and improve its reproducibility and impact.", "reviewer_3": "The authors have sequenced the genomes of 544 strains of the fungus Aspergillus flavus. This large dataset was combined with previously published genomes, as well as metabolomics and transcriptomics data. Using various computational approaches, new insights were obtained regarding evolution of biosynthetic gene clusters. Of most interest, the data was used to identify novel genes putatively involved in various aspects of secondary metabolism, and this was confirmed using gene knockout and metabolic phenotyping.\n\nThe manuscript is generally well written, but aspects of the data presentation should be considerably improved. Specifically, there are a lot of figures that are unclear, either due to tiny fonts, mislabeling, incorrect ordering, unlabeled axes, or incomplete explanations for the average reader. See specific comments below.\n\nThe pangenome was determined based on “567 representative genomes” (line 252). Unless I missed it, it is not explained in the corresponding section of the material and methods how these 567 genomes were representatively selected. This has implications for the use of the rarefaction curve in Figure 3b, for example. If the 567 genomes comprised the 544 newly-sequenced Chinese genomes plus some additionals, then this is hardly a representative set.\n\nFigure 3d shows the enrichment of KEGG terms in the accessory gene families. Enrichment suggests that there is over-representation in the accessory genes compared to something else (core, or all other genes?). The figure appears to just show gene counts, though. Not enrichment.\n\nThe metabolomic analysis was performed in only 1 condition, so it could mean that some BGCs were not active under the tested condition (and likely has an effect on the data presented in Figure 2bcd). For such a large set of strains it is not required to test additional conditions, but it could be discussed somewhere.\n\nLines 331-332. Figure 7 is referenced before even figure 4. Probably the reference is wrong, because I don’t see how the authors draw these important conclusions from this figure (nor can I find the actual figure that they likely intended to reference).\n\nFigure 4 is not easy to interpret. The axes in 4b and 4c are unlabeled, making it very unclear what is depicted here and what conclusions can be drawn. It is very unclear to me what is depicted in 4ghi. Note that in line 467 the panels 4ghi are erroneously referred to as 4fgh, adding to the confusion.\n\nFigure 5. The labels for 5e and 5d are swapped…\n\nLine 431. “dispatched” is an odd word choice. Maybe just “used”?\n\nLines 196-198. “We observed that the genomes of these A. favus strains (W-clade 8) harboured homologues of 156 virulence genes….” What about the other strains? And were these homologues of virulence genes over-represented in the clinical strains versus the environmental strains?\n\nIn lines 142-177 (and elsewhere, including figures). There is inconsistency in the naming of the populations: population B, PopB, pop B\n\nSome of the text in the figures (e.g. Figure 2) is way too small, even in the original files. Especially in Figure 2b.\n\nLine 313. Typos “In addition, by”\nLine 488 Typo “functional”\nLine 584. Typo “and”\nLine 595. Typo “environmental”"}, "author_response": "Response to Reviewer #1\nReviewer #1 comments:\nThis paper describes genome sequencing of a large number of Aspergillus flavus isolates from\nacross China and a detailed comparative genome analysis of this species through comparing the\nnew sequences with earlier sequenced genomes. The authors also perform a detailed metabolome\nanalysis of samples obtained from the isolates and the study represents the most comprehensive\nanalysis of a single fungal species known today! Considering the importance of A. flavus as a\npathogen and a known producer of mycotoxins the study is extremely valuable, and a particular\ninteresting finding is the detailed analysis of variations in secondary metabolite production among\nthe individual strains due to both genetic variations as well as mutations in genes of the\nbiosynthetic gene clusters. The paper reads very well and presents the results in a clear fashion\nand I do not have any comments for improvements. I would just like to congratulate the authors\nwith such a comprehensive paper.\nResponse: We would like to express our sincere gratitude to the reviewer for their positive and\nencouraging feedback.\n\nResponse to Reviewer #2\nReviewer #2 comments:\n1. This study presents a large-scale, multi-omics analysis of Aspergillus flavus, integrating\ngenomics, transcriptomics, and untargeted mass spectrometry-based metabolomics to investigate\nphylogeographic patterns and regulatory drivers of specialized metabolism. The manuscript\nprovides valuable insights into fungal secondary metabolite diversity and proposes connections\nbetween biosynthetic gene clusters (BGCs), environmental origin, and metabolic output. Given\nthe growing interest in microbial metabolomics in ecological and evolutionary research, the work\nis timely and potentially impactful.\nResponse: We would like to express our sincere gratitude for your constructive comments, which\nhighlighted the significant aspects of our research. Your comments have been encouraging and have\nmotivated us to continue our efforts to improve our study.\n2.However, several critical aspects of the mass spectrometry methodology (particularly related to\nsignal-to-noise thresholds, gap filling practices, metabolite identification confidence, and dataset\nmerging) require clarification and improvement to enhance the reproducibility and\ninterpretability of the results. Major Revisions.\nResponse: Thank you for your very helpful comments, which have enabled us to significantly improve\nthe quality of the manuscript. We agree with the reviewers’ comments and requests and have added\nanalyses and results to support our feature finding, annotation workflow, and results. We clarified and\nimproved the key deficiencies in the mass spectrometry method raised by the reviewers. These indicated\nissues have been addressed in accordance with the reviewers' comments, and we respond to the points\none by one in detail below.\nL910. The authors report on using an MS2 noise threshold of 100 intensity units (1E2) in the\nMZmine feature processing workflow. This threshold is unusually permissive for high-resolution\nOrbitrap data, where baseline noise in complex biological matrices (e.g., fungal extracts) typically\nexceeds this value. The inclusion of such low-intensity fragment ions increases the likelihood of\nincorporating noise-derived or non-informative peaks, which can degrade the quality of MS/MS\nspectra. This setting undermines the reliability of spectral similarity scores and increases the risk\nof false-positive annotations, spurious molecular networking connections, and diminished\nreproducibility. The authors are strongly encouraged to reprocess the data using a more stringent\nMS2 noise threshold (≥1E3–1E4), as recommended for Orbitrap-based untargeted metabolomics\n(Mahieu and Patti, 2017).\nResponse: Thanks to the reviewer for raising this point. We agree with the reviewer's comments and\nrecognize the concern that too low a threshold for MS2 noise may lead to inaccurate metabolite\nidentification. We follow the reviewers' suggestions and raised the MS2 noise threshold from E2 to 6E3.\nTo better justify the chosen parameter settings, the final MS2 noise threshold (6E3) selection was\ndetermined by comparing the number of core and unique metabolites identified with the GNPS platform\nthrough Feature-Based Molecular Networking (FBMN) analysis across eight preset MS2 thresholds\n(1E0; 1E2; 1E3; 2E3; 4E3; 6E3; 8E3; 1E4) in the MS2 peak detection step, mostly falling within the\nrecommended Orbitrap-based untargeted metabolomics MS2 noise threshold range (≥1E3–1E4)1. To\nfacilitate reviewers in viewing the MS/MS spectra and quantitative matrix tables extracted at different\nMS2 noise thresholds, we have uploaded analysis results to Google Drive\n(https://drive.google.com/drive/folders/1dD04vMPQ8pYAHq87mubIPPowUvpJdn-\nh?dmr=1&ec=wgc-drive-globalnav-goto). In this parameter optimization analysis, except for the MS2\nnoise threshold, all other parameters were kept equal. The metabolome batch processing workflow for\n544 samples in MZmine 4 is shown in Re_Figure 1a. We used optimized mass spectrometry acquisition\nparameters in the data acquisition method2, including the AGC value, the high-energy collision-induced\ndissociation (HCD). The intensity filtering parameter (MS1=1E5, in this study) between the MS1\nacquisition and the MS2 acquisition part was set. Only those precursor ions with MS1 features greater\nthan the threshold will be fragmented into MS2 product ions by the high-energy collision-induced\ndissociation (HCD). Thanks to the high selectivity and sensitivity of HPLC-Orbitrap Fusion-MS, we\nobtained relatively clean MS/MS spectra. In the FBMN analysis, the cosine similarity score threshold\nwas set to 0.6, and the minimum number of MS2 spectra matching data points between the sample and\nreference spectra was 6. Other detailed analysis procedures and parameters in mzmine 4.29 are\nexplained in the Methods section.\nRe_Figure 1| a. Mzmine 4 analyzes batch processing in a streamlined analysis mode, and the gap-\nfilling step includes four filtering steps: the intensity tolerance (20%), mz/ tolerance (sample-to-sample),\nretention time, and minimum scans (data points) filtering steps.\nWe observed a decrease in metabolic features containing MS/MS spectra with increasing MS2\nnoise threshold (Re_Figure 2a). For example, we observe a decrease of annotated metabolites (n=371-\n335=36) when increasing the MS2 noise threshold from 1E2 to 1E4, and as the MS2 noise threshold\nincreases, the number of identified metabolites further decreases (Re_Figure 2b). In particular, the\npreviously set MS2 noise threshold (1E2) had ~15% ((371-317)/371=~15%) identified metabolites that\nmay carry a high risk of false positives, as these were not consequently annotated (Re_Figure 2c). 317\nannotated metabolites were identified under all MS2 noise threshold conditions without gap-filling\ntreatment, accounting for ~85% (317/371). Based on the above threshold comparison analysis results,\nwe made a trade-off between the number of identified metabolites and reliability and decided to use the\n6E3 threshold.\nRe_Figure 2| a. The number of extracted metabolic features/MS2 spectra by Feature-Based Molecular\nNetworking (FBMN) analysis across seven preset MS2 thresholds, with or without a gap-filling step. b.\nThe number of identified metabolites by Feature-Based Molecular Networking (FBMN) analysis across\nseven preset MS2 thresholds with or without a gap-filling steps. c. Venn diagram showing the\nintersection and union of metabolites identified by FBMN analysis at seven MS2 noise threshold\nsettings without gap-filling step. c. Venn diagram showing the intersection and union of metabolites\nidentified by FBMN analysis at seven MS2 noise threshold settings with a gap-filling step.\nAfter determining the most critical MS2 noise threshold in metabolite identification, we identified\nmetabolite MS/MS spectra by aligning with the MS/MS spectra in reference database based on three\nmetabolite identification platforms (annotation Level 2). These platforms include: (1) FBMN from the\nGNPS metabolomics analysis environment. The FBMN analysis results after increasing the MS/MS\nextraction threshold to 6E3 were performed on the GNPS platform and the results were stored in GNPS\nwebsite (https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=1d02472cebf54e63999ee7c0d7891a5d); (2)\nCompound Discovery, an online or offline metabolomics analysis software from the instrument\nmanufacturer Thermo Fisher, which compares to mzCloud, a database built using 32,330 standards by\nThermo Fisher (Data as of October 2025 for https://www.mzcloud.org/). The metabolite identification\nanalysis results of Compound Discovery are listed in the supplementary material table; and (3)\nMetDNA3, which primarily identifies primary metabolites in this study. The annotation results of\nMetDNA3 are presented in the supplementary material table 10; Furthermore, on top of these large mass\nspectral libraries, for the aflatoxin biosynthesis pathway, representing key metabolites of interest in this\nstudy, we previously developed a targeted metabolite monitoring method. In our previous study, 12\nMS/MS spectra of standards for compounds in the aflatoxin biosynthesis pathway were collected3. We\nuploaded 355 metabolites MS/MS spectral alignment mirror plots to demonstrate the reliability of the\nmetabolites annotation results through corresponding matches of higher-abundant mass fragments\ntogether with high cosine similarity scores of the spectral match. The mirror plots are stored in\nhttps://figshare.com/articles/dataset/_b_554AFs_1E56E3_no_gap-filling_Mirror_plot_b_/30434335.\nSince Figshare website does not provide a convenient way to view mirror plots, we have selected eight\nmetabolites and displayed their mirror plots below (Re_Figure 3).\nRe_Figure 3| The mirror plots of 8 metabolites. The black portion at the top is the MS/MS spectrum\nfrom the sample, and the green portion at the bottom is the MS/MS spectrum from the GNPS library.\nThe mirror graph includes the cosine similarity score.\nL920. The use of gap filling to address missing features further compounds the concern raised\nabove. Although gap filling is widely accepted to mitigate missing values caused by stochastic MS2\nsampling in DDA workflows, it can introduce false-positive intensities, particularly when applied\nto features near the noise threshold. In this study, the permissive MS2 threshold amplifies this\nrisk, increasing the probability that baseline fluctuations or background ions are retrospectively\nintegrated as genuine peaks. The authors should clarify how gap-filled features were treated in\ndownstream analyses and whether any filtering based on reproducibility or detection frequency\nwas applied.\nResponse: We agree with the reviewer's comments and appreciate the suggestions to improve the\nquality of our manuscript. As the reviewer was concerned about, using gap filling with a low MS2 noise\nthreshold to address missing features exacerbates the issues raised regarding false-positive annotations,\nspurious molecular networking connections, and reduced reproducibility (impacting qualitative analysis\nresults), as well as leading to false-positive intensities (impacting quantitative analysis results). As\nindicated above, following the reviewer's suggestion, we increased the MS2 noise threshold to 6E3 to\naddress these concerns. We evaluated the effects of a gap-filling step on the qualitative and quantitative\nmetabolite analysis results. We will next describe the comparative evaluation results in detail, clarifying\nhow gap-filled features were treated in downstream analyses and the frequency with which they were\ndetected.\nFirstly, we added a set of analyses on the effect of gap-filling on metabolite qualitative analysis\nwhen performing MS2 noise threshold parameter optimization analysis experiments. We found that the\nMS2 noise threshold should be increased to at least 2E3 to affect the change of the MS/MS spectrum\nfor this dataset. The gap-filling treatment reduced the number of MS1 features, MS/MS spectra, and\nidentified metabolites, and as the MS2 noise threshold increased, the trend of the number of identified\nmetabolites was similar to that without the gap-filling step (Re_Figure 2a,b). Importantly, the difference\nin the number of MS/MS spectra reduced by adding the gap-filling step is much greater than the\ndifference in the number of MS/MS spectra reduced by changing the MS2 noise threshold (Re_Figure\n2a,b,c,d). We also observed that almost all metabolites identified (99%) in the gap-filling treatment\nwere also identified in the treatment without gap-filling (Re_Figure 4). In the export step parameter\nsettings, we only export metabolic features with MS2, so the number of MS1 features and the number\nof MS/MS spectra are the same (Re_Figure 2a). We speculate that adding a gap-filling step increases\nthe intensity tolerance, mz/tolerance (sample-to-sample), retention time, and minimum scans (data\npoints) filtering steps (Re_Figure 1a), thereby reducing the number of MS/MS spectra we obtain and\nthus leading to a decrease in the number of identified metabolites. So, adding a gap-filling step increases\nthe data quality of the MS/MS spectra, but it may also filter out the specialized metabolites, such as\nAflatoxin G2, Aflatoxin M1, and Aflatrem, that are produced by a small number of strains we are\ninterested in this study (Re_Figure 4). In summary, it is a matter of choice: with additional filters in\n“Feature finder” steps, we may (and do) end up with fewer features, but they may be of higher data\nquality for the feature matrix. Also, the gap filling may help to get better statistics downstream, but it\nmay also yield artifacts due to the added minimal values (Re_Figure 2a,b). We upload the MS/MS\nalignment mirror plot of these metabolites in\nhttps://figshare.com/articles/dataset/_b_554AFs_1E56E3_no_gap-filling_Mirror_plot_b_/30434335.\nRe_Figure 4| Adding a gap-filling step can lead to the loss of useful MS/MS spectra information,\nthereby reducing the number of metabolites of interest identified, for example, Aflatrem.\nIn addition, a dilemma arises from comparing the Aspergillus metabolites in the Natural Products\nAtlas (https://www.npatlas.org/), CMMAspergillusDB\n(https://github.com/albertogilf/ceuMassMediator/tree/master/CMMAspergillusDB), and Reaxys\n(https://www.reaxys.com), three microbial natural product databases, with the current microbial natural\nproduct platform, GNPS, which utilizes MS/MS spectra. It was found that a large proportion of\nAspergillus metabolites lack MS/MS spectra. This situation has been reflected in several other studies\nof fungal secondary metabolites. For example, in the recent article on an A. flavus population published\nin PNAS4, putative MS/MS spectra of 13 secondary metabolites were provided in the supplementary\nmaterials. However, only 4 metabolites in this study belonged to the annotation level 2 (annotated by\nGNPS), and 9 metabolites did not show the comparison results of the MS/MS spectra of the samples\nand the MS/MS spectra of the standards/ reference database (mirror plot). Instead, they were identified\nby searching the MS1 accurate mass (±5 ppm) and isotope pattern information in the PubChem database\n(annotation level 3). PubChem lacked MS/MS spectral information for these metabolites. Another effort\nrecently studied the diversity of the cyclopiazonic acid family of mycotoxins in A. flavus by UHPLC\nTriple-TOF HRMS5. In this study, the authors collected 30 CPA-type alkaloids previously isolated and\nidentified from Aspergillus genus strains by natural product researchers. They compiled their molecular\nformulas and exact mass (MS1, ±5 ppm) data and then identified 22 members of the cyclopiazonic acid\nfamily from the A. flavus strains extracts using UHPLC Triple-TOF HRMS. The authors also provided\nMS/MS spectra and putative structural fragments for 22 detected secondary metabolites. We decided to\nuse a similar approach as followed by other natural product researchers and used MS1 (±3 ppm) and\nspecies origin information to annotate by matching an in-house Aspergillus metabolite list (n= 6796)\nextracted from CMMAspergillusDB6, NPAtlas7, and Reaxys. In addition, some MS/MS spectra also\nexist in the literature, but they are not included in the GNPS library. We also performed in-silico\nannotation of the MS/MS spectra through Sirius8,9, MolDiscovery10and DEREPLICATOR+11, making\nuse of the assembled metabolite list. Finally, we compared MS/MS spectra in the literature with the\ncharacteristic ions in the MS/MS spectra of our samples for metabolites annotation. We define these\nidentified metabolites as annotation level 3 in this study. Although they are not as reliable as level 2\nmetabolites, we believe they still have a sufficient degree of credibility to make them worth reporting.\nTools for performing MS1 matching tasks can prepare in-house databases, and then the tools of MZmine\n412, TidyMass13, and Compound Discovery 3 (Thermo Fisher, USA) be used to annotate experimental\ndata. For example, the following is a schematic diagram of the process of performing MS1 metabolite\nannotation using the metID component in TidyMass (Re_Figure 5). We then verified the metabolites\nwith annotation level 3 in Xcalibur to ensure their presence in the quality control samples and strain\nsamples. Finally, 195 specialized (secondary) metabolites (annotation level 3) were mapped using\nCompound Discovery 3 and an in-house Aspergillus metabolites list (n= 4799), and we then manually\ninspected them in Xcalibur 4.0 and selected 46 annotation level 3 and 14 annotation level 2 from GNPS\nto study their differential abundance across samples. This metabolite abundance matrix is provided in\nthe Supplementary Table 15, and the MS/MS spectra of these metabolites are provided in the\nSupplementary Material Figures 18-26.\nRe_Figure 5| A schematic diagram of the process of performing three metabolite annotation levels using\nthe metID component in TidyMass13.\nIn terms of quantification, after weighing the pros and cons, we decided to perform quantification\nin two parts. We set the MS2 noise threshold at 6E3 and conduct metabolome analysis to obtain a gap-\nfilled features matrix for downstream analysis (Supplementary Materials Table 14). For the part of\nmetabolites identified by GNPS, MetDNA3 and Compound Discovery (n=306+87+85=478, >= level\n2), a 30% detection frequency filter was applied. It is well known that an 80% detection frequency filter\nrule is popular in focusing on intergroup differences in metabolomics studies, for example, clinical\nmetabolomics, where researchers expect more similar samples/sample matrices to ensure consistency\nbetween biological replicates. However, given that we are also interested in specialized metabolites that\nare produced by a small number of strains from different ecological environments, we put this filter\nlower, i.e., 30% at the beginning of the current study. After this detection frequency filtering (>=30%\nsamples) and manual deduplication, 260 metabolites remained (>= level 2) for subsequent multi-group\ndifferential analysis to find metabolites specific to different populations. However, among these 260\nmetabolites (>= level 2), only a small number of specialized metabolites (n=6) isolated from A. flavus\ncan be matched one-to-one with the predicted and annotated biosynthetic gene cluster (BGCs) of the A.\nflavus genome (predicted secondary metabolites) by antiSMASH (fungal version) with matches to the\nMIBiG database and 16 secondary metabolite biosynthetic pathways14 that have been experimentally\nelucidated in A. flavus. These are Aflatoxin B , Aflatoxin B , Cyclopiazonic acid, O-\n1 2\nMethylsterigmatocystin, Sterigmatocystin, and Ferrichrysin. However, although we identified Aflatoxin\nG1, Aflatoxin G2, Alternariol, Penicillin G, and Aflatrem using GNPS (without gap-filling treatment),\nthey were filtered out due to their low occurrence in the samples (<30% samples). It can be seen that\nusing only gap-filling treatment will lead to the loss of some secondary metabolite information of\ninterest. For annotation levels 2 and 3 of A. flavus secondary metabolites with a detection frequency of\nless than 30%, we manually semi-quantified the target metabolites' intensity using the quantitation\nbrowser in the instrument manufacturer's Xcalibur 4.0 software (similar to a targeted metabolomics\nanalysis strategy). This peak intensity detection method is more accurate than high-throughput non-\ntargeted peak detection, such as the one used in MZmine. If we find that the starting point or end point\nof some mass spectrometry peaks is not calculated correctly, we will manually select the starting and\nending points of each mass spectrometry peak to capture a more accurate peak area. This metabolite\nintensity matrix is provided in the supplementary table 15.\nNext, according to the input file format requirements of MetaboAnalyst, we added the sample\nmetadata information to the metabolite peak table matrix. 260 metabolite matrices filtered by 30%\ndetection frequency (MS2 noise threshold=6E3 and with gap-filling treatment, annotation level 2) and\n60 manually semi-quantitative secondary metabolite matrices (MS2 noise threshold=6E3 and without\ngap-filling treatment, annotation level 2 or 3) will be imported into the metabolomics statistical analysis\nplatform MetaboAnalyst15 to perform metabolomics data matrix preprocessing and multivariate\nstatistical analysis. Detailed data preprocessing steps: 1. The data filtering step is omitted because we\nhave already filtered the 260 metabolite data matrices. The 60 manually semi-quantitative secondary\nmetabolites matrices do not require filtering. 2. The Normalization step includes: . Sample\nnormalization, select Normalization by sum; . Data transformation, select Log transformation (base\n10); . Data scaling, select Pareto scaling. The data normalization results are shown in Re-Figure 6.\nThe data normalization results show that the intensities of 260 annotated metabolites changed from a\nskewed distribution before normalization to a nearly normal distribution after normalization. 3.\nSubsequently, multivariate statistical analysis has been performed. We chose the PLS-DA differential\nanalysis algorithm to perform multi-group comparative analysis of different clades to determine the\ndifferentiation patterns of different metabolites. We compared the intensity distribution of Aflatoxin B\n1\nin different clades in the gap-filled data matrix with an MS2 noise threshold of 6E3 (Re-Figure 7ab)\nand the distribution of Aflatoxin B detected manually using Xcalibur semi-quantitative analysis (Re-\n2, 4, 5, and 6 was high, especially Chinese clade 4, which was the highest and a high-toxin-producing\nsubgroup. The two data processing methods had high consistency.\nFinally, the secondary metabolites data matrix (supplementary table 15) used in the modified\nFigure 2b, Figure 4a, and Figure 6a of the manuscript was manually checked using Xcalibur to ensure\nthat the metabolites had reliable intensity values for semi-quantification. The metabolite intensity data\nmatrix used in Figure 7 of the manuscript uses the gap-filled data matrix (supplementary table 14) was\npreprocessed and then subjected to multivariate statistical analysis, for example, PLS-DA.\nRe_Figure 6| The data normalization results of 260 metabolite data matrices.\nRe_Figure 7| a. Comparative analysis of different clades using the unnormalized intensity data of\nAflatoxin B from 260 annotation level 2 metabolites data matrices processed with an MS2 noise\n1\nthreshold of 6E3 and gap filling treatment. b. Comparative analysis of different clades using the\nnormalized intensity data of Aflatoxin B from 260 annotation level 2 metabolites data matrices\n1\nprocessed with an MS2 noise threshold of 6E3 and gap filling treatment. c. Comparative analysis of\ndifferent clades using the unnormalized intensity data of Aflatoxin B from 60 annotation level 2 or 3\n1\nsecondary metabolites data matrices processed with a manually checked secondary metabolites matrix.\nd. Comparative analysis of different clades using the normalized intensity data of Aflatoxin B from 60\n1\nannotation level 2 or 3 secondary metabolites data matrices processed with a manually checked\nsecondary metabolites matrix.\nL924. The manuscript describes the fusion of MS1-only data from Drott et al. with the authors’\nown MS1/MS2-rich dataset. This reviewer finds this practice debatable. Without MS/MS data, no\nstructural confirmation is possible, and without retention time alignment or adduct correction,\naccurate mass alone is insufficient for robust feature matching. This approach violates standard\nand well-published practices in cross-dataset untargeted metabolomics integration and raises\nsignificant concerns about false feature merging and annotation bias (Dunn et al., 2011;\nBroadhurst et al., 2018). The authors should, at minimum, report the m/z tolerance window used\nfor matching, whether adduct species were deconvoluted, and how annotation uncertainty was\nhandled during matrix fusion.\nResponse: We are grateful to the reviewer for pointing out this issue. The dataset created by our co-\nauthor, Dr. Drott, was reanalyzed according to the reviewer's comments, and the two datasets were fused\naccording to the minimum requirements1,16 recommended by the reviewer. Previously, we only utilized\nquantitative data containing MS1, but did not include the files containing MS/MS spectrum information\nfor qualitative analysis. We have now uploaded the 24 qualitative files containing the MS/MS mass\nspectra to the GNPS repository. The updated LC-MS/MS 24 qualitative files (12 positive and 12\nnegative files) are now available through the MassIVE repository (ID no. MSV000087134,\nhttps://doi.org/doi:10.25345/C54226). They used a target method of PRM mode acquired the data on\nthe Thermo Scientific Q Exactive Orbitrap mass spectrometer. For the 94 USA strain samples, MS1\ninformation had been acquired for quantification, and then qualitative files were collected to acquire the\nMS/MS spectral information for a targeted set of compounds. The previous PNAS-published effort\nstudied 13 experimentally verified pathways for specialized (or secondary) metabolites, but did not\nstudy those pathways that have been isolated and identified in A. flavus but whose corresponding\nunknown BGCs have not yet been resolved. The identification of 13 metabolites was completed by Tom\nás and Drott and they provided MS/MS spectra information of 13 secondary metabolites in the PNAS\narticle4 (Table S4 and SI Appendix, Figs. S5–S29 in this article). Four of the specialized metabolites\n(SMs) belong to annotation level 2, and the other nine metabolites belong to annotation level 3.\nPreviously, we merged the metabolite quantification table with our data matrix based on the MS1\naccurate mass (±5 ppm, m/z tolerance window) (annotation level 3). However, this data fusion method\ncomes with risks, as pointed out by the reviewers. Therefore, we performed MS1 and MS2 detection\nusing the mzmine 4 workflow and FBMN analysis on the updated files containing MS/MS data with\nthe same parameters, aiming to merge the two datasets based on MS1 accurate mass (±5 ppm, m/z\ntolerance window) and MS/MS spectrum comparison results (annotation level 2). We only retained the\nMS1 metabolic features that had MS/MS spectra. From the mzmine analysis results, we can see that\nDrott et.al’s dataset contains MS/MS spectral information on hundreds of metabolites, which is\nconsistent with their goal of targeting the analysis of a few SMs (Re_Figure 8). However, compared\nwith the number of MS1 spectra reported in their PNAS article4 (Fig 2A in this article), the number of\nmetabolic features (MS1) containing MS/MS spectra is small. Therefore, it is difficult to perform FBMN\njoint analysis with our dataset, rich in MS/MS spectra information. We therefore performed a separate\nFBMN analysis and only matched 3 metabolites with the GNPS reference library, including aflatoxin\nB , cyclopiazonic acid, and aflatrem (GNPS job result link:\n1\nhttps://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=3532f01f1f8a47c9b2e36eca093b89bc). This is\nconsistent with the results of annotation level 2 shown in Table S4 in their PNAS4 article (annotation\nlevel 2), except for the metabolite Ditryptophenaline. Furthermore, we did not find the\nDitryptophenaline MS/MS spectrum in the GNPS library. Although we identified aflatrem in our dataset,\nthe metabolite was present at low concentration levels (MS1 peak intensity less than 1E5) under the\nculture conditions we used, and it was not detected in most samples. The dataset created by Drott used\ntwo different culture media. In order to comply with the minimum standards proposed by the reviewers,\nwe therefore selected two A. flavus secondary metabolites including Aflatoxin B (cosine similarity\n1\nscore=0.99) and Cyclopiazonic acid (cosine similarity score=0.96) for data fusion (±5ppm, m/z\ntolerance window, and with MS/MS spectrum alignment mirror plot ) to verify that the presence/absence\nvariation of the BGC gene in the strains from China and the United States in our manuscript can only\npartially explain the changes in specialized metabolites. The results still support our conclusions\n(updated Figure 4a and Figure S8). The fact that the pattern of BGC variation can only partially explain\nthe changes in SMs has thus been demonstrated in the two BGC-SM pairs for aflatoxin B and\nRe_Figure 8| Metabolic features and number of MS/MS spectra extracted from mass spectrometry data\nfrom Drott’s database under different MS1 and MS2 noise thresholds.\nL932. The authors state that metabolite identification was performed using an in-house mycotoxin\nlibrary of over 260 compounds, with confirmation based on accurate mass (±5 ppm), MS/MS\nspectra, retention time, and isotope pattern. However, no evidence is provided (such as GNPS\nmolecular networking results, mirror plots, or cosine similarity scores) to support the claim of\nLevel 1 or Level 2 identifications. This omission limits the reviewer’s ability to verify structural\nassignments. The authors should provide: GNPS job links (not currently listed in L945–953),\nrepresentative mirror plots showing MS/MS spectral matches between standards and samples,\nand cosine similarity metrics for key annotated compounds. Without these supporting materials,\nthe claim of confirmed identification using authentic standards remains unsubstantiated.\nResponse: We thank the reviewer for pointing out this key omission. As requested by the reviewers, we\nprovide the GNPS job links (6E3 MS2 noise threshold results:\nhttps://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=9bcbcfc1ba7945a887c8b4a418265b97) and mirror\nplots generated by GNPS in our dataset in the manuscript (The figshare links:\nhttps://figshare.com/articles/dataset/_b_554AFs_1E56E3_no_gap-filling_Mirror_plot_b_/30434335).\nThe in-house mycotoxin database contains the following mycotoxins produced by the species\nAspergillus flavus: Aflatoxin B , Aflatoxin B , Aflatoxin G , Aflatoxin G , Aflatoxin M , Aflatoxin M,\n1 2 1 2 1 2\nO-Methylsterigmatocystin, 11-hydroxy O-methylsterigmatocystin (HOMST), Sterigmatocystin, 5-\nmethoxysterigmatocystin, Versicolorin A, Versicolorin B, versiconol, Tenuazonic acid,\nDitryptophenaline, and Kojic acid. We selected four aflatoxin pathway-related metabolites and\npresented their MS/MS alignment results (Re-Figure 9, Level 1). Our laboratory is a national biotoxin\nreference laboratory; and therefore, common biotoxin standards in agricultural products/foods have\nbeen collected, and we established a mycotoxins mass spectrometry database. Some of the standards are\ncommercially available, and some are donated by collaborative laboratories. For a few metabolites with\nan annotation level of 1, in addition to MS1 exact mass number (±5 ppm, m/z tolerance window) and\nMS/MS spectrum comparison, retention time comparison was also added. We present chromatographic\nretention time comparison plots of compounds in the aflatoxin biosynthetic pathway as an example to\ndemonstrate this. We can see that the retention times of these metabolites in the standards and actual\nsamples show almost no drift, which indicates that our chromatographic separation system is robust and\nhighly reproducible (Re-Figure 10). We also then selected two metabolites for peak alignment and\nfound that the peaks exhibited excellent reproducibility, demonstrating the robustness of our\nchromatographic method when analyzing different samples (Re-Figure 11). The robustness and\nreproducibility of peak intensity are crucial for quantitative results.\nRe_Figure 9| The mirror plot of MS/MS mass spectra by matching standards and samples MS/MS mass\nspectra. (a). Mirror plot of Aflatoxin B . (b). The mirror plot of 11-hydroxy O-methylsterigmatocystin\n1\n(HOMST). (c). The mirror plot of versicolorin B. (d). The mirror plot of versiconol. The black part\nabove represents the MS/MS mass spectra of the standards from our lab. The green section below\nrepresents the MS/MS mass spectra of the A.flavus samples.\nRe_Figure 10| The chromatograms of six representative metabolites in (a) standards and (b) A.flavus\nextracts. We can see that the retention times of these metabolites in the standards and actual samples\nshow almost no drift, except for Aflatoxin G (0.49s RT drift), which indicates that our chromatographic\n1\nseparation system is robust and highly reproducible.\nRe_Figure 11| The chromatograms' peak alignment of (a) aflatoxin B and (b) versicolorin B in Blank,\nand high reproducibility of our chromatographic separation system.\nWe also manually downloaded the mirror plots of the 355 metabolites from the GNPS website\ninto SVG format files one by one, and then uploaded them to the figshare platform for reviewers to\ncheck (The figshare links: https://figshare.com/articles/dataset/_b_554AFs_1E56E3_no_gap-\nfilling_Mirror_plot_b_/30434335). We hope this mirror plot dataset will increase reviewers' and readers'\nconfidence in metabolite identification. The 355 mirror plots of identified metabolites included the A.\nflavus secondary metabolites identified in the in-house library and the GNPS-library, such as aflatoxin\nB , aflatoxin B , aflatoxin G , aflatoxin G , O-methylsterigmatocystin, sterigmatocystin, cyclopiazonic\n1 2 1 2\nacid, aflatrem, Alternariol, Penicillin G, and Polanrazine b. The mirror plot contains MS/MS spectral\nmatches between standards and samples, and cosine similarity scores. We selected mirror plots for six\nA. flavus secondary metabolites as examples to demonstrate the reliability of our secondary metabolite\nannotation results from GNPS (Re-Figure 12 and Extended Data Fig. 5, Level 2). These 355\nmetabolites are unique metabolites, and the adduct peaks have been deduplicated (for example, when a\ncompound contains [M+H]+ and [M+Na]+, only the metabolic characteristics of [M+H]+ are retained).\nIn addition, we have compiled the GNPS job IDs generated by different MS2 noise thresholds and with\nor without gap-filling treatment optimization experiments into a table for reviewers to check (Re-Table\n1). MetDNA3 only outputs MS/MS similarity score results, and their database is not open source;\ntherefore, we cannot obtain mirror plots. In addition, we provide chromatographic retention times and\nMS/MS chromatograms of standards and A.flavus mycelial extracts of six precursor metabolites in the\naflatoxin biosynthetic pathway, including 5-methoxysterigmatocystin, Norsolorinic acid, Averantin,\nAverufin, Versicolorin B, and Versicolorin A for the first time (Re-Figure 13-18). Finally, we present\nMS/MS spectra information for Ditryptophenaline, cAATrp, and speradine A in our dataset and in\nDrott's4 and Valdet's5 articles (Re-Figure 19-21). We will submit the MS/MS spectra of these\nmetabolites to the GNPS library.\nRe_Figure 12| The mirror plot of MS/MS mass spectra by matching samples and GNPS library MS/MS\nmass spectra. (a). Mirror plot of Aflatoxin B . (b). The mirror plot of Aflatoxin B . (c). The mirror plot\n1 2\nof Aflatoxin G . (d). The mirror plot of Aflatoxin G . (e). The mirror plot of Sterigmatocystin. (f). The\n1 2\nmirror plot of Aflatrem. The black part above represents the MS/MS mass spectra of the samples of\nA.flavus extracts. The green section below represents the MS/MS mass spectra of the GNPS library.\nRe-Table 1. GNPS job analysis results at different MS2 noise thresholds and with or without gap-filling\ntreatment optimization experiments.\nGNPS analysis task name GNPS analysis results id links\n554AF.posi.1E5E0.mzmine4-20251015-gafilling.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=4a0b5cd205b5415187459bf0562cb379\n554AF.posi.1E5E0.mzmine4-20251015-nongafilling.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=c4efcb700fd94df2ba49bd82773089bc\n554AF.posi.1E5E2.mzmine4.gapfilling_FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=5f94976d716247e5b48499a38b23b74d\n554AF.posi.1E5E2.mzmine4.non.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=c44ec92e6d0e4fbeb5fd5046e06eca3e\n554AF.posi.1E5E3.mzmine4.non.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=cd0764bbf371468e9224ee776fa1e098\n554AF.posi.1E5E3.mzmine4.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=30df7feb3cc54f5ab1e0419d3fe9ebd1\n554AF.posi.1E52E3.mzmine4.20250925.nongapfilling.0.6.6-0 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=59610df36cd8426b95258c17b15fec55\n554AF.posi.1E52E3.mzmine4.20250925.gapfilling-FBMN.0.6.6 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=7277f7b15d654c0ebd6212fff1599da6\n554AF.posi.1E54E3.mzmine4.20250925.nongapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=ac76b79fafe743629e6f89c035285ef9\n554AF.posi.1E54E3.mzmine4.20250925.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=9d0c47ea61164865a2ea2abbbb69c6ee\n554AF.posi.1E56E3.mzmine4.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=9bcbcfc1ba7945a887c8b4a418265b97\n554AF.posi.1E56E3.mzmine4.non.gapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=1d02472cebf54e63999ee7c0d7891a5d\n554AF.posi.1E58E3.mzmine4.gapfilling.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=29fd7042f5da4220b0b740b1a1c28ea6\n554AF.posi.1E58E3.mzmine4.non.nongapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=5ac87ecae815443ab077517a9b8a49ce\n554AF.posi.1E5E4.mzmine4-20250925.gapfilling.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=9a89a5fa694647aa853137bdf8263a22\n554AF.posi.1E5E4.mzmine4-20250925.nongapfilling.FBMN.0.6.6_V2 https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=902298bdfc5d4ddd9c3b05fe1411f2ad\nRe_Figure 13| Chromatographs, retention time and MS/MS spectra of Norsolorinic acid (NOR) in\nsamples and standards (The negative ion mode has a higher MS1 response intensity).\nRe_Figure 14| Chromatographs, retention time and MS/MS spectra of Averantin in samples and\nstandards (Negative mode).\nRe_Figure 15| Chromatographs, retention time and MS/MS spectra of Averufin in samples and\nstandards (Negative mode).\nRe_Figure 16| Chromatographs, retention time and MS/MS spectra of Versicolorin B in samples and\nstandards (Negative mode).\nRe_Figure 17| Chromatographs, retention time and MS/MS spectra of Versicolorin A in samples and\nstandards (Negative mode).\nRe_Figure 18| Chromatographs, retention time and MS/MS spectra of 5-Methoxysterigmatocystin in\nsamples and standards.\nRe_Figure 19| Chromatographs, retention time and MS/MS spectra of Ditryptophenaline in our dataset\nand Drott’s articles4 (Positive mode).\nQC-1 #6866 RT:14.22 AV:1 NL:4.13E5\nF:FTMS + c ESI d Full ms2 271.0381@hcd40.00 [75.0000-282.0000]\n130.06438\n100\n90\n80\n70\ne\nc\nn 60\na\nd\nn\nu\nb A 50\ne\nv\nita\nle 40\nR\n30\n20\n270.27737\n10 116.10645 228.99994\n243.01569\n88.07515 102.09091 158.15312 198.08365 213.00520 257.03174\n0\n80 100 120 140 160 180 200 220 240 260 280\nm/z\nRe_Figure 20| The MS/MS spectra of cyclo-acetoacetyl-L-tryptophan (cAATrp) in our dataset (above)\nand Valdet’s article5 (below).\nP3-QC-11 #6529 RT:12.44 AV:1 NL:2.47E5\nF:FTMS + c ESI d Full ms2 367.1648@hcd40.00 [80.0000-378.0000]\n367.16473\n100\n160.07561\n90\n80\n349.15421\n70\ne\nc\nn 60\na\nd n 182.08107\nu 226.12244\nb A 50\ne\nv ita le 40 266.11725 307.14365\nR\n30\n190.08612 325.15417\n20\n214.12250 238.12300 279.14893\n10 98.06020 140.07033\n130.06516\n0\n80 100 120 140 160 180 200 220 240 260 280 300 320 340 360\nm/z\nRe_Figure 21| The MS/MS spectra of speradine A in our dataset (above) and Valdet’s article5 (below).\nMinor comments:\nL853. Change polarized metabolites to polar metabolites.\nResponse: We have corrected this now.\nL855. Define “proper amount”, is there any range of beads used to sample homogenization?\nResponse: Thanks, yes, we usually add 4~6 beads. The description of \"proper amount\" has been revised\nto “4~6 beads” in the manuscript.\nL863. Did the author use extraction blanks? The authors state that only methanol is used as a\nsolvent blank.\nResponse: Thank you very much for the reviewer's questions. During the extraction experiment, we set\nup a single tube containing two extraction solvents as an extraction blank control. This tube did not\ncontain any A. flavus mycelium sample, but instead contained 4-6 pieces of bread and two internal\nstandard solvents. We use this extraction blank control to check the background information. For\nexample, whether it is a substance inherent in the solvent or a substance introduced from bread. We have\nnow added this information to the manuscript in lines 939-941. The methanol is used as a solvent blank\nin running the LC-HRMS sequence.\nL870-873. Use names instead of molecular formulas to define water, methanol, formic acid and\nammonium formate.\nResponse: We sincerely thank the reviewer for the revision suggestions. We have replaced all the\nmolecular formula definitions of water, methanol, formic acid and ammonium formate with name\ndefinitions in the whole text.\nL873-877. Gradient must be addressed based on organic mobile phase (B) rather than water (A).\nThis section needs to be corrected.\nResponse: We would like to express our gratitude to the reviewer for their comments. As requested,\nwe have made corrections to the original article that reflect the request.\nL1077. Change porsitive for positive.\nResponse: We thank the reviewer for pointing out this error, which has been corrected.\nL1079. Massive data is not available.\nResponse: Thank you to the reviewer for pointing out this issue. Currently, the Massive data will be\nreleased after the manuscript acceptance. We have provided the reviewer with an access link and\npassword.\n[redacted]\nInclude any potential limitations related to mass spectrometry and metabolomics data analysis in\nthe discussion section. Also, the authors should ensure that mass spectrometry-related\nterminology is consistently and correctly formatted throughout the manuscript. For example, the\nmass-to-charge ratio should be written as m/z (in italics and with a lowercase \"m\") following\nstandard conventions. Please review the text to correct any inconsistencies in the formatting of\nthis and related terms.\nResponse: We would like to thank the reviewer for their valuable comments. In the Discussion section,\nas requested, we have added the potential limitations related to mass spectrometry and metabolomics\ndata analysis, as well as our future effort directions. Specifically: “Under the currently used medium\nculture conditions, metabolite structures of nine of the sixteen experimentally resolved biosynthetic\ngene clusters (BGCs) could be identified in the metabolome data. We note that a single culture condition\n(universal PDA medium here) is unlikely to be sufficient to describe the comprehensive specialized\nmetabolic diversity of fungi, as a significant proportion of BGCs are typically inactive and require\nspecific conditions to be activated, expressed, and biosynthesized17,18. Our pan-transcriptome analysis\nof 28 strains from different subpopulations under the same culture conditions validated that a large\nproportion of BGC genes were expressed at low levels or not at all (Supplementary Table 6).\nAlternatively, limitations in extraction solvents (the protocol used here does not include dedicated\nlipidome extraction solvents) and separation conditions (e.g., HILIC column for polar metabolites) may\nhave prevented the detection of additional metabolites. In terms of mass spectrometry methods, adopting\nboth DDA (Data-Dependent Acquisition) and DIA (Data-Independent Acquisition) data acquisition\nmodes could facilitate the identification of more metabolites. In metabolomics data analysis, newly\ndeveloped AI-driven methods such as MassCube19, MS2Query20, Flash Entropy Search21, and the\nSpectral Denoising Search22 may enhance both the number and accuracy of metabolite annotations,\nsurpassing traditional cosine similarity-based alignment algorithms. Currently, many fungal metabolites\nthat have been isolated and identified, however, still lack high-quality MS/MS spectral information in\npublic databases, which limits the number of specialized metabolites we can identify.”\nWe also appreciate the reviewers' suggestions and have consistently and correctly formatted the\nmass spectrometry-related terminology terms such as “m/z”, throughout the manuscript, as requested by\nthe reviewers.\nConclusion\nThe integration of omics layers in this study is conceptually strong, and the manuscript addresses\na highly relevant biological question. However, multiple issues in the mass spectrometry data\nprocessing pipeline and metabolite annotation framework require careful revision. Addressing\nthese concerns will significantly strengthen the manuscript’s technical foundation and improve its\nreproducibility and impact.\nResponse: We thank the reviewers for their expert review of our manuscript and for their positive\nfeedback and affirmation of its merits. We optimized and adjusted the parameter thresholds as requested\nby the reviewers and reanalyzed the untargeted metabolomics data. Multiple issues in the mass\nspectrometry data processing pipeline and metabolite annotation framework have been carefully revised\none by one. This greatly improved the quality of our manuscript and the clarity of our results.\n\nResponse to Reviewer #3\nReviewer #3 (Remarks to the Author):\n------------------------------------------------------------------------------------------------------------\nThe authors have sequenced the genomes of 544 strains of the fungus Aspergillus flavus. This large\ndataset was combined with previously published genomes, as well as metabolomics and\ntranscriptomics data. Using various computational approaches, new insights were obtained\nregarding evolution of biosynthetic gene clusters. Of most interest, the data was used to identify\nnovel genes putatively involved in various aspects of secondary metabolism, and this was\nconfirmed using gene knockout and metabolic phenotyping.\nResponse: We would sincerely like to thank the reviewers for their review of our manuscript, as well\nas for their positive feedback and affirmation of its quality.\nAnother way to view the high-resolution revised manuscript:\nDue to network limitations and files that are too large, we are unable to upload versions of\nthe rivised manuscript containing ultra-high-resolution figures to the NC revision systems.\nOur revised manuscript, which includes a revision history and a clean version have been\nalso uploaded to WUR Drive\n[redacted]\n---------------------------------------------------------------------------------------------------------------------------\nThe manuscript is generally well written, but aspects of the data presentation should be\nconsiderably improved. Specifically, there are a lot of figures that are unclear, either due to tiny\nfonts, mislabeling, incorrect ordering, unlabeled axes, or incomplete explanations for the average\nreader. See specific comments below.\nResponse: We sincerely thank the reviewers for their valuable suggestions regarding the aspects of the\ndata presentation in our manuscript. We have meticulously revised the unclear figures, which were either\ndue to tiny fonts, mislabeling, incorrect ordering, unlabeled axes, or incomplete explanations in the\nmanuscript, in accordance with the reviewers' requirements. Subsequently, we carefully addressed each\nspecific issue raised by the reviewers.\nThe pangenome was determined based on “567 representative genomes” (line 252). Unless I\nmissed it, it is not explained in the corresponding section of the material and methods how these\n567 genomes were representatively selected. This has implications for the use of the rarefaction\ncurve in Figure 3b, for example. If the 567 genomes comprised the 544 newly-sequenced Chinese\ngenomes plus some additionals, then this is hardly a representative set.\nResponse: We sincerely thank the reviewers for pointing out these concerns. The pan-genome\nconstructed from the genomes of the previous 567 strains did indeed suffer from insufficient sample\nrepresentation. Therefore, following the reviewers' suggestions, we reassembled and annotated 410 new\nstrains from previously published data using the same pipeline, and then constructed a new pan-genome\ncomprising 977 strains. The genomes of about 75 strains in 1051 strains were filtered out due to issues\nsuch as low assembly quality. These 977 strains cover 9 countries, including 75 clinical strains23 and\n902 environmental strains isolated from soil, peanuts, corn, cotton, and other sources4,24. We hope that\nthe current edition of the pangenome is representative of the available A. flavus genome globally. We\nadded a column for \"Pangenome samples\" to the strain metadata table in the supplementary materials\ntable 1 to label the samples used for pangenome construction. In fact, we constructed three versions of\nthe A.flavus pangenome, starting with 259 strains from China (Re-Figure 22a), increasing to 567 (Re-\nFigure 22b), and now reaching 977 strains (Re-Figure 22c). In the second edition of the pan-genome\nanalysis, we did indeed primarily analyze samples from China (n=539), and also selected 28 strains\ngenomes from different geographical environments from NCBI. This version of the pangenome does\nnot include clinical samples and strains from a wider range of environmental sources. We therefore re-\nassembled and annotated 410 strains from previously published data, then constructed a new pan-\ngenome comprising 977 strains. Following this pangenome construction, we performed rigorous\nbacterial sequence decontamination analysis by aligning the NR database. We did not perform a rigorous\ndecontamination step during the pangenome construction process of 567 strains version. Therefore, the\ntotal number of genes in the previously submitted pangenome was 16,416, while the actual number after\ndecontamination was 15,035. Compared to the pangenome constructed using 567 strains, the\npangenome constructed using 977 samples increased in size by 2.58 Mb and added 647 genes. We found\nthat as the genome size increases, the pangenome size and the number of predicted genes also increase,\nbut the increase is relatively small, unlike in prokaryotes. It is evident that although the number of\ngenome samples increased by 410, the number of genes increased only slightly. Orthology inference\nrevealed that the A. flavus pangenome is composed of 15,628 protein-coding genes, which is 3,181\ngenes more than the reference genome of A.flavus NRRL3357 (GCA_009017415.1)25. These non-\nredundant orthogroups could be subdivided into core (in all isolates, n = 9,584), accessory (in 5-95%\nof the isolates, n = 5,085) and ‘unique’ (in <5% of the isolates, n = 959) genomes (Fig. 3a) after rigorous\nremoval of bacterial and other prokaryotic sequences by aligning the NR database. Here we need to\ncorrect an analysis error in our previously submitted manuscript. Among them, the conserved core\ngenome accounts for 61.3%, the accessory genome accounts for 32.5%, and the unique genes account\nfor the remaining 6.2% of genes. The pangenome exhibits a closed trend as the population size increases\n(Fig. 3b). Further results regarding the pangenome were updated in the manuscript.\nRe_Figure 22| Pangenome simulation curves constructed with three different numbers of A. flavus\ngenomes. a. Pangenome simulation curves constructed with 259 A. flavus genomes. b. Pangenome\nsimulation curves constructed with 567 A. flavus genomes. c. Pangenome simulation curves constructed\nwith 977 A. flavus genomes.\nFigure 3d shows the enrichment of KEGG terms in the accessory gene families. Enrichment\nsuggests that there is over-representation in the accessory genes compared to something else (core,\nor all other genes?). The figure appears to just show gene counts, though. Not enrichment.\nResponse: Thanks to the reviewers for pointing out this issue. As the reviewers noted, we previously\nvisualized the number of genes annotated in the Kyoto Encyclopedia of Genes and Genomes (KEGG)\npathway rather than an enrichment analysis. Therefore, according to the reviewers' requirements, we\nperformed KEGG enrichment analysis on the core, accessory, and unique genes after updating the\npangenome. KEGG enrichment analysis confirmed that genes responsible for genetic information\nprocessing are dominant in the core genome, alongside genes related to primary metabolism and cellular\nprocesses. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis confirmed that\ngenes responsible for genetic information processing are dominant in the core genome, alongside genes\nrelated to primary metabolism and cellular processes in the top 20 pathway (Re_Figure 23c). 81% of\naccessory genes and 72% of unique genes were annotated in biosynthesis of primary and specialized\nmetabolites(Re_Figure 23d,e).\nRe_Figure 23| Super-pangenome analysis reveals A.flavus population conserved and unique genomic attributes. (a)\nPresence/absence matrix heatmap of 17,114 orthologous genes identified from 977 representative A. flavus isolates\ngenomes. The pangenome contains core (orthogroups present in all isolates), accessory (orthogroups present in 5–95%\nof the isolates) and unique (orthogroups present in less than 5% of the isolates) genomes. (b) Rarefaction curve of gene\nfamily variation in the pan-genome and core genome as the number of genomes increases. (c). KEGG enrichment analysis\nof the core gene set. (d) KEGG enrichment analysis of the accessory gene set. (e) KEGG enrichment analysis of the\nunique gene set. (f) Cluster analysis of the distribution of gain and loss of accessory genes recapitulates clade-specific\nevolutionary trends. (g) Clade-specific patterns emerge in the expression profiles of accessory genes in different W-\nclades.\nThe metabolomic analysis was performed in only 1 condition, so it could mean that some BGCs\nwere not active under the tested condition (and likely has an effect on the data presented in Figure\n2bcd). For such a large set of strains it is not required to test additional conditions, but it could be\ndiscussed somewhere.\nResponse: Thank you to the reviewer for this suggestion. As suggested by the reviewers, we have\nhighlighted this limitation in the Discussion section. Indeed, under the currently used medium culture\nconditions, metabolite structures of nine of the sixteen experimentally resolved biosynthetic gene\nclusters (BGCs) could be identified in the metabolome data. We note that a single culture condition\n(universal PDA medium here) is unlikely to be sufficient to describe the comprehensive specialized\nmetabolic diversity of fungi, as a significant proportion of BGCs are typically inactive and require\nspecific conditions to be activated, expressed, and biosynthesized17,18. Our pan-transcriptome analysis\nof 28 strains from different subpopulations under the same culture conditions validated that a large\nproportion of BGC genes were expressed at low levels or not at all (Supplementary Table 6).\nLines 331-332. Figure 7 is referenced before even figure 4. Probably the reference is wrong,\nbecause I don’t see how the authors draw these important conclusions from this figure (nor can I\nfind the actual figure that they likely intended to reference).\nResponse: We would like to thank the reviewers for pointing out this issue. This is likely an error that\narose after multiple rounds of revisions. The connections between the preceding and following sentences\nwere weak, indicating a logical flaw. As a result, we have adjusted the order of this content in lines 331-\n332 of the previous review version (In the current modified version, this corresponds to lines 344-350).\nThis statement is indeed not suitable for description in this section of results, as it needs to be compared\nwith subsequent results and should be done in the discussion section. Therefore, we deleted them here.\nFigure 4 is not easy to interpret. The axes in 4b and 4c are unlabeled, making it very unclear what\nis depicted here and what conclusions can be drawn. It is very unclear to me what is depicted in\n4ghi. Note that in line 467 the panels 4ghi are erroneously referred to as 4fgh, adding to the\nconfusion.\nResponse: Thank you to the reviewer for pointing out this shortcoming. After receiving comments from\nthe reviewer, we realized that Figure 4 did indeed have issues with unclear labeling and data\nvisualization confusion. Therefore, we made significant revisions based on the reviewers' comments.\nDue to limited figure space, displaying the BGC core gene variation and corresponding metabolite\nintensity change for 12 specialized metabolites would become blurred. Figure 4a has been simplified.\nWe use aflatoxin BGC as an example to demonstrate that BGC variation can only partially explain\nmetabolic differentiation. Considering reviewer #2's minimum requirement for merging two mass\nspectrometry datasets, we only retained aflatoxin B (cosine score=0.99) and cyclopiazonic acid (cosine\n1\nscore=0.96) with annotation level 2 and high peak intensity in both datasets to demonstrate that BGC\nvariation can only partially explain metabolic differentiation. The results are shown in Figure 4a\n(aflatoxin B ) and supplementary material Figure S8 (cyclopiazonic acid). The labels were added to\na more commonly used visualization (Re_Figure 24). Previously, we used fold change values, which\nmay have been somewhat difficult to understand. Furthermore, we removed Figure 4ghi because the\ngene KEGG enrichment analysis in Figure 4d already contained information on differential expression\nin different A.flavus populations across the three metabolic pathways that Figure 4ghi (Gene set\nenrichment analysis (GSEA) of aflatoxin, tryptophan, and fatty acid biosynthesis pathway) was intended\nto demonstrate. Previously, we used Gene Set Enrichment Analysis (GSEA) simply to use another\npathway enrichment method to further validate the reliability of the pathways enriched in Figure 4d.\nWith these revisions, we hope to have resolved the issues that confused the reviewers. Alternatively, we\nadded the results of pathway enrichment analysis of differentially expressed genes in C-clade6/clade3\n(Re_Figure 24g) and C-clade6/clade4 (Re_Figure 24f) to emphasize that differences in transcriptional\nlevels among different subpopulations may be an important driver in metabolic differentiation.\nRe_Figure 24| (a) Experimentally characterized BGC core gene presence/absence variation (PAV) only partially\nexplains specialized metabolites' abundance profiles. For example, in aflatoxin BGC, 34% of W-clade 1 and 14% of W-\nclade 6 strains showed multiple peripheral gene loss events, 44% strains without core gene loss events still showed no/low\naflatoxin production. (b) The differentially expressed gene ranking dotplot of 58 BGCs of 28 strains in Chinese\nclade4/clade1. (c) Results of pathway enrichment analysis of differentially expressed genes in C-clade4/clade1. (d)\nResults of pathway enrichment analysis of differentially expressed genes in C-clade4/clade3. (e) Results of pathway\nenrichment analysis of differentially expressed genes in C-clade3/clade1. (f) Results of pathway enrichment analysis of\ndifferentially expressed genes in C-clade6/clade3. (g) Results of pathway enrichment analysis of differentially expressed\ngenes in C-clade6/clade4.\nFigure 5. The labels for 5e and 5d are swapped…\nResponse: We thank the reviewers for pointing out this error. It has been corrected in the manuscript\nsubmission.\nLine 431. “dispatched” is an odd word choice. Maybe just “used”?\nResponse: Thank you to the reviewer for pointing out this issue. Based on the reviewer’s suggestion,\nwe have replaced it with \"used\".\nLines 196-198. “We observed that the genomes of these A. favus strains (W-clade 8) harboured\nhomologues of 156 virulence genes….” What about the other strains? And were these homologues\nof virulence genes over-represented in the clinical strains versus the environmental strains?\nResponse: We thank the reviewers for their interest in this distribution. Among the 156 reported\nvirulence genes in Aspergillus fumigatus26, we identified 149 virulence genes present in the Aspergillus\nflavus population through large-scale homology clustering using OrthoFinder27. After obtaining the\ndistribution matrix of virulence gene copy number distribution data matrix, we performed cluster\nanalysis and multiple group difference tests. The differential distribution of virulence gene copy number\namong strains was determined by the Wilcox test between groups (P-value = 0.05). Overall, we found\nthe vast majority of the virulence genes (n=137, 137/149=92%) in both clinical and environmental\nstrains, and the phenomenon of homologues of virulence genes over-represented in clinical strains\nversus environmental strains was not widespread(n=12, 12/149=8%). (Re_Figure 25a, and\nSupplementary Material Figure S1 in the revised manuscript). It is true that a few virulence genes\nhave higher copy numbers in clinical strains than in environmental strains, for example, the homologues\nvirulence gene of Afu5g03790(Iron transport multicopper oxidase fetC in Aspergillus fumigatus) (Fig.\n1a,d). In addtion, the homologue virulence gene of Afu8g00230(Chain A, Verruculogen synthase in\nAspergillus fumigatus Af293) is mainly distributed in W-Clade 2, mainly including environmental non-\naflatoxigenic clade strains (Fig. 1a in the revised manuscript). The distribution of the top 12 virulence\ngenes that differed in different subpopulations is shown in Re_Figure 25b.\nRe_Figure 25| The heatmap shows the copy number variation distribution of (a) 149 A. fumigatus\nvirulence homologs and (b) the top 12 differential distribution of virulence homologs in 1003 A. flavus\ngenomes.\nIn lines 142-177 (and elsewhere, including figures). There is inconsistency in the naming of the\npopulations: population B, PopB, pop B\nResponse: We thank the reviewer for pointing out this issue. The entire manuscript was checked and is\nnow using the format Pop A, Pop B, Pop C, Pop D consistently.\nSome of the text in the figures (e.g. Figure 2) is way too small, even in the original files.\nEspecially in Figure 2b.\nResponse: We would like to thank the reviewers for highlighting this issue. We have revised Fig. 2b in\nline with reviewers' comments. The main revisions include the following: i) due to limitations in the\navailable plot space, the number of metabolites in the different metabolic patterns has been reduced to\nmore clearly demonstrate the metabolic differentiation patterns of the A. flavus population; ii) in addition,\ngiven the large number of genomes and metabolites, and the limited space in the main figure of the\nmanuscript, we followed the reviewers' comments and tried our best to clearly visualize only the most\nimportant information in the main figure of the article; and iii) the phylogenetic tree in Figure 2b has\nbeen reshaped to more clearly show the sample names (n>500) (Re_Figure 26). We hope they meet the\nreviewers' requirements.\nRe_Figure 26| (a) Phylogenetic tree, field sampling sites and phylogeographic patterns of A. flavus strains from China,\nand geographical distribution pattern of aflatoxin-producing capacity. 97% nodes in the phylogenetic tree have a\nbootstrap value greater than 0.98. (b) The phylogenetic tree from the population of A. flavus in China. The heatmap\ncontains the five A. flavus metabolic patterns mapping on the phylogenetic tree. Low aflatoxin-producing capacity clades\noften produce other metabolites instead. (c) -cyclopiazonic acid (-CPA) molecular family mainly found in clades 1\nand 3 (see Online Methods for metabolite annotation details). (d) The highest abundance of the aflatoxin molecular family\nwas observed in Chinese clades 2, 4, 5, and 6, and the lowest abundance was observed in Chinese clades 1 and 3.\nLine 313. Typos “In addition, by”\nResponse: Thank you. We have now corrected this in the manuscript (Line 330 in the currently revised\nversion).\nLine 488 Typo “functional”\nResponse: Thank you. We have now corrected this in the manuscript (Line 409 in the currently revised\nversion).\nLine 584. Typo “and”\nResponse: Thank you. We have now corrected this in the manuscript (Line 637 in the currently revised\nversion).\nLine 595. Typo “environmental”\nResponse: Thank you. We have now corrected this in the manuscript (Line 649 in the currently revised\nversion).\nReference\n1. Dunn, W.B. et al. Procedures for large-scale metabolic profiling of serum and plasma using\ngas chromatography and liquid chromatography coupled to mass spectrometry. Nature\nProtocols 6, 1060–1083 (2011).\n2. Barbier Saint Hilaire, P. et al. Evaluation of the High-Field Orbitrap Fusion for Compound\nAnnotation in Metabolomics. Analytical Chemistry 90, 3030–3035 (2018).\n3. Xie, H.L. et al. Monitoring Metabolite Production of Aflatoxin Biosynthesis by Orbitrap Fusion\nMass Spectrometry and a D-Optimal Mixture Design Method. Analytical Chemistry 90,\n14331–14338 (2018).\n4. Drott, M.T. et al. Microevolution in the pansecondary metabolome of Aspergillus flavus and\nits potential macroevolutionary implications for filamentous fungi. Proceedings of the\nNational Academy of Sciences 118, e2021683118 (2021).\n5. Uka, V. et al. Unravelling the Diversity of the Cyclopiazonic Acid Family of Mycotoxins in\nAspergillus flavus by UHPLC Triple-TOF HRMS. Toxins (Basel) 9(2017).\n6. Gil-de-la-Fuente, A. et al. Aspergillus Metabolome Database for Mass Spectrometry\nMetabolomics. Journal of Fungi 7(2021).\n7. van Santen, J.A. et al. The Natural Products Atlas 2.0: a database of microbially-derived natural\nproducts. Nucleic Acids Res 50, D1317–D1323 (2022).\n8. Dührkop, K. et al. SIRIUS 4: a rapid tool for turning tandem mass spectra into metabolite\nstructure information. Nature Methods 16, 299–302 (2019).\n9. Stravs, M.A., Dührkop, K., Böcker, S. & Zamboni, N. MSNovelist: de novo structure generation\nfrom mass spectra. Nature Methods 19, 865–870 (2022).\n10. Cao, L. et al. MolDiscovery: learning mass spectrometry fragmentation of small molecules. Nat\nCommun 12, 3718 (2021).\n11. Mohimani, H. et al. Dereplication of microbial metabolites through database search of mass\nspectra. Nat Commun 9, 4035 (2018).\n12. Schmid, R. et al. Integrative analysis of multimodal mass spectrometry data in MZmine 3.\nNature Biotechnology 41, 447–449 (2023).\n13. Shen, X. et al. TidyMass an object-oriented reproducible analysis framework for LC–MS data.\nNature Communications 13(2022).\n14. Uka, V. et al. Chemical repertoire and biosynthetic machinery of the Aspergillus flavus\nsecondary metabolome: A review. Comprehensive Reviews in Food Science and Food Safety\n19, 2797–2842 (2020).\n15. Pang, Z. et al. MetaboAnalyst 5.0: narrowing the gap between raw spectra and functional\ninsights. Nucleic Acids Res 49, W388–W396 (2021).\n16. Broadhurst, D. et al. Guidelines and considerations for the use of system suitability and quality\ncontrol samples in mass spectrometry assays applied in untargeted clinical metabolomic\nstudies. Metabolomics 14(2018).\n17. Brakhage, A.A. Regulation of fungal secondary metabolism. Nature Reviews Microbiology 11,\n21–32 (2012).\n18. Keller, N.P. Fungal secondary metabolism: regulation, function and drug discovery. Nature\nReviews Microbiology 17, 167–180 (2019).\n19. Yu, H. et al. MassCube improves accuracy for metabolomics data processing from raw files to\nphenotype classifiers. Nature Communications 16(2025).\n20. de Jonge, N.F. et al. MS2Query: reliable and scalable MS2 mass spectra-based analogue\nsearch. Nature Communications 14(2023).\n21. Li, Y. & Fiehn, O. Flash entropy search to query all mass spectral libraries in real time. Nature\nMethods 20, 1475–1478 (2023).\n22. Kong, F. et al. Denoising Search doubles the number of metabolite and exposome annotations\nin human plasma using an Orbitrap Astral mass spectrometer. Nature Methods 22, 1008–\n1016 (2025).\n23. Hatmaker, E.A. et al. Population structure in a fungal human pathogen is potentially linked to\npathogenicity. Nature Communications 16(2025).\n24. Gangurde, S.S. et al. Aspergillus flavus pangenome (AflaPan) uncovers novel aflatoxin and\nsecondary metabolite associated gene clusters. BMC Plant Biology 24, 354 (2024).\n25. Drott, M.T. et al. The Frequency of Sex: Population Genomics Reveals Differences in\nRecombination and Population Structure of the Aflatoxin-Producing Fungus Aspergillus\nflavus. mBio 11(2020).\n26. Rhodes, J. et al. Population genomics confirms acquisition of drug-resistant Aspergillus\nfumigatus infection by humans from the environment. Nature Microbiology 7, 663–674 (2022).\n27. Emms, D.M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative\ngenomics. Genome Biology 20, 238 (2019)."}, {"round": 2, "reviewer_comments": {"reviewer_2": "This letter is to inform you that after carefully reviewing the current and corrected version of the manuscript entitled \"Large-Scale Multi-Omics Profiling Reveals Environmental and Evolutionary Drivers in Phylogeographic Diversity of Fungal Specialized Metabolism\", I’m very glad to express that the authors have addressed my major concerns (metabolomics component of the study) in a thorough, transparent, and technically rigorous manner. This new version of the manuscript substantially improves the reliability and reproducibility of metabolomics analysis, which was the central focus of my earlier critique.\n\nBelow, I summarize my assessment:\n1. The authors corrected the overly permissive MS2 noise threshold (originally 1E2) and performed a systematic optimization across a wide range of thresholds. The selection of 6E3 is well justified based on comparative feature counts, metabolite overlap, and annotation stability. There is no more concern about this as presented in the current version.\n2. The authors performed a detailed comparison with and without gap filling, demonstrating: the impact on MS/MS feature counts, the risk of losing uncommon specialized metabolites, the effect on downstream quantification, and how gap-filled versus manually quantified data were separated in analysis. They now clearly describe how gap-filled features are treated and filtered, and this response is acceptable and demonstrates a solid understanding of the tradeoffs.\n3. Regarding the dataset fusion with MS1-only used from Drott et al. The authors state now using MS/MS files and even uploading 24 qualitative files to the GNPS repository. The integration is now performed according to minimal acceptable standards for cross-dataset untargeted metabolomics analysis. This adequately resolves my concern.\n4. In case of the metabolite identification confidence levels (Level 1, Level 2), the authors provided the corresponding GNPS job links, as well as retention time comparisons with authentic standards, cosine similarity metrics, chromatographic evidence, and the identification framework is now transparent and appropriately supported. This is a major improvement and fully satisfies my request.\n5. The authors now acknowledge important limitations, including single culture conditions, inactive BGCs, and incomplete fungal MS/MS coverage in public databases, among others. These additions provide a more realistic frame of the metabolomics results.", "reviewer_3": "The authors have sufficiently (and in high detail) addressed the comments that I raised."}, "author_response": "Response to Reviewer #2\nReviewer #2 (Remarks to the Author):\nThis letter is to inform you that after carefully reviewing the current and corrected version of the\nmanuscript entitled \"Large-Scale Multi-Omics Profiling Reveals Environmental and Evolutionary\nDrivers in Phylogeographic Diversity of Fungal Specialized Metabolism\", I’m very glad to\nexpress that the authors have addressed my major concerns (metabolomics component of the study)\nin a thorough, transparent, and technically rigorous manner. This new version of the manuscript\nsubstantially improves the reliability and reproducibility of metabolomics analysis, which was the\ncentral focus of my earlier critique.\nBelow, I summarize my assessment:\n1. The authors corrected the overly permissive MS2 noise threshold (originally 1E2) and\nperformed a systematic optimization across a wide range of thresholds. The selection of 6E3 is\nwell justified based on comparative feature counts, metabolite overlap, and annotation stability.\nThere is no more concern about this as presented in the current version.\n2. The authors performed a detailed comparison with and without gap filling, demonstrating:\nthe impact on MS/MS feature counts, the risk of losing uncommon specialized metabolites, the effect\non downstream quantification, and how gap-filled versus manually quantified data were separated\nin analysis. They now clearly describe how gap-filled features are treated and filtered, and this\nresponse is acceptable and demonstrates a solid understanding of the tradeoffs.\n3. Regarding the dataset fusion with MS1-only used from Drott et al. The authors state now using\nMS/MS files and even uploading 24 qualitative files to the GNPS repository. The integration is\nnow performed according to minimal acceptable standards for cross-dataset untargeted\nmetabolomics analysis. This adequately resolves my concern.\n4. In case of the metabolite identification confidence levels (Level 1, Level 2), the authors\nprovided the corresponding GNPS job links, as well as retention time comparisons with authentic\nstandards, cosine similarity metrics, chromatographic evidence, and the identification framework\nis now transparent and appropriately supported. This is a major improvement and fully satisfies\nmy request.\n5. The authors now acknowledge important limitations, including single culture conditions,\ninactive BGCs, and incomplete fungal MS/MS coverage in public databases, among others. These\nadditions provide a more realistic frame of the metabolomics results.\nResponse: We would like to express our sincere gratitude to the reviewer#2 for their encouraging\nfeedback. Your comments and suggestions significantly improved the quality of our article and\nmade the Mass Spectrometry analysis results more robust.\n\nResponse to Reviewer #3\nReviewer #3 (Remarks to the Author):\nReviewer #3 (Remarks to the Author):\nThe authors have sufficiently (and in high detail) addressed the comments that I raised.\nResponse: We are very grateful to Reviewer#3 for their approval of our revised manuscript."}]

Decision letter