Genetic Responses and Aflatoxin Inhibition during Co-Culture of Aflatoxigenic and Non-Aflatoxigenic Aspergillus flavus

Aflatoxin is a carcinogenic mycotoxin produced by Aspergillus flavus. Non-aflatoxigenic (Non-tox) A. flavus isolates are deployed in corn fields as biocontrol because they substantially reduce aflatoxin contamination via direct replacement and additionally via direct contact or touch with toxigenic (Tox) isolates and secretion of inhibitory/degradative chemicals. To understand touch inhibition, HPLC analysis and RNA sequencing examined aflatoxin production and gene expression of Non-tox isolate 17 and Tox isolate 53 mono-cultures and during their interaction in co-culture. Aflatoxin production was reduced by 99.7% in 72 h co-cultures. Fewer than expected unique reads were assigned to Tox 53 during co-culture, indicating its growth and/or gene expression was inhibited in response to Non-tox 17. Predicted secreted proteins and genes involved in oxidation/reduction were enriched in Non-tox 17 and co-cultures compared to Tox 53. Five secondary metabolite (SM) gene clusters and kojic acid synthesis genes were upregulated in Non-tox 17 compared to Tox 53 and a few were further upregulated in co-cultures in response to touch. These results suggest Non-tox strains can inhibit growth and aflatoxin gene cluster expression in Tox strains through touch. Additionally, upregulation of other SM genes and redox genes during the biocontrol interaction demonstrates a potential role of inhibitory SMs and antioxidants as additional biocontrol mechanisms and deserves further exploration to improve biocontrol formulations.

Non-tox 17, Tox 53 and their co-cultures produced different quantities of aflatoxin B 1 after growing in liquid medium for different time points (30,72 and 96 h) as indicated by significant interactions (F 4,29 = 207, p-value < 0.0001). Tox 53 started producing significant quantities of aflatoxin at 72 h of growth (Table 1). Very limited aflatoxin (<2 ppb) was detected in the biocontrol interaction samples consisting of Tox 53 and Non-tox 17 co-cultures, suggesting the presence of Non-tox 17 severely limited aflatoxin production by Tox 53. Additionally, aflatoxin degradation by Non-tox 17 may have resulted in lower aflatoxin [41], despite the addition of citrate buffer to limit aflatoxin degradation [39,40,43]. Non-tox 17 alone did not produce aflatoxin, thereby confirming its non-aflatoxigenic phenotype.

Fungal Biomass and Total RNA
Tox 53, Non-tox 17 and their co-cultures produced different amounts of mycelial biomass at 30 and 72 h (F 2,21 = 58.0, p-value < 0.0001). For each mono-and co-culture, there was more mycelial mass after 72 h (Figure 1). At both 30 and 72 h culture ages, Tox 53 produced less mycelia than Non-tox 17 and the co-cultures. Very little Tox 53 tissue was harvested at 30 h and the least squares estimate was not different from 0 (t 21 = 0.38, p-value = 0.71). In contrast to the amount of mycelial tissue harvested, the differences between Non-tox 17, Tox 53 and their co-cultures in amount of total RNA extracted did not vary between 30 and 72 h (F 2,18 = 1.82, p-value = 0.19) and culture age alone did not affect total RNA (F 1,18 = 2.54, p-value = 0.13). This suggested there were no differences in RNA extraction efficiencies or in RNA production within older tissue. Isolate type affected the total RNA (F 2,18 = 32.64, p-value < 0.0001) as less total RNA was extracted from Tox 53 than Non-Tox 17 and co-cultures.

RNA Sequencing
Between 11 and 30 million paired-end reads were sequenced for Non-tox 17, Tox 53 and co-cultures (Table 2). On average 37% of sequence reads mapped to genes with single-nucleotide differences between Non-tox 17 and Tox 53. In co-culture, only 3% of the unique reads were assigned to Tox 53, which would indicate very little presence or gene expression by Tox 53. That 3% was slightly more than 1% of reads misaligned to Tox 53 in Non-tox 17 pure cultures, suggested as little as 2% of the reads uniquely aligned to Tox 53 in co-culture. oxins 2021, 13, x FOR PEER REVIEW RNA (F2,18 = 32.64, p-value < 0.0001) as less total RNA was extracted from Tox Tox 17 and co-cultures.

RNA Sequencing
Between 11 and 30 million paired-end reads were sequenced for Nonand co-cultures (Table 2). On average 37% of sequence reads mapped to gene nucleotide differences between Non-tox 17 and Tox 53. In co-culture, on unique reads were assigned to Tox 53, which would indicate very little pre expression by Tox 53. That 3% was slightly more than 1% of reads misalig in Non-tox 17 pure cultures, suggested as little as 2% of the reads uniquely a 53 in co-culture.

Equating Growth with RNA Production
The observed proportion of reads that uniquely aligned to Tox 53 in co-culture was very low (approximately 0.03 or 3%). Therefore, the observed proportion of Tox 53 reads in co-culture was compared to the expected proportions of Tox 53 based on differences in growth (biomass) or RNA production between Tox 53 and Non-tox 17. The expected proportion of Tox 53 biomass in co-culture was based on mono-cultures and calculated as: p 53 biomass = (Tox 53 biomass (mg)) ÷ total biomass (Tox 53 biomass (mg) + Non-tox 17 biomass (mg)) Total biomass was an estimate of co-cultures' total biomass that assumed Tox 53 or Non-tox 17 do not influence the growth of either isolate. Expected proportion of Tox 53 was also calculated using RNA (µg/mg mycelium). Multicategorical data analysis (i.e., multiple contingency tables) compared the observed proportion of reads that uniquely aligned to Tox 53 in co-culture to the expected proportion based on Tox 53 biomass or RNA ( Figure 2). There was a significant interaction between the proportion of Tox 53 as determined by reads, biomass, total RNA and 30 or 72 h time points (F 4,22 = 9288, p-value < 0.0001). At 30 h, 3% of the reads were not significantly less than would be expected based on the growth of Tox 53, but at 72 h there were significantly fewer reads aligned to Tox 53 than would be expected based on both biomass and RNA production ( Figure 2). This indicated that co-culturing Tox 53 with Non-tox 17 decreased both RNA transcription and growth of Tox 53. Growth medium was buffered with citrate to maintain pH~4 [39,40,43] and avoid acidification from fungal growth which reduces aflatoxin production and fungal growth, suggesting the reduced Tox 53 growth and transcription during co-culture is unlikely solely from acidification by Non-tox 17. multiple contingency tables) compared the observed proportion of aligned to Tox 53 in co-culture to the expected proportion based o RNA ( Figure 2). There was a significant interaction between the pro determined by reads, biomass, total RNA and 30 or 72 h time points < 0.0001). At 30 h, 3% of the reads were not significantly less than wou on the growth of Tox 53, but at 72 h there were significantly fewer rea than would be expected based on both biomass and RNA productio dicated that co-culturing Tox 53 with Non-tox 17 decreased both RN growth of Tox 53. Growth medium was buffered with citrate to main and avoid acidification from fungal growth which reduces aflatoxin gal growth, suggesting the reduced Tox 53 growth and transcription unlikely solely from acidification by Non-tox 17.

Differential Gene Expression
Based on principle component analysis, the biological reps cluste however, 72 h co-cultures had the most variation (Figure 3a). Express between 30 and 72 h cultures and for each time point the co-cultures c Non-tox 17. The Non-tox 17 mono-cultures and co-cultures expresse Figure 2. Proportion of RNA sequence reads uniquely aligned to A. flavus Tox 53 and Non-tox 17 in co-culture vs. the expected proportions based on biomass and RNA production of isolates grown apart. Proportions were compared using a generalized linear mixed model with the logit link and binomial distribution. Different letters denote significance as per least squares means (odds) comparisons (α < 0.05).

Differential Gene Expression
Based on principle component analysis, the biological reps clustered closely together; however, 72 h co-cultures had the most variation ( Figure 3a). Expression patterns differed between 30 and 72 h cultures and for each time point the co-cultures clustered closely with Non-tox 17. The Non-tox 17 mono-cultures and co-cultures expressed between 1500 and 2000 genes more than Tox 53 with similar amounts of gene downregulation in Tox 53 ( Figure 3b). However, very few genes differed in expression between the co-cultures and Non-tox 17 mono-cultures because Non-tox 17 growth dominated Tox 53 growth in co-cultures.

Functional Analysis of Differentially Expressed Genes
Like overall gene expression, gene ontology functional categories differed between Tox 53 and both Non-tox 17 and co-cultures, and very few categories differed between cocultures and Non-tox 17, likely due to Non-tox 17′s growth dominating the co-cultures ( Table 3). Hundreds of genes involved in oxidation/reduction reactions and encoding proteins localized to apoplastic or extracellular spaces were more differentially expressed than any of the other categories between Tox 53 and Non-tox 17. Genes that were typically upregulated in Non-tox 17 were also downregulated in Tox 53. Those trends were mostly consistent between co-cultures and Tox 53 as well. More than 50 zinc finger transcription factor genes were upregulated in Non-tox 17 and co-cultures compared to Tox 53. Most other gene categories were only differentially regulated at 30 or 72 h. Most of the secondary metabolite genes involved in aflatoxin, sterigmatocystin and cyclopiazonic acid production were all downregulated in Non-tox 17 and co-cultures compared to Tox 53. Nontox 17 does not have any of the genes from these mycotoxin biosynthesis pathways. A few genes from these pathways were upregulated in co-cultures indicating there was some Tox 53 growth present; however, most pathway genes were not expressed at a detectable level. Genes from the biosynthesis pathways of the putative asperfuranone and characterized imizoquin secondary metabolites were also upregulated in Non-tox 17 compared to Tox 53.

Functional Analysis of Differentially Expressed Genes
Like overall gene expression, gene ontology functional categories differed between Tox 53 and both Non-tox 17 and co-cultures, and very few categories differed between co-cultures and Non-tox 17, likely due to Non-tox 17 s growth dominating the co-cultures ( Table 3). Hundreds of genes involved in oxidation/reduction reactions and encoding proteins localized to apoplastic or extracellular spaces were more differentially expressed than any of the other categories between Tox 53 and Non-tox 17. Genes that were typically upregulated in Non-tox 17 were also downregulated in Tox 53. Those trends were mostly consistent between co-cultures and Tox 53 as well. More than 50 zinc finger transcription factor genes were upregulated in Non-tox 17 and co-cultures compared to Tox 53. Most other gene categories were only differentially regulated at 30 or 72 h. Most of the secondary metabolite genes involved in aflatoxin, sterigmatocystin and cyclopiazonic acid production were all downregulated in Non-tox 17 and co-cultures compared to Tox 53. Non-tox 17 does not have any of the genes from these mycotoxin biosynthesis pathways. A few genes from these pathways were upregulated in co-cultures indicating there was some Tox 53 growth present; however, most pathway genes were not expressed at a detectable level. Genes from the biosynthesis pathways of the putative asperfuranone and characterized imizoquin secondary metabolites were also upregulated in Non-tox 17 compared to Tox 53.

Gene Expression in the Aflatoxin Biosynthesis Cluster
The Non-tox 17 isolate does not have aflatoxin cluster genes, explaining the low expression levels indicated in Table 4. Some genes were expressed at a higher level in co-culture compared to Non-tox 17 alone, indicating limited growth of Tox 53 in co-culture. However, since there were fewer than 10 reads per gene at 72 h, there was very little expression of aflatoxin cluster genes, which is supported by the lack of aflatoxin production in co-culture. More genes were differentially expressed, and greater fold differences occurred at 72 h suggesting the lack of detectable aflatoxin at 30 h was due to very low expression of some aflatoxin cluster genes.

Genes Highly Upregulated in Biocontrol Non-Tox 17 Compared to Tox 53
To understand what specific genes may contribute to Non-tox 17 s ability to outcompete Tox 53 and reduce its aflatoxin production during the biocontrol interaction, genes with the greatest upregulation (Log 2 -fold change ≥ 8) in Non-tox 17 compared to Tox 53  (Table 5a). Differential gene expression between Non-tox 17 and Tox 53 was similar to that of co-culture and Tox 53 alone, likely due to limited growth of Tox 53 in co-culture. The upregulated genes in Non-tox 17 compared to Tox 53 represent a diversity of potential functions including oxidation/reduction reactions, peroxisome production, metabolism, and protein-protein interactions. These functions are consistent with the predominant gene functions identified by functional enrichment analysis. Interesting, most of the highly expressed genes were on chromosome 5. AFLA_060320 and AFLA_060350, and AFLA_095290 and AFLA_095300 were co-located, but there was no similar trend for surrounding genes to be differentially expressed in those regions of the chromosome (Table S1). However, nearby genes AFLA_062960, 062980, and 062990 were upregulated and are in a secondary metabolite cluster 20 as predicted by SMURF [44,45]. Most of the remaining genes in cluster 20 were also upregulated in both Non-tox 17 and co-cultures, suggesting a potential role for differential secondary metabolism during the biocontrol interaction (Table 5b). Gene cluster 20 is predicted to produce asperfuranone [45,46] which was only found to be enriched at 72 h. Here, the genes were all differentially expressed at 30 and 72 h. When Tox 53 grew alone, no reads aligned to genes in the latter portion of SMURF cluster 20, starting with AFLA_062940. Since only AFLA_062800-AFLA_062880 are required for asperfuranone production [46], the remaining genes in cluster 20 are likely a separate secondary metabolite cluster only produced by Non-tox 17.  Values are color-scaled, blue is less than zero and red is greater than zero. A darker shade indicates a more negative or positive value and scaled to the maximum and minimum values in table. 2 Gene ID is the AFLA_# gene number designated in the A. flavus NRRL 3357 reference genome. 3 Chr is the chromosome where aflatoxin cluster genes are located. Values are color-scaled with a darker red tint indicating a more positive value and scaled to the maximum value in table. 2 Gene ID is the AFLA_# gene number designated in the A. flavus NRRL 3357 reference genome. Gene ID fonts color-coded in blue, red or gray are in close proximity on chromosomes. 3 Chr is the chromosome where each gene is located. 4 SM is the secondary metabolite cluster where a gene is located. The 1st number is predicted by antiSMASH, and the 2nd number is predicted by SMURF. Dashes indicate gene is not located in a predicted secondary metabolite cluster.

Genes Upregulated in Response to Tox 53
To further identify genes that may contribute to Non-tox 17 s ability to limit aflatoxin production during the biocontrol interaction, genes upregulated in response to Tox 53 were identified. Since co-cultures were dominated by Non-tox 17, genes that were further upregulated in Non-tox 17 in response to the presence of Tox 53 should have greater than 2 log 2 -fold differential expression in co-cultures than both Tox 53 and Non-tox 17 alone (Table S2). Only 10 genes fit these criteria and the fold changes only ranged from 2-3.2 when comparing co-cultures to Non-tox 17. If the criteria were loosened to greater than 1, 12 genes were upregulated compared to Non-tox 17 and Tox 53 in co-culture. Of the 12 genes, only AFLA_016350 (predicted NAD (P)H-dependent reductase) was expressed at a higher level in co-culture compared to Non-tox 17 at 30 h.
A closer inspection of the fold change values between co-cultures and Tox 53 revealed several genes that had slightly higher fold changes than Non-tox 17 alone vs. Tox 53, suggesting there could be a higher number of reads from Non-tox 17 in co-culture compared to Non-tox 17 alone, despite some relative expression levels indicating a lack of significance between co-culture and Non-tox 17. Reads per kilobase per millions of reads mapped (RPKM) values are shown for genes with higher RPKM values in co-culture than both Non-tox 17 and Tox 53 mono-cultures (Table 6a). Additionally, generalized linear models and post hoc least squares means (log odds) comparisons separated treatments based on normalized read counts per gene per total reads (proportion of total reads) ( Table 6), like DESeq2 differential gene expression methodology on read counts without using their data smoothing algorithms [47]. Without smoothing an additional 17 genes, that were expressed at slightly higher levels in co-culture vs. Tox 53 comparisons than in the individual Non-tox 17 vs. Tox 53 comparisons, were found to have significantly more reads mapped to the co-culture than both Tox 53 and Non-tox 17 alone, suggesting that co-culturing the two isolates induced expression of several genes in Non-tox 17. Like differential expression using the fold changes from DESeq2, the greatest differential expression was at 72 h and most genes were expressed in larger abundance in co-cultures based on RPKM at 72 h.
±5.6 *** 8.5, 46 FAD-dependent oxygenase   1 Reads per kilobase per millions reads mapped (RPKM) were calculated for each 30 and 72 h Non-tox 17, Tox 53 and co-cultures reported as means and standard error. RPKM are measures of read abundance, which account for the total reads in a single replicate and individual gene length. This allows for better comparisons between treatments and genes rather than directly reporting reads. RPKM values for a single gene within 30 or 72 h with different numbers of asterisks *, ** or *** are statistically different based on α = 0.05. Values are color-scaled red for each gene and darker shades indicate more positive values scaled to the maximum within an individual gene; 2 Gene ID is the AFLA_# gene number (i.e., 029700 = AFLA_029700) designated in the A. flavus NRRL 3357 reference genome. Gene ID fonts colorcoded in orange, purple, blue, red or gray are in close proximity on chromosomes. The red gene numbers are in secondary metabolite clusters. Orange genes were not predicted by a secondary metabolite program but are involved in kojic acid synthesis; 3 Chr is the chromosome where each gene is located; 4 SM is the secondary metabolite cluster where a gene is located. The 1st number is the SM cluster predicted by antiSMASH, the 2nd number is the SM cluster predicted by SMURF.
Many genes induced in response to Tox 53 functions are membrane transport proteins, metal-or heme-binding proteins (potentially chelators), involved in oxidation/reduction reactions and metabolism. Several genes were part of a secondary metabolite cluster predicted by both antiSMASH [48] (cluster 8.5) and SMURF (cluster 46) [45], which encodes for the PKS responsible for orsellinic acid, a precursor to many polyketides like lecanoric acid and the pigments F-9775A and B [49] as well as meroterpenoids like LL-Z1272β (Ilicicolin B) [50]. The most highly-expressed genes in this cluster were the efflux pump and an Omethyl transferase that could convert the orsellinic acid precursor to 3,5-dimethylorsellicinc acid, itself a precursor to ausitinol in A. nidulans [49]. AFLA_096040 and AFLA_096060 are two of the three genes found in the kojic acid biosynthesis pathway [51]. Despite co-culture having a 1.3 log 2 -fold change from Non-tox 17 alone, AFLA_096040 (an oxidoreductase) had the greatest RPKM value of all genes upregulated in the co-culture, from both Non-tox 17 and Tox 53. The RPKM value for AFLA_096040 was 16X greater than AFLA_096060, suggesting that even though there was a smaller log 2 -fold change difference (1.3 vs. 2.2) between co-culture and Non-tox 17, there would be more mRNA molecules for AFLA_096040. Similarly, when comparing the RPKM values of highly upregulated genes in Non-tox 17 and co-culture compared to Tox 53 to RPKM values of genes upregulated in co-culture than both Tox 53 and Non-tox 17, only five of the 13 (38%) genes had RPKM values greater than 50 when selected based on greater than 8-fold changes (Table 6b). Conversely, 14 of the 29 (48%) genes had RPKM values greater than 50 despite low log 2 -fold changes (1-3.2) or lack of significant difference from DeSeq2 analysis between co-cultures and both Non-tox 17 and Tox 53. This suggests that when selecting influential genes, both abundance and relative abundance should be considered.

Differential Expression of Imizoquin Biosynthesis Genes
Imizoquin biosynthesis was predicted to be enriched in Non-tox 17 compared to Tox 53; however, during close inspection of differential expression between Tox 53 and Non-tox 17 and co-cultures, none of the genes in imizoquin biosynthesis (imq) were highly differentially expressed (Table 3). Only 4 of 11 genes (AFLA_064230-064330) in the imq cluster [52] were found to be upregulated in both Non-tox 17 and co-cultures compared to Tox 53 at 72 h with log 2 -fold changes ranging between 1.8 and 4.8, (Table S1). However, upon comparing RPKM values there were differences in gene expression in between Tox 53, Non-tox 17 and co-cultures (Table 7). At both 30 and 72 h there was very little expression of genes in the imq cluster by Tox 53. However, it was found that at 30 h there is substantial expression of genes in Tox 53 from a secondary metabolite gene cluster (AntiSMASH cluster 1.1) adjacent to the imq cluster that may be associated with production of a toxic gliotoxin-like metabolite, likely aspirochlorine (AFLA_064340-AFLA_064610, acl) [53]. In several instances, there was less gene expression in co-cultures than Non-tox 17 though still greater than Tox 53, suggesting that imizoquin and aspirochlorine production is slightly attenuated in response to Tox 53.  ±2.8 ** 1.1, 21 1 Reads per kilobase per millions reads mapped (RPKM) were calculated for each 30 and 72 h Non-tox 17, Tox 53 and co-culture, reported as means and standard error. RPKM are measures of read abundance, which account for the total reads in a single replicate and individual gene length. This allows for better comparisons between treatments and genes rather than directly reporting reads. RPKM values for a single gene within 30 or 72 h with different numbers of asterisks * or ** are statistically different based on α ≤ 0.05. Values are color-scaled red for each gene and a darker shade indicates more positive values that are scaled to the maximum within an individual gene. 2 Gene ID is the AFLA_# gene number (i.e., 029700 = AFLA_029700) designated in the A. flavus NRRL 3357 reference genome. Gene ID font color red and blue correspond to the imizoquin and aspirochlorine gene clusters, respectively. 3 Chr is the chromosome where each gene is located. 4 SM is the secondary metabolite cluster where a gene is located. The 1st number is the SM cluster predicted by antiSMASH, and the 2nd number is the SM cluster predicted by SMURF.

Discussion
The current study investigated differences in gene expression during the biocontrol interaction between a non-aflatoxigenic Aspergillus flavus isolate (Non-tox 17) and a toxigenic isolate (Tox 53) which severely limits aflatoxin production by Tox 53. In support of the prevailing theory that competitive exclusion occurs by direct replacement of Tox with Non-tox isolates, Tox 53 grew much less both in mono-culture and co-culture with Non-tox 17. Since there was less Tox 53 RNA in co-cultures than expected based on its growth characteristics, in addition to significantly reduced transcription of both aflatoxin and CPA cluster genes, the data suggest Non-tox 17 is limiting gene expression and growth of Tox 53 as previously observed [32,33,37]. Non-tox 17 and other biocontrol isolates need to be further evaluated for their anti-fungal activity during co-culture. Limited growth of Tox 53 resulted in similar gene expression profiles between co-cultures and Non-tox 17. Expression of genes encoding proteins presumptively functioning in redox reactions, transcription factors and secreted proteins differed between Non-tox 17 and Tox 53 suggesting their possible roles in fungal growth and aflatoxin inhibition or degradation. Genes in select secondary metabolite clusters were either upregulated in Non-Tox 17 (asperfuranone and imizoquin) or further upregulated when co-cultured with Tox 53 (kojic acid and orsellinic acid). We are currently investigating if these secondary metabolites play a role in inhibition of aflatoxin production through both touch inhibition and recently reported contactless inhibition by chemicals secreted in culture filtrates from Non-tox (e.g., Non-tox 17) biocontrol isolates [37][38][39][40]. Several genes with statistical differences between samples but a log 2 -fold change less than 2 had very high RPKM (>100-1000) values, whereas genes with the highest log 2 -fold changes had RPKM values typically under 50. This suggests that using log 2 -fold changes can identify genes with high differential expression that are not expressed at high levels, therefore, RPKM values should also be considered to determine if differential expression of a gene will contribute more transcripts and potentially become more biologically influential. Based on our observations, biocontrol strains such as Non-tox 17 likely lower aflatoxin contamination by a combination of outcompeting and displacing Tox 53 and producing secondary metabolites, which may alter the redox state and extracellular environment or otherwise inhibit important cellular processes.
The majority of differentially expressed genes in the Non-tox 17 mono-culture and during co-culture were involved in oxidation and reduction reactions. It is hypothesized that aflatoxin is produced to minimize oxidative stress from the host plant's oxidative burst that occurs during fungal invasion or drought stress [36,54,55]. Several genes in the aflatoxin biosynthesis pathway are sources of reactive oxygen species (ROS) [54] and mutants and natural non-aflatoxigenic A. flavus and A. parasiticus strains are more sensitive to growth medium amended with H 2 O 2 [54,55]. Aflatoxin production is induced by H 2 O 2 and it was suggested that during aflatoxin synthesis, antioxidative enzymes scavenge H 2 O 2 from the environment and sequester ROS in vesicles, thereby alleviating oxidative stress in the fungus [54][55][56]. Alternatively, aflatoxin production may be a source of oxidative stress to the fungus due to a buildup of ROS, and it was shown that toxigenic isolates have greater glutathione S-transferase activity at the onset of aflatoxin production in comparison with Non-tox isolates [57,58]. Glutathione S-transferase activity should mollify oxidative stress resulting in a decrease in aflatoxin production [57,58]. Interestingly, most corn isolates are Non-tox or low toxin producers [42], provide the majority of biomass during co-infection of kernels with Tox isolates [33], and survive greater ROS defense responses from plants [36]. This suggests Non-tox isolates have alternative mechanisms to alleviate oxidative stress which may explain why we observed that most differentially expressed genes are involved in oxidation and reduction reactions. NRRL 21882, the Non-tox isolate in AflaGuard ® , differentially expressed more genes involved in oxidation and reduction than Tox isolates from H 2 O 2 -induced oxidative stress [56]. Further studies are needed to determine if Nontox isolates alter the redox environment, resulting in decreased aflatoxin production and invasion of plant tissue by Tox isolates.
In addition to limited growth of Tox 53 during co-culture with Non-tox 17, there was also reduced expression of aflatoxin biosynthesis pathway genes. Multiple Non-tox isolates downregulated aflR, aflJ, omtA, ordA, pksA, and vbs when co-cultured with Tox isolates [59]. During co-culture, it is impossible to rule out that inhibition of aflatoxin production is only due to outcompeting the Tox isolate by the Non-tox isolate since here Tox 53 grew substantially less than Non-tox 17. However, cell-free Non-tox media filtrates from A. flavus, including Non-tox 17 and A. oryzae, inhibited aflatoxin production [37][38][39][40]60] or degraded aflatoxin [41]. Genes in the early and middle portions of the aflatoxin biosynthesis pathway were downregulated in NRRL 3357 in response to A. oryzae filtrates [60]. The aflatoxin biosynthetic pathway-specific co-activator, aflS, was substantially downregulated, but there was not significantly less expression of the transcriptional activator aflR [60]. Contrary to our findings, there was greater expression of imizoquins and cyclopiazonic acid upon exposure to only culture filtrates [60]. These results indicate that Non-tox isolates may lower aflatoxin production by both displacement and inhibition of aflatoxin production through production of chemicals capable of downregulating expression of critical aflatoxin biosynthetic pathway genes.
Expression of several secondary metabolite cluster genes was either upregulated more in Non-tox 17 compared to Tox 53 and/or further upregulated in response to Tox 53 during co-culture. Some of these may be candidate compounds that interfere with aflatoxin production during the biocontrol interaction. Genes involved in kojic acid synthesis had the greatest RPKM values during co-culture. Kojic acid is commonly found in soy sauce and miso, and functions as an antioxidant that inhibits browning due to polyphenol oxidases in potatoes, apples and mushrooms [61]. It is also used in the cosmetic industry to lighten skin by inhibiting melanization [61]. During the biocontrol interaction, kojic acid may serve as an antioxidant resulting in less aflatoxin production by Tox isolates. Under elevated H 2 O 2 -induced oxidative stress, kojA expression increased in NRRL 3357 and NRRL 21,882 (AflaGuard), while other Tox and Non-tox isolates demonstrated normal levels of kojA expression [56]. In this manuscript, 30 and 72 h Non-tox 17 fungal cultures produced more transcripts than one-week-old cultures in Fountain et al. [56], suggesting transcription of genes in kojic acid synthesis may diminish with culture age, or Non-tox 17 produces much more kojic acid transcripts than other A. flavus isolates. Although the RPKM values were less, genes in the predicted orsellinic acid biosynthesis cluster (antiSMASH cluster 8.5, SMURF 46) [45] were also upregulated in response to Tox 53. The orsellinic acid gene in A. nidulans was turned on when the fungus physically interacted with the bacterium Streptomyces rapamycinicus [62], resulting in production of orsellinic acid and its derivatives: lecanoric acid, F-9775A, and F-9775B. A similar phenomenon could be occurring in our experiments (e.g., increased expression of the orsellinic acid biosynthesis cluster genes occurring after contact with other fungal hyphae). Predicted asperfuranone genes (SMURF cluster 20) [45,46] were expressed more by Non-tox 17 and co-cultures than Tox 53. Asperfuranone inhibits growth of small lung cancer cells and induces apoptosis [63], suggesting that asperfuranone could potentially inhibit growth of Tox 53. Finally, imizoquin cluster genes [52] were expressed at higher levels by Non-tox 17 at 30 and 72 h compared to Tox 53; co-cultures expressed intermediate levels. Imizoquins were downregulated in response to an isolate of Ralstonia solanacearum that produced a lipopeptide, which induced chlamydospore production in A. flavus [52,64]. Loss of imizoquin production delays spore germination and increases sensitivity to H 2 O 2 -induced oxidative stress [52] suggesting it is involved in spore germination and can act as an antioxidant. Continued expression of imizoquin cluster genes by Non-tox 17 may reduce aflatoxin production in Tox 53 by reducing oxidative stress. Future metabolomic studies will be used (1) to determine if kojic acid, orsellinic acid, asperfuranone, and imizoquins are produced by Non-tox 17 alone and in co-culture, and (2) to understand how they regulate growth and aflatoxin production of A. flavus.
Non-tox A. flavus isolates are widely used as biocontrol agents to effectively manage aflatoxin contamination of peanuts, corn, cottonseed and pistachios [15][16][17][18][19][20][21]. Although the biocontrol has been shown to work primarily via direct replacement of Tox isolates with Non-tox isolates [17,[25][26][27][28], as was confirmed in this manuscript, it is important to understand how Non-tox isolates molecularly and biochemically inhibit growth and toxin production of Tox A. flavus. Secondary metabolites previously found to be regulated in response to other microorganisms also produced different numbers of transcripts. Kojic acid and imizoquins, along with different individual genes, potentially alter aflatoxin production by serving as antioxidants. The greater antioxidant activity provided by kojic acid, imizoquins and other oxidation/reduction genes potentially gives the Non-tox a competitive advantage when infecting crops. Asperfuranone potentially acts in the biocontrol interaction by inhibiting growth. Future directions include determining if these chemicals are produced during the biocontrol interaction and assess their effects on A. flavus growth. If A. flavus chemicals (i.e., secondary metabolites) inhibit aflatoxin production, biocontrols should be evaluated for production of the most inhibitory chemicals, and then engineered to overproduce those chemicals or developed into a spray treatment mimicking the presence of Non-tox A. flavus.

Fungal Isolates
Aspergillus flavus Non-tox isolate 17, also named 07-S-3-1-6 (SRRC1588), was isolated from Louisiana corn field soil in 2007 [42] and is highly inhibitory to aflatoxin production [39,40]. Tox isolate 53 (SRRC1669) was isolated from Louisiana-grown, surfacesterilized corn kernels in 2003 [34], is highly toxigenic, and belongs to vegetative compatibility group RRS4 [42] originally isolated from corn kernels throughout Louisiana and along the Mississippi River in the US [65]. Tox 53 demonstrated the importance of physical interaction for toxin inhibition during a previous biocontrol interaction [34]. Both isolates are deposited in an accessible culture collection at the USDA-ARS Southern Regional Research Center (SRRC) in New Orleans Louisiana.

Biocontrol Interaction Cultural Experimental Design
RNA and aflatoxin were extracted from Tox 53 and Non-tox 17 isolates grown alone in mono-culture and together in co-cultures to understand how their gene expressions and aflatoxin production change in response to one another during the biocontrol interaction (i.e., co-cultures). Fresh conidia from 5-day old Tox 53 and Non-tox 17 cultures (grown in dark at 30 • C on 5% V8 juice, 2% agar, pH 5.2 medium) were suspended in glucose salts liquid medium (GS), with citrate buffer (pH 4) [32,34,43,66], at 1 × 10 5 conidia/mL. Citrate buffer was added to maximize aflatoxin production, limit aflatoxin degradation and limit detrimental effects to fungal growth by low pH [39,40,43,66]. Buffered GS medium consisted of one part water, 2 parts 2.5 X salts and 2 parts 2.5 X glucose, mixed after autoclaving. One liter of 2.5 X salts consisted of 8.76 g (NH 4 )SO 4 (Table 8). To obtain enough tissue for RNA extraction at 30 h, a biological replicate consisted of 9 Petri-dish cultures. For all other time points, only one Petri-dish was used per biological replicate. The isolate mono-cultures, co-cultures and time points grew in separate boxes (five Petri dishes per box) within a darkened 25 • C incubator to minimize potential effect of differential volatile production between isolates on aflatoxin production [66,67]. RNA was extracted from 30 and 72 h cultures and aflatoxin was extracted from 30, 72 and 96 h cultures.

Aflatoxin Extraction and Quantification
At 30 h, 100 µL of liquid medium from each dish (n = 9) per biological rep was added to HPLC grade methanol. At 72 and 96 h, 500 µL of medium from each biological replicate (a single Petri dish) was added to 500 µL of HPLC grade methanol. Extracts were filtered through 200 mg basic alumina (58Å, 60-mesh powder, 11503-A1, Alfa Aesar, Tewksbury MA) in 1.5 mL polypropylene columns with 20 µm polyethylene frits [68]. 10 µL of each sample was separated in a Waters e2695 Separations Model HPLC (Waters Corp., Milford, MA) using a Nova-Pak C18 4 µm, 3.9 mm × 150 mm column held at 38 • C with an isocratic solvent system (37.5 Methanol: 62.5 water at a 0.8 mL/min) coupled to a PHRED photochemical reactor cell (Aura Industries Inc., New York, NY, USA). After separation and photolytic derivatization, a 2475 FLR Detector (Waters Corp., Milford, MA, USA) was used to detect and quantify aflatoxin B 1 (365 nm Ex, 440 nm Em) [69,70]. Run time was 17 min with aflatoxin B 1 peak eluting at~13.5 min. Empower software (Waters Corp., Milford, MA, USA) was used to integrate the aflatoxin B 1 peak. Aflatoxin quantity was calculated based on a calibration curve calculated from 4 replicates of standards with 1, 5, 50, 500 and 1000 ng/mL aflatoxin B 1 [70]. Aflatoxin B 1 minimum level of detection was <0.05 ppb and minimum quantification from standard curve was 1 ppb.

Whole Fungal Mycelia Harvest and RNA Extraction
At 30 and 72 h, mycelia and medium were removed from the Petri dishes and centrifuged at 8000× g for 5 min at 4 • C. Thirty-hour tissues from nine plates per biological rep were pooled and centrifuged a second time for 5 min. Excess medium was removed by carefully blotting mycelia on chromatography paper. The tissue was added to a pre-weighed microcentrifuge tube (to calculate wet weight) and flash frozen with liquid nitrogen. RNA extraction was performed according the manufacture's guidelines for the Spectrum™ Plant Total RNA Kit (STRN250, Sigma-Aldrich, St. Louis, MO, USA) and the On-column Dnase I Digestion Set (DNASE70, Sigma-Aldrich, St. Louis, MO, USA) with a couple of modifications. All tissue from a single biological replicate was ground directly in lysis buffer (100 mg mycelia/500 µL lysis buffer). A few 30 h cultures had less than 100 mg, which were still ground in 500 µL lysis buffer. For each sample, 500 µL was retained for RNA extraction. Binding buffer was increased to 750 µL due to inefficient RNA extraction from the residual medium.

RNA Sequencing and Analysis
Three RNA extracts per experimental condition were sequenced by NC State University's Genomic Sciences Laboratory using an Illumina NextSeq 500, which generated 150 bp paired-end reads. Sequencing reads were submitted to NCBI's Sequence Read Archive and can be accessed under BioProject ID PRJNA764255. Sequence reads were trimmed to remove adapters and low-quality sequences using BBDuk [71]. Sequencing reads were mapped to the A. flavus NRRL 3357 genome (JCVI-afl1-v2.0 assembly, (https: //www.ncbi.nlm.nih.gov/assembly/GCF_000006275.2/#/st, accessed on 8 April 2019) using STAR v2.6.1 [72]. Reads mapped to exons were counted using featureCounts v1.6.0 [73] followed by differential expression testing of normalized reads using a generalized linear model with log link and a negative binomial distribution within DESeq2 [47]. Genes were removed if they did not have at least 10 reads in 3 or more samples. Genes were considered differentially expressed if the pairwise comparison by DESeq2 software p-value was less than 0.05 and if there was a log 2 -fold change greater than 2 [47]. To make the principal component analysis (PCA) plot, regularized log counts were produced with the DESeq2 s rlog function and the option "blind = TRUE" was set [47]. These were used as input to the plotPCA function in DESeq2 [47]. In order to quantify the fraction of RNA-seq reads contributed by each strain, variants were called using Freebayes [74]. Variants that were different between Non-tox 17 and Tox 53 were used in a custom python script utilizing the pysam library (https://github.com/pysam-developers/pysam, accessed on 8 April 2019)) and SAMtools [75] to assign reads from the mixed cultures to each strain. Functional enrichment analysis was performed with the enrichment function in BC3NET R package [76], which uses a one-sided Fisher's Exact test with the Benjamini and Hochberg adjustment [77].
Excel version 2102 (Microsoft corp., Redmond, WA) was used to sort pairwise log 2fold differential gene expression testing from DESeq2 for each pairwise comparison of Non-tox 17 vs. Tox 53, Co-culture vs. Tox 53, and Co-culture vs. Non-tox 17 at 30 and 72 h. Genes that were overexpressed in biocontrol isolate Non-tox 17 were selected if the log 2 -fold change was ≥8. Genes that were further upregulated in Non-tox 17 during co-culture were selected if Co-culture vs. Tox 53 and Co-culture vs. Non-tox 17 log 2 -fold changes were >1. Additionally, further upregulated genes were selected if the differences between Co-culture vs. Tox 53 and Non-tox 17 vs. Tox 53 were log 2 -fold changes at least 1. Since the latter selection criterion was not statistically different based on DESeq2 analysis of normalized reads, generalized linear models were calculated to compare gene expression for each of those genes using the logit (log odds, i.e., (proportion reads (proportion (p) reads aligned to gene X/(p reads not aligned to gene X)) link for binomial data with SAS version 9.4 (SAS Institute, Cary, North Carolina). The fixed effects were culture type (Non-tox 17, Tox 53 and Co-culture) and culture age (30 and 72 h). The response variable was reads/total reads. Treatments were separated by post hoc comparison of odds with a difference of least squares means at α ≤ 0.05. Excel was also used to calculate reads per kilobase per million mapped reads (RPKM) for genes selected by sorting. RPKM for gene X = (1 × 10 9 ) (read mapped to gene X)/(gene X length bp) (total reads mapped) [47,78].

Other Data Analysis
Generalized linear models estimated multivariate analysis of variance to compare biomass, total RNA and aflatoxin B 1 between treatments using SAS. To address issues with normality, aflatoxin values were log transformed. In each model, fixed effects were either isolate growing alone or in co-culture, extraction time, and their interaction. Means were separated by post hoc comparison with a difference of least squares means at α ≤ 0.05. To determine if the number of reads which uniquely aligned to Non-tox 17 and Tox 53 during co-culture was similar to the expected ratio based on biomass and RNA production of each isolate growing separately, generalized linear models estimated multiple categorical data analysis (i.e., multiple contingency tables) using logit link and binomial distribution with SAS. Log odds (p Tox 53/p Non-Tox 17) were calculated within the model by inputting the events (either number of unique reads, biomass or total RNA of the Non-tox) and dividing by trials (total number of reads, sum of biomass and total RNA of Non-Tox 17 and Tox 53 isolates). Odds were separated by post hoc comparison with a difference of least squares means at α ≤ 0.05.

Informed Consent Statement: Not applicable.
Data Availability Statement: Sequencing reads were submitted to NCBI's Sequence Read Archive and can be accessed under BioProject ID PRJNA764255 and can be accessed at https://www.ncbi.nlm. nih.gov/sra. The remaining data is contained within the article and Supplementary Materials.