Responses of Aspergillus flavus to Oxidative Stress Are Related to Fungal Development Regulator, Antioxidant Enzyme, and Secondary Metabolite Biosynthetic Gene Expression

The infection of maize and peanut with Aspergillus flavus and subsequent contamination with aflatoxin pose a threat to global food safety and human health, and is exacerbated by drought stress. Drought stress-responding compounds such as reactive oxygen species (ROS) are associated with fungal stress responsive signaling and secondary metabolite production, and can stimulate the production of aflatoxin by A. flavus in vitro. These secondary metabolites have been shown to possess diverse functions in soil-borne fungi including antibiosis, competitive inhibition of other microbes, and abiotic stress alleviation. Previously, we observed that isolates of A. flavus showed differences in oxidative stress tolerance which correlated with their aflatoxin production capabilities. In order to better understand these isolate-specific oxidative stress responses, we examined the transcriptional responses of field isolates of A. flavus with varying levels of aflatoxin production (NRRL3357, AF13, and Tox4) to H2O2-induced oxidative stress using an RNA sequencing approach. These isolates were cultured in an aflatoxin-production conducive medium amended with various levels of H2O2. Whole transcriptomes were sequenced using an Illumina HiSeq platform with an average of 40.43 million filtered paired-end reads generated for each sample. The obtained transcriptomes were then used for differential expression, gene ontology, pathway, and co-expression analyses. Isolates which produced higher levels of aflatoxin tended to exhibit fewer differentially expressed genes than isolates with lower levels of production. Genes found to be differentially expressed in response to increasing oxidative stress included antioxidant enzymes, primary metabolism components, antibiosis-related genes, and secondary metabolite biosynthetic components specifically for aflatoxin, aflatrem, and kojic acid. The expression of fungal development-related genes including aminobenzoate degradation genes and conidiation regulators were found to be regulated in response to increasing stress. Aflatoxin biosynthetic genes and antioxidant enzyme genes were also found to be co-expressed and highly correlated with fungal biomass under stress. This suggests that these secondary metabolites may be produced as part of coordinated oxidative stress responses in A. flavus along with antioxidant enzyme gene expression and developmental regulation.

The infection of maize and peanut with Aspergillus flavus and subsequent contamination with aflatoxin pose a threat to global food safety and human health, and is exacerbated by drought stress. Drought stress-responding compounds such as reactive oxygen species (ROS) are associated with fungal stress responsive signaling and secondary metabolite production, and can stimulate the production of aflatoxin by A. flavus in vitro. These secondary metabolites have been shown to possess diverse functions in soil-borne fungi including antibiosis, competitive inhibition of other microbes, and abiotic stress alleviation. Previously, we observed that isolates of A. flavus showed differences in oxidative stress tolerance which correlated with their aflatoxin production capabilities. In order to better understand these isolate-specific oxidative stress responses, we examined the transcriptional responses of field isolates of A. flavus with varying levels of aflatoxin production (NRRL3357, AF13, and Tox4) to H 2 O 2 -induced oxidative stress using an RNA sequencing approach. These isolates were cultured in an aflatoxin-production conducive medium amended with various levels of H 2 O 2 . Whole transcriptomes were sequenced using an Illumina HiSeq platform with an average of 40.43 million filtered paired-end reads generated for each sample. The obtained transcriptomes were then used for differential expression, gene ontology, pathway, and co-expression analyses. Isolates which produced higher levels of aflatoxin tended to exhibit fewer differentially expressed genes than isolates with lower levels of production. Genes found to be differentially expressed in response to increasing oxidative stress included antioxidant enzymes, primary metabolism components, antibiosis-related genes, and secondary metabolite biosynthetic components specifically for aflatoxin, aflatrem, and kojic acid. The expression of fungal development-related genes including aminobenzoate degradation genes and conidiation regulators were found to be regulated in response to increasing stress. Aflatoxin biosynthetic genes and antioxidant enzyme genes were also found to be

INTRODUCTION
The contamination of crops with aflatoxin, a carcinogenic secondary metabolite of the facultative plant parasite Aspergillus flavus (Guo et al., 1996), is a threat to human health, global food safety and security (Williams et al., 2010;Guo et al., 2012;Torres et al., 2014;Andrade and Caldas, 2015). Aflatoxin contamination of staple and dietary supplemental crops such as maize and peanut result in both losses in crop value in international trade due to restrictions on aflatoxin content (Matumba et al., 2015;Wu, 2015), and negative impacts in human and animal health (Williams et al., 2004(Williams et al., , 2010Kew, 2013). These concerns are the impetus for investigations into the biology of this organism and its interactions with host plants related to aflatoxin contamination (Diener et al., 1983(Diener et al., , 1987Amaike and Keller, 2011;Guo et al., 2012;Fountain et al., 2014).
The aflatoxin biosynthetic pathway has been well characterized in A. flavus and in other aflatoxigenic species of Aspergillus such as A. parasiticus, and sterigmatocystin producing species such as A. nidulans (Amaike and Keller, 2011). Aflatoxin biosynthesis is encoded by a cluster of 25 genes which has been highly conserved among Aspergillus spp. and has been well characterized (Yu et al., 2004). While the biosynthetic mechanisms involved in aflatoxin production have been well characterized, little is known regarding the biological role of aflatoxin in A. flavus or other Aspergillus spp. Secondary metabolites produced by soil-dwelling fungi exhibit various biological activities including fungivory resistance, stress tolerance, and quorum sensing (Reverberi et al., 2008(Reverberi et al., , 2010Roze et al., 2013).
Recent studies have shown that reactive oxygen species (ROS) and their reactive products such as peroxidized lipids (oxylipins) are required for the production of aflatoxin and can stimulate aflatoxin production if applied in vitro (Jayashree and Subramanyam, 2000). Induction of oxylipin and ROS accumulation in A. flavus mycelia through peroxisome proliferation has also been linked with increased aflatoxin production and antioxidant enzyme activity (Reverberi et al., 2012). Similarly, several studies have also been performed examining the effects of antioxidants on the growth and aflatoxin production of Aspergilli. For example, phenolic compounds such as caffeic acid tannic acid derived from tree nuts have been shown to inhibit aflatoxin production in A. flavus . Other synthetic phenolic compounds such as butylated hydroxyanisole (BHA) and propyl paraben (PP) have also been found to have a similar effect as a function of medium pH and water activity (Nesci et al., 2003;Passone et al., 2005). Treatment with BHA as also been shown to inhibit sclerotial differentiation in A. flavus. Oxidative stress responsive signaling mechanisms have been found to be involved in the regulation of aflatoxin production such as the stress responsive transcription factors AtfB, AP-1, and VeA (Reverberi et al., 2008;Sakamoto et al., 2008;Roze et al., 2011Roze et al., , 2013Hong et al., 2013a,b;Baidya et al., 2014). Also, VeA along with VelB and LaeA form the Velvet protein complex to regulate both secondary metabolite production and reproductive development in A. flavus and other Aspergillus spp. (Park et al., 2012). The consumption and/or generation of ROS have also been found to co-localize to vesicles known as aflatoxisomes where the final phases of aflatoxin production are carried out (Chanda et al., 2009(Chanda et al., , 2010Roze et al., 2015).
Because of the close association of ROS and aflatoxin production, it has been hypothesized that aflatoxin production may serve as a component of oxidative stress alleviation mechanisms employed by Aspergillus spp. in addition to antioxidant enzymes, altered carbon metabolism, and the production of other secondary metabolites (Narasaiah et al., 2006;Roze et al., 2013;Fountain et al., 2014Fountain et al., , 2016. The tolerance of A. flavus and A. parasiticus isolates to oxidative stress has been shown to be correlated with their levels of aflatoxin production. For example, Roze et al. (2015) showed that conidia of isolates with higher levels of aflatoxin production also exhibited greater viability when cultured in ROS-amended medium.
This correlation between ROS and aflatoxin production has also lead to the hypothesis that host plant-derived ROS and oxylipins may function in the host-parasite interaction between A. flavus or other Aspergilli and their host plants. Host plant resistance to aflatoxin contamination have been identified, and this resistance has been found to be heavily influenced by environmental stresses, particularly drought stress (Diener et al., 1983(Diener et al., , 1987Williams, 2006;Guo et al., 2008;Holbrook et al., 2009;Jiang et al., 2012;Kebede et al., 2012;Pandey et al., 2012;Fountain et al., 2014). Interestingly, these aforementioned ROS and oxylipins have also been shown to accumulate in the tissues of host plants during drought stress, and their levels have been correlated with aflatoxin contamination resistance (Gao et al., 2009;Yang et al., 2015Yang et al., , 2016. Host plant tissue antioxidant enzyme activity and capacity has also been found to be correlated with reduced A. flavus growth and aflatoxin production in model species such as buckwheat (Chitarrini et al., 2014).
In our previous studies, isolates of A. flavus were found to exhibit different degrees of oxidative stress tolerance which appeared to correlate with their aflatoxin production capability suggesting that aflatoxin production may contribute to stress tolerance . However, the observed stress tolerance of the toxigenic isolates cultured in aflatoxin nonconducive medium was reduced yet comparable indicating that factors in addition to aflatoxin production also contributed to the observed differences. In order to better understand the differences in isolate-specific responses to oxidative stress, and to further explore the potential role of aflatoxin production in stress alleviation in A. flavus, we examined the global transcriptional responses of several isolates of A. flavus to increasing oxidative stress, which has been summarized in an overview-type publication . Here, we present a detailed analysis of changes in the transcriptomes of different toxigenic A. flavus isolates with distinguished aflatoxin production capabilities to increasing oxidative stress in an aflatoxin conducive culture medium. By examining the oxidative stress responses of A. flavus, the biological role of aflatoxin production in stress responses and in competition with other soil-dwelling organisms may be better understood.

MATERIALS AND METHODS
The methodologies presented here are adapted from Fountain et al. (2016) and describe the data generation and analysis procedures for the toxigenic isolates examined in this study.

Isolate Culture Conditions
The A. flavus toxigenic isolates utilized in this study were obtained as previously described . The isolates AF13, NRRL3357, and Tox4 were initially cultured on V8 agar (20% V8, 1% CaCO 3 , 3% agar) at 32 • C for 5 days. Conidia were collected from the cultures by washing the plates with sterile 0.1% (v/v) Tween 20. This conidial suspension (∼4.0 × 10 6 conidia/mL) was used to inoculate liquid cultures of yeast extract sucrose (YES; 2% yeast extract, 1% sucrose) medium containing various concentrations of hydrogen peroxide (H 2 O 2 , 3% stabilized solution). For AF13 and Tox4, the YES medium was supplemented with 0, 10, and 25 mM H 2 O 2 representing a control, moderately high, and high levels of stress, respectively. For NRRL3357, the YES medium was supplemented with 0, 10, and 20 mM H 2 O 2 due to the lower concentration of H 2 O 2 the isolate could tolerate . The experiment was performed in 125 mL Erlenmeyer flasks containing 50 mL H 2 O 2 amended YES medium and 100 µL conidial suspension plugged with a sterile cotton ball. The cultures were incubated at 32 • C for 7 days in the dark with two biological replicates for each isolate/treatment combination. Mycelia were then harvested from each culture, immediately frozen in liquid nitrogen, and stored at −80 • C.

RNA Isolation
The harvested mycelia were homogenized in a chilled mortar and pestle to a fine powder, and then used for total RNA isolation. Total RNA was isolates using an RNeasy Plant Mini Kit with DNase digestion according to the manufacturer's instructions (Qiagen, Hilden, Germany). The isolated RNA was then quantified using a Nano-Drop ND1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and the RNA integrity numbers (RINs) for each sample were validated using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Samples with an RIN ≥ 5 were used for RNA sequencing.

Library Construction and Illumina Sequencing
The cDNA libraries for each sample and biological replicate were generated from 1 µg of total RNA. A total of 18 libraries were generated for the toxigenic A. flavus isolate samples and used for transcriptome sequencing using an Illumina TruSeq RNA Sequencing Kit according to the manufacturer's instructions (Illumina, San Diego, CA, USA). Following quantitation using a Qubit 2.0 fluorometer (Thermo Scientific), and validation using an Agilent 2100 Bioanalyzer (Agilent), the libraries were used for cluster generation using a cBot (Illumina), and pairedend sequencing using a HiSeq 2500 platform according to the manufacturer's instructions (Illumina).

Bioinformatics Analysis
Initial quality checks on raw sequencing reads were performed using FastQC v0.11.2 with low quality reads being removed using Trimmomatic v0.32, along with rRNA contamination following alignment with the SILVA database. The remaining, high quality reads were used for differential expression analysis using the tuxedo protocol. Briefly, alignment of filtered reads to the A. flavus NRRL3357 reference genome (GCF_000006275.2) was done using tophat2 v2.0.13 and bowtie2 v2.2.4. Cufflinks v2.2.1 and cuffdiff were used to assemble the transcripts and determine transcript abundance in terms of Fragments Per Kilobase of exon per Million fragments mapped (FPKM). A gene was considered significantly differentially expressed when |log 2 (fold change)| ≥ 2 with an adjusted p ≤ 0.05. Tophat2 alignments were used to perform RABT (reference annotation based transcript) assembly for both genes and isoforms. The assemblies were then compared and merged using cuffmerge and used for analysis.
The assembled transcripts were annotated using a standalone blast 2.2.30+ and analyzed through Blast2GO and KEGG for gene ontology (GO) and pathway analysis, respectively. Heatmaps of select secondary metabolites biosynthesis genes and principal components analysis (PCA) of the gene expression profiles of the isolates was performed with multiple experiment viewer (MeV) v4.9.0. For promoter analysis of select genes, upstream gene sequences containing conserved transcription factor binding domains were obtained using FungiDB (http://fungidb.org/) followed by sequence analysis using MEME 4.11.2 (Bailey et al., 2009;Stajich et al., 2012). Co-expression analysis was performed using the weighted correlation network analysis (WGCNA) package in R v3.3.0 (Langfelder and Horvath, 2008).

Quantitative RT-PCR Validation
In order to validate the RNA sequencing results, we selected several genes for expression analysis on a subset of samples using quantitative RT-PCR (qPCR). Using total RNA remaining following library construction, we synthesized cDNA using a Taqman reverse transcriptase kit (Thermo-Fisher) according to the manufacturer's instructions. Reverse transcription was carried out using a PTC-200 thermal cycler (Bio-Rad, Hercules, CA, USA) with the following cycling parameters: 25 • C for 10 min, 48 • C for 30 min, and 95 • C for 5 min. The qPCR was then performed using a 20 µl reaction volume containing: 1X SYBR Green PCR Master Mix (Thermo-Fisher), 0.4 µM forward primer, 0.4 µM reverse primer, and 25 ng cDNA template. In this study, β-tubulin was used as a housekeeping control, and primer sequences for this and the other genes examined can be found in Table S1. The reactions were carried out using an ABI 7500 platform (Thermo-Fisher) using the following cycles: 50 • C for 2 min, 95 • C for 10 min, and 40 amplification cycles of 95 • C for 15 s, and 60 • C for 1 min. Dissociation curves were performed at the end of each cycle to identify possible primer dimerization or off-target amplification. With the obtained threshold cycle (Ct) values, relative gene expression was calculated using a modified Livak method where Relative Expression = 2 Ct and C t = C t Housekeeping Gene (β-tubulin)-C t Target (Livak and Schmittgen, 2001).

Transcriptome Sequencing
In order to examine the transcriptional responses of aflatoxigenic isolates of A. flavus to oxidative stress, we cultured three isolates (AF13, NRRL3357, and Tox4) in aflatoxin production-conducive medium amended with H 2 O 2 at different concentrations. Whole transcriptome sequencing of these tissues yielded a total of 1.32 billion raw sequencing reads from 18 cDNA libraries with an average of 73.17 million reads per library. Of the 18 libraries sequenced, one contained a relatively low number reads and was excluded from the analysis. Initial quality filtration resulted in the removal of 9-10% of reads followed by further removal of rRNA contamination. Following filtration, a total of 687.33 million filtered reads were obtained with an average of 40.43 million reads per library. For all the libraries, an average of 92.64% of filtered reads mapped to the A. flavus NRRL3357 reference genome. For the individual isolates, an average of 91.66, 94.52, and 91.58% of filtered reads from AF13, NRRL3357, and Tox4, respectively, mapped to the reference genome.

Differential Expression Is Correlated with Oxidative Stress Tolerance
Following transcriptome assembly and alignment, differential expression analyses were performed to examine the oxidative stress responses of the isolates. Genes were considered significantly differentially expressed if they exhibited a |log 2 (fold change)| ≥ 2 and adjusted p ≤ 0.05. The Tox4 isolate exhibited 4 and 29 DEGs when comparing the control with 10 mM H 2 O 2 and 25 mM H 2 O 2 treatments, respectively, and 57 when comparing the two treatments ( Table 1). The AF13 isolate exhibited 6 and 122 DEGs when comparing the control with 10 mM H 2 O 2 and 25 mM H 2 O 2 treatments, respectively, and 85 when comparing the two treatments ( Table 1). The NRRL3357 isolate exhibited 53 and 117 DEGs when comparing the control with 10 mM H 2 O 2 and 20 mM H 2 O 2 treatments, respectively, and 112 when comparing the two treatments ( Table 1).
The numbers of significant DEGs for the isolates were then correlated with the previously observed maximum H 2 O 2induced oxidative stress tolerance levels observed for each isolate in our previous study . The number of DEGs exhibited a strong, negative correlation with the previously observed tolerance levels ( Figure S1) when comparing the control and the 10 mM H 2 O 2 treatment (r = −0.979), and when comparing the two treatment conditions (r = −0.958). The correlation was not as strong, however, when comparing the control and the 20 or 25 mM H 2 O 2 treatment (r = −0.658). Interestingly, the isolates which exhibited the lower numbers of differentially expressed genes, especially under 10 mM H 2 O 2 , also tended to produce higher levels of aflatoxin (Table 1; Fountain et al., 2015). However, while the isolate with the highest level of H 2 O 2 tolerance, Tox4, did exhibit the fewest number of DEGs, it produces similar levels of aflatoxin to AF13 which exhibited slightly lower H 2 O 2 tolerance and almost four times the number of DEGs at 25 mM H 2 O 2 ( Table 1). In addition, comparison of the overall expression profiles shows a clear segregation by isolate indicative of the differences in isolate-specific gene expression patterns due to background expression differences between isolates and responses to increasing stress (Figure 2, Table S2).

Secondary Metabolite Gene Expression Is Regulated in Response to Oxidative Stress
Genes involved in the biosynthesis of aflatoxin were differentially expressed in the toxigenic isolates with distinguished aflatoxin production capabilities in response to increasing levels of oxidative stress ( Figure 3A). The most significant changes in expression mainly came in the isolate with the least tolerance to high levels of H 2 O 2 than the other isolates included in this study, NRRL3357, which also produces only a moderate level of aflatoxin under elevated oxidative stress . In NRRL3357, the genes encoding Omethyltransferase A (omtA), versicolorin B synthase (vbs), and versicolorin dehydrogenase/ketoreductase (ver-1) were exhibited log 2 (fold change) > 3.61 when comparing 0 and 10 mM H 2 O 2 treatments, and > 4.85 when comparing 0 and 20 mM H 2 O 2 (Dataset S3). Additional genes encoding polyketide FIGURE 1 | Venn diagrams of genes expressed with an FPKM ≥ 2. Genes exhibiting an FPKM ≥ 2 in at least one isolates and treatment were considered to be expressed. Genes expressed within all isolates and treatments (A), and between treatments within the NRRL3357 (B), AF13 (C), and Tox4 (D) isolates were compared using Venn diagrams generated with Venny 2.1.0. synthase (pksA), norsolorinic acid ketoreductase (nor-1), and a cytochrome p450 monooxygenase (aflV/cypX) also exhibited log 2 (fold change) > 2.32 when comparing the 0 and 20 mM H 2 O 2 treatments (Dataset S3). In AF13, the cypX, omtA, and ver-1 genes along with the genes encoding the versicolorin B desaturase p450 monooxygenase (verB), and the regulatory transcription factor aflS/aflJ were upregulated when comparing the 0 and 20 mM H 2 O 2 treatments (Dataset S4). None of these genes were found to be differentially expressed when comparing 0 and 10 mM H 2 O 2 (Dataset S4).
the p450 monooxygenase moxY, oxidoreductase B (ordB), and the hypothetical protein hypA were expressed with an FPKM ≥ 2 only in AF13 and NRRL3357 (Dataset S2). The remaining aflatoxin pathway genes were found to be expressed in AF13, NRRL3357, and Tox4 (Dataset S2). In addition to aflatoxin biosynthesis genes, genes involved in the biosynthesis of two additional secondary metabolites, aflatrem and kojic acid, were found to be differentially expressed in response to increasing oxidative stress ( Figure 3B). In NRRL3357, the aflatrem biosynthesis genes encoding prenyl transferase (atmC) and geranylgeranyl diphosphate synthase (atmG) (Nicholson et al., 2009) were significantly upregulated in when comparing 0 and 10 mM H 2 O 2 (Dataset S3). The atmG gene was however, not upregulated when comparing 0 and 20 mM H 2 O 2 , though it is significantly reduced when comparing 10 and 20 mM H 2 O 2 indicating a dose-dependent response (Dataset S3). In contrast, in AF13 the dimethylallyl tryptophan synthase gene (atmD), atmC, and atmG were downregulated when comparing 0 and 25 mM H 2 O 2 while no aflatrem biosynthesis genes were significantly differentially expressed when comparing 0 and 10 mM H 2 O 2 (Dataset S4). Tox4, however, showed no significant differences in aflatrem biosynthesis gene expression when comparing the control with either treatment, though comparing 10 and 25 mM H 2 O 2 did show a significant reduction in atmG expression (Dataset S5).
For kojic acid biosynthetic genes (Terabayashi et al., 2010), the kojA and kojR/T genes were constitutively expressed in all the isolates and did not vary significantly in expression in response to increasing stress ( Figure 3C). The synaptic vesicle transporter SVOP, which is also involved in the biosynthesis of kojic acid (Terabayashi et al., 2010), was found to be differentially expressed in the examined isolates. In NRRL3357, the SVOP gene was upregulated by a log 2 (fold change) of 3.63 when comparing 0 and 10 mM H 2 O 2 , but was not significantly different between 0 and 20 mM (Dataset S3). For AF13, the SVOP gene was also upregulated by a log 2 (fold change) of 2.47 when comparing 0 and 25 mM H 2 O 2 , but was not significantly different between 0 and 10 mM H 2 O 2 (Dataset S4). For Tox4, the SVOP gene was not significantly different between the control and the treatments, but was found to be downregulated when comparing 10 and 25 mM H 2 O 2 indicating a slight reduction in response to higher levels of stress (Dataset S5). The expression of kojA and SVOP along with select aflatoxin and aflatrem genes were further confirmed with qPCR which were generally correlated with the FPKM data obtained by RNA sequencing (Figure 4).

Stress Responsive and Related Biochemical Pathways Regulated by Increasing Stress
In addition to the production mechanisms regulating secondary metabolite biosynthesis, additional stress responsive and associated biochemical pathway genes were also regulated in the isolates in response to increasing stress (Figure 5). Among the stress responsive genes, antioxidant and oxidase enzyme-encoding genes were among the most commonly regulated. In the NRRL3357 and AF13, thioredoxin peroxidase and thioredoxin reductase along with several cytochrome p450 monooxygenase genes were upregulated in response to increasing levels of oxidative stress (Datasets S3, S4). Interestingly, in Tox4 there were no antioxidant or oxidase genes regulated by 10 mM H 2 O 2 , and only one monooxygenase gene was regulated at 20 mM (Dataset S5).
In addition to these antioxidant and oxidase genes, genes related to fungal catabolism and antibiosis tolerance were also regulated by stress, particularly in AF13. Comparing 0 and 10 mM H 2 O 2 , two multiple drug resistance protein-encoding genes were upregulated and expressed at a high level (FPKM's of 1636.6 and 7697.7, respectively; Dataset S4). These genes were further upregulated when comparing 0 and 25 mM H 2 O 2 (FPKMs of 2627.0 and 15,066.5, respectively), along with a third expressed at a lower level (Dataset S4). When comparing 0 and 25 mM H 2 O 2 , AF13 also exhibited the FIGURE 4 | Relative expression and FPKM associations for select secondary metabolite genes. Real-time PCR was used to determine the relative expression of six select secondary metabolite genes: pksA and vbs (aflatoxin), atmG and atmQ (aflatrem), and kojA and SVOP (kojic acid). Relative expression is indicated by the blue bars scaled to the left y-axis. The RNA sequencing-derived FPKM data is indicated by the red line scaled to the right y-axis. The plots are organized into three columns corresponding to data obtained from NRRL3357, AF13, Tox4, respectively, with each plot containing data from the 0, 10, and 20/25 mM H 2 O 2 treatments in increasing order. downregulation of a number of catabolic enzyme-encoding genes including acetyl xylan esterase, alpha-1,3-glucanase, aspergillopepsin 2 and F, class III and V chitinases, chitosanase, and penicillolysin/deuterolysin (Dataset S4). NRRL3357 also exhibited a similar downregulation in the chitin catabolic genes encoding beta-N-acetylhexosaminidase (nagA) and beta-N-hexosaminidase when comparing 0 and 20 mM H 2 O 2 (Dataset S3). Iron metabolic genes were also regulated in NRRL3357 and AF13 in response to increasing stress. In NRRL3357, the siderophore biosynthesis acetylase gene aceI was upregulated when comparing 0 and 20 mM H 2 O 2 , and in AF13, a putative transferrin receptor gene was downregulated when comparing 0 and 25 mM H 2 O 2 (Datasets S3, S4). NRRL3357 also exhibited an upregulation in cinnamoyl-CoA reductase gene expression when comparing 0 and 20 mM H 2 O 2 (Dataset S3). Tox4 exhibited regulation of two phosphate signaling genes in response to increasing stress, but none of the aforementioned genes (Dataset S5).

Fungal Development-Related Genes Are Regulated by Oxidative Stress
Genes involved in the development of A. flavus and other Aspergillus spp. were also differentially expressed among the isolates. In NRRL3357, the gene encoding the Cis 2 His 2 (C 2 H 2 ) transcription factor flbC was upregulated along with a conidiation-specific family protein gene and a duf221 domain protein gene (Dataset S3). AF13 also exhibited an upregulation of a conidial hydrophobin gene in response to increasing stress (Dataset S4). In addition, genes involved in benzoate degradation were also regulated in response to stress in the isolates. NRRL3357 and AF13 both exhibited the upregulation of genes involved in benzoate degradation including a multiple inositol polyphosphate phosphatase and a benzoate 4monooxygenase (Datasets S3, S4). Conversely, Tox4 exhibited a downregulation of an amidohydrolase family protein gene and a multiple inositol polyphosphate phosphatase gene (Dataset 5).

bZIP Transcription Factors and Promoter Occurrence among Differentially Expressed Genes
The orthologs of two bZIP transcription factor genes, atfA and atf21, were previously shown to regulate the expression of genes involved in both the biosynthesis of aflatoxin and sterigmatocystin, and oxidative stress responsive mechanisms (Roze et al., 2011). Here, atf21 and atfA exhibited marginal increases in expression in NRRL3357 and AF13 under high levels of H 2 O 2 stress ( Figure S2). The ortholog of atf21 in A. parasiticus, atfB, has been well characterized and has been shown to bind a conserved promoter motif, AGCC(G/C)T(G/C)(A/G), upstream of several aflatoxin biosynthetic genes (Roze et al., 2011). Examination of the 1kb upstream regions of the significant DEGs observed in NRRL3357 showed that nine genes possess this predicted atf21 binding domain including a putative cytochrome p450 monooxygenase gene, aflC, and the siderophore biosynthesis acetylase aceI (S2).

Gene Co-expression Is Correlated with Mycelial Dry Weight under Oxidative Stress
In order to identify groups of genes whose expression is directly correlated with isolate oxidative stress tolerance, a co-expression analysis was performed for all isolates and treatments and correlated with previously obtained fungal dry weights (Figure 6). Co-expressed genes were grouped based on expression profile similarities into color-coded modules, which were further separated by hierarchical clustering. Gene network analysis revealed that the largest modules in terms of gene content were the turquoise, blue, and brown modules ( Figure 6A). Module-trait association analysis was then performed to correlate fungal dry weights under each treatment. The blue, purple, and turquoise modules were found to be associated with isolate dry weight under oxidative stress ( Figure 6B). Statistical eigengene expression and fungal dry weight were also found to be significantly correlated with the purple (r = 0.83, p = 0.006) and turquoise (r = 0.73, p = 0.02) module genes showing a strong positive correlation with isolate dry weight and the blue module genes showing a strong negative correlation (r = −0.93, p = 0.0002) with isolate dry weight ( Figure 6C). Gene ontology (GO) enrichment analysis of these modules showed that the blue module was enriched for simple sugar metabolism terms, the purple module for regulation of gene expression and metabolic processes, and the turquoise module for translational mechanisms and additional metabolic processes (Dataset S6). The blue and turquoise modules were also found to contain the aflatoxin biosynthetic genes while the purple module contained several MFS transporter and antioxidant genes (Dataset S6).

DISCUSSION
The production of secondary metabolites in soil dwelling fungi has been shown to be regulated by a number of factors, particularly environmental stress (Calvo and Cary, 2015). Aflatoxin production by A. flavus has been shown to be stimulated by a number of compounds in vitro, with reactive oxygen species (ROS) being of particular interest given their prevalence in host plant defense signaling and environmental stress, particularly drought stress, responses (Jayashree and Subramanyam, 2000;Bhattacharjee, 2005;Roze et al., 2013;Yang et al., 2015Yang et al., , 2016. Such ROS have also been shown to be required for the production of aflatoxin with treatment with antioxidant compounds resulting in reduced aflatoxin production and fungal growth (Nesci et al., 2003;Passone et al., 2005;Mahoney et al., 2010). However, the practical role of aflatoxin production in A. flavus biology is still not well understood. In order to better understand to biological role of aflatoxin production in A. flavus, and its role in oxidative stress responses, we examined the transcriptomes of three aflatoxigenic field isolates exposed to increasing levels of H 2 O 2 -derived oxidative stress.
The isolates which had previously exhibited greater levels of H 2 O 2 tolerance and aflatoxin production , Tox4 and AF13, exhibited a fewer numbers of significantly differentially expressed genes in comparison to the isolate with less tolerance, NRRL3357 ( Table 1). This trend also resulted in strong correlations between previously observed stress tolerance and the numbers of DEGs ( Figure S1). Interestingly, the AF13 isolate exhibits both a similar number and composition of DEGs at 25 mM H 2 O 2 in comparison to NRRL3357 at 10 mM ( Figure S1; Datasets S3, S4). The requirement of additional stress to elicit a similar response is commensurate with the observed stress tolerance and aflatoxin production capabilities of AF13 and NRRL3357 . Tox4, however, may possess stress responsive mechanisms which at more vigorously earlier accounting for the reduced DEG counts observed in this study. Together these results indicate, as has been previously suggested in the literature, that aflatoxin production may contribute to oxidative stress tolerance in Aspergillus spp. (Narasaiah et al., 2006;Roze et al., 2013Roze et al., , 2015Fountain et al., 2014Fountain et al., , 2016. Examining the expression of aflatoxin biosynthetic genes showed that the NRRL3357 isolate exhibited the up-regulation and expression of a greater number of these genes than the other isolates examined in the study (Table 1, Figure 3A). An increase in expression, though to a lesser extent, was also observed for AF13 under higher levels of stress, but not for Tox4 ( Figure 3A). This was surprising given that aflatoxin biosynthetic gene expression is typically higher with periods in fungal development in vitro that coincide with the highest levels of aflatoxin production, approximately 2 to 6 days in conducive media (Davis et al., 1966). The high level of expression of aflatoxin genes in NRRL3357 at 7 days under high levels of stress coupled with the reduced fungal biomass and conidiation observed for this isolate and not with AF13 or Tox4 in our preliminary study , suggests that a delay in fungal growth and development may contribute to prolonged aflatoxin gene expression in addition to stimulation provided by oxidative stress. In addition, these differential expression patterns of aflatoxin biosynthetic genes may also be indicative of post-transcriptional regulation of aflatoxin production and oxidative stress responses at either the RNA or protein level. For example, there has been a low degree of correlation observed between gene expression and proteomics data for traits such Mutual correlations between corresponding modules is higher than most module associations with dry weight, however the purple and turquoise modules exhibit a strong positive correlation with dry weight. In contrast, the blue module shows a strong negative correlation. Eigengene adjacency, which reflects correlation in the heatmap, increases from blue to red in the heatmap.
as temperature stress (r = 0.14) suggesting post-transcriptional regulation of protein accumulation in A. flavus (Bai et al., 2015). Also, given that aflatoxin biosynthetic gene expression can be detected as early as 8 h in culture, there is ample opportunity for such regulation to have a significant effect on later responses to stress and warrants further investigation (Price et al., 2006). However, factors in addition to aflatoxin production have been shown to be at play which influence A. flavus tolerance to oxidative stress. For example, Reverberi et al. (2012) examined the effects of excessive peroxisome proliferation in a mutant isolate expressing the P33 gene. They found that the P33-containing isolate exhibited elevated ROS and oxylipin accumulation commensurate with excess peroxisomes, and countered this with increased aflatoxin production and higher antioxidant enzyme activity compared to the wildtype isolate. The mutant isolate also exhibited similar fungal growth but reduced conidiation compared to the wildtype. Elevated expression of antioxidant and aflatoxin genes at the time point examined in our study would imply that continued stress is still being experienced by NRRL3357, and not to the same extent in AF13 or Tox4. This suggests that NRRL3357 is less able to sequester ROS and remediate its damage than the other isolates. Our previous observation that culturing these isolates in an aflatoxin production non-conducive medium results only in a slight reduction in stress tolerance also suggests that additional oxidative stress remediation mechanisms are at play in all isolates in addition to aflatoxin production .
Two bZIP transcription factor genes, atf21 and atfA, whose homologs in other Aspergillus spp. have been shown to integrate oxidative stress responses and sterigmatocystin/aflatoxin biosynthesis regulation (Sakamoto et al., 2008(Sakamoto et al., , 2009Lara-Rojas et al., 2011;Roze et al., 2011Roze et al., , 2013Hong et al., 2013a,b), also displayed some marginally significant increases in expression in NRRL3357 and AF13 ( Figure S2). Atf21, whose ortholog in A. parasiticus, atfB, has been well characterized (Roze et al., 2011), was found to bind a conserved motif in the promoter regions of several aflatoxin biosynthetic genes. This same motif was found in the upstream promoter regions of nine DEGs in NRRL3357 under high levels of stress, and include aflatoxin production, oxidative stress response, and siderophore biosynthetic genes ( Figure S2). This suggests that these transcription factors may actively regulate gene expression following aflatoxin production initiation (Roze et al., 2011(Roze et al., , 2013 and further supports the case for their involvement in stress tolerance. Also, the co-expression analysis identified several modules of genes exhibiting similar expression profiles which were closely associated with fungal biomass under oxidative stress, including the blue, purple, and turquoise modules ( Figure 6C). These modules contain genes enriched for carbon metabolism and gene expression regulation, but also contain antioxidant enzyme-encoding genes and the aflatoxin biosynthetic genes (Dataset S6). This further supports the observed association of aflatoxin production and oxidative stress tolerance in A. flavus.
The link between aflatoxin production and fungal development has been widely documented in the literature to involve G-protein mediated signaling along with numerous transcription factors and regulatory proteins (Amaike and Keller, 2011). Here, velvetA (VeA), a key developmental regulator influencing sclerotial and conidia formation in response to light levels (Bayram et al., 2008), was constitutively expressed regardless of stress levels while conidiation-related genes such as flbC, a conidiation-specific family protein gene, and a duf221 domain protein gene were also upregulated under stress conditions (Ben-Ami et al., 2010;Kwon et al., 2010). This suggests as has been previously observed Roze et al., 2015) that oxidative stress can stimulate conidiation in Aspergillus spp. Benzoate and amino benzoate degradation genes were also regulated in response to stress (Datasets S3-S5). These compounds have been shown to reduce aflatoxin production and slow fungal growth and development when applied in vitro providing further indication that fungal development may be affected by oxidative stress exposure (Chipley and Uraih, 1980).
Chitin catabolism-related genes tended to be down-regulated under oxidative stress (Datasets S3, S4) which may reduce cell wall degradation and enhance cell well integrity to protect against oxidative stress (Fuchs and Mylonakis, 2009;Jain et al., 2011;Fountain et al., 2016). Also, antibiosis-related genes encoding multiple drug resistance proteins were also upregulated, particularly in AF13, in response to stress. This may provide additional protection for this isolate in competition with other soilborne or plant parasitic microorganisms, especially under environmental stress (Daguerre et al., 2014). NRRL3357 also exhibited an upregulation in cinnamoyl-CoA reductase gene expression in response to stress (Dataset S3). This is noteworthy since it was only fairly recently found that phenylpropanoid metabolic genes were present in the genomes of fungi including the Aspergilli (Seshime et al., 2005). Several members of this class of compounds have been found to influence aflatoxin production by A. flavus in vitro, and have been shown to function as antioxidants which could enhance stress tolerance and competitive capability of the isolates (Chipley and Uraih, 1980). This may indicate that phenylpropanoids or flavonoids produced in host plants along with ROS during drought stress may also be able to influence the production of aflatoxin by A. flavus and warrants further study (Alvarez et al., 2008;Payton et al., 2009). Genes encoding antioxidant enzymes such as thioredoxin peroxidase and thioredoxin reductase were also significantly regulated in response to increasing stress. In the co-expression analysis, the purple module, which exhibited the highest positive correlation with fungal biomass, also contained antioxidant enzyme encoding genes such as glutathione-Stransferase (Dataset S6).
Two additional secondary metabolite biosynthetic gene clusters were also induced by oxidative stress, aflatrem and kojic acid. Aflatrem is a tremerogenic mycotoxin which causes stagger's syndrome in cattle and neurological effects in mammals (Nicholson et al., 2009). Aflatrem is also proposed to function in fungivory prevention for soilborne fungi (TePaske and Gloer, 1992). The biosynthesis of aflatrem is carried out by the indolediterpenoid (isoprenoid) pathway in Aspergillus spp. and shares with aflatoxin biosynthesis the commonality of cytochrome p450 monooxygenase activity (Nicholson et al., 2009;Calvo and Cary, 2015). This along with the activities of various oxidases has led to the proposition that aflatoxin production may function in the alleviation of oxidative stress by promoting antioxidant enzyme production and conidial oxidative stress tolerance through secondary ROS production (Roze et al., 2015;Fountain et al., 2016). We also propose that this process may result in the fixation of excess molecular oxygen (O 2 ) and detoxified ROS into the aflatoxin molecules which are then secreted from the cell as an indirect means of relieving oxidative stress . Alone or in combination, these two proposals would also suggest that isolates which produce higher levels of aflatoxin would exhibit greater oxidative stress tolerance, a trend observed here and in our preliminary study, but which will require further study to confirm . Given the prevalence of cytochrome p450 monooxygenase activity in the aflatrem biosynthetic mechanism, it is possible that a similar phenomenon may occur which is also in need of further investigation.
Kojic acid is a bi-or tri-dentate iron chelating antioxidant compound produced by several species of Aspergilli including A. sojae (Chang et al., 2010). The compound and derivatives of it are useful for the treatment of iron overload in people with thalassemia who accumulate dangerously high levels of iron in their blood stream following multiple transfusions (Nurchi et al., 2016). The iron chelating properties of kojic acid may also, as in the case of cyclopiazonic acid (CPA), may help soil-borne fungi in the remediation of iron starvation under stress conditions by fixing additional iron cations and outcompeting other microbes (Chang et al., 2009). This may also be affected by oxidative stress in A. flavus since siderophore biosynthetic genes were also found to be regulated in AF13 and NRRL3357. Also, the chelation of excess iron cations, particularly Fe 2+ , may prevent Fenton reaction-derived ROS such as hydroxyl radicals (OH .− ) from forming and resulting in additional cellular damage (Fenton, 1894). Therefore, the production of kojic acid may also supplement A. flavus oxidative stress tolerance. Also, these same hydroxyl radicals have also been shown to influence aflatoxin production with treatment of A. flavus mycelia with DMSO, a hydroxyl radical scavenger, resulting in significantly reduced aflatoxin production and sclerotial differentiation (Grintzalis et al., 2014).
Together, these results provide an insight into a potentially coordinated secondary metabolite response to oxidative stress in A. flavus. While the biological function of the aflatoxin compound in fungal biology remains unclear, it seems possible that the aflatoxin production process may contribute to oxidative stress tolerance. However, in the absence of aflatoxin production capability, the production of other secondary metabolites may compensate for any reduction in stress tolerance, and may be of interest in the biology of atoxigenic biological control isolates. This role of oxidative stress in stimulating aflatoxin production in A. flavus is made all the more interesting by the observed accumulation of ROS in drought stressed host plant tissues , and the correlation of elevated antioxidant enzyme and compound accumulation in hosts resistant to aflatoxin contamination (Chitarrini et al., 2014). Together, this may be a contributing factor to the observed correlation between drought stress susceptibility and aflatoxin contamination in host plant species such as maize and peanut (Holbrook et al., 2009;Kebede et al., 2012). The role of oxidative stress in both secondary metabolite production and fungal primary metabolism may also provide a basis for identifying host defense attributes which contribute to reduced aflatoxin contamination in the form of selectable biomarkers and in the selection and characterization of novel resistance genes and genetic markers identified through quantitative trait loci (QTL) studies for marker assisted selection (MAS) and cultivar resistance improvement.

AUTHOR CONTRIBUTIONS
JF: Performed the experiments, assisted in data analysis, and wrote the manuscript; PB: Performed the data analysis; MP and SN: Assisted with manuscript preparation and data analysis; VK: Performed the sequencing; AJ and AC: Assisted in RNA preparation for sequencing. LY, RL, and RK: Contributed to the design and project discussion; RV and BG: Conceived and supervised the project, secured funding, and revised and submitted the manuscript.

ACKNOWLEDGMENTS
We thank Billy Wilson, Hui Wang, Xiaohong Guo, Xiangyun Ji, Gaurav Agarwal, and Hui Song for technical assistance in the laboratory. This work is partially supported by the U.S. Department of Agriculture Agricultural Research Service (USDA-ARS), the Georgia Agricultural Commodity Commission for Corn, the Georgia Peanut Commission, the Peanut Foundation, and AMCOE (Aflatoxin Mitigation Center of Excellence, Chesterfield, MO, USA). This work has also been undertaken as part of the CGIAR Research Program on Grain Legumes and the USAID University Linkages Program between USDA-ARS and ICRISAT. ICRISAT is a member of CGIAR Consortium. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. The USDA is an equal opportunity provider and employer.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2016.02048/full#supplementary-material Figure S1 | Linear regression analysis of the numbers of differentially expressed genes relative to isolate H 2 O 2 tolerance. The numbers of significantly, differentially expressed genes were plotted relative to the observed maximum H 2 O 2 level tolerated by the isolates observed in Fountain et al. (2015). Significant correlations were observed when comparing 0 and 10 mM H 2 O 2 (R 2 = 0.9580), and 10 and 20/2 5 mM H 2 O 2 (R 2 = 0.9174). The correlation between 0 and 20/25 mM H 2 O 2 was not as significant (R 2 = 0.4332) possibly due to a similarity between NRRL3357 and AF13 responses to high levels of stress. Numerical data can be found in Table 1. Figure S2 | FPKMs of genes encoding the bZIP transcription factors atf21 and aftA. The expression of the bZIP transcription factor genes aft21 (A) and aftA (B) was marginally, significantly regulated in response to high levels of stress in NRRL3357 and AF13. While in the study, a p-value and a q-value ≤ 0.05 with a log 2 (fold change) ≥ 2 were considered to be significant, these genes were numerically and marginally expressed higher in response to stress. Green bars in the plot represent 0 mM H 2 O 2 , blue bars 10 mM H 2 O 2 , and yellow bars 20/25 mM H 2 O 2 . The bracket displays the p-and a q-values observed between the 0 and 20/25 mM H 2 O 2 treatments.
Promoter analysis with MEME identified nine genes (C) in NRRL3357 which were differentially expressed between 0 and 20 mM H 2 O 2 that contain the conserved motif (D) recognized by the atf21 ortholog, atfB, in A. parasiticus (Roze et al., 2011).

Additional Information
Data availability: Raw and processed RNA sequencing data can be accessed at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) with accession number SRP094082 under BioProject PRJNA348383.