Carbon Dioxide Mediates the Response to Temperature and Water Activity Levels in Aspergillus flavus during Infection of Maize Kernels.

Aspergillus flavus is a saprophytic fungus that may colonize several important crops, including cotton, maize, peanuts and tree nuts. Concomitant with A. flavus colonization is its potential to secrete mycotoxins, of which the most prominent is aflatoxin. Temperature, water activity (aw) and carbon dioxide (CO₂) are three environmental factors shown to influence the fungus-plant interaction, which are predicted to undergo significant changes in the next century. In this study, we used RNA sequencing to better understand the transcriptomic response of the fungus to aw, temperature, and elevated CO₂ levels. We demonstrate that aflatoxin (AFB₁) production on maize grain was altered by water availability, temperature and CO₂. RNA-Sequencing data indicated that several genes, and in particular those involved in the biosynthesis of secondary metabolites, exhibit different responses to water availability or temperature stress depending on the atmospheric CO₂ content. Other gene categories affected by CO₂ levels alone (350 ppm vs. 1000 ppm at 30 °C/0.99 aw), included amino acid metabolism and folate biosynthesis. Finally, we identified two gene networks significantly influenced by changes in CO₂ levels that contain several genes related to cellular replication and transcription. These results demonstrate that changes in atmospheric CO₂ under climate change scenarios greatly influences the response of A. flavus to water and temperature when colonizing maize grain.

1. Introduction

Aspergillus flavus is a saprophytic fungus that infects several crops of agronomic importance, including corn, cotton, peanuts and tree nuts. Prior to harvest, A. flavus may infect the fruiting bodies or seeds in crops and produce several toxic secondary metabolites, including the polyketide derived aflatoxins (AFs), cyclopiazonic acid and aflatrem. In the USA and other industrialized countries, the establishment of contamination thresholds by regulatory agencies and close monitoring of crops have minimized the direct impact on human health. However, economic losses remain significant. Regarding AF contamination in the USA alone, estimates of losses are between $163 million for maize crops to $500 million annually for maize, peanuts and other crops [1,2]. In lower middle income countries (LMCs) where the regulatory controls either do not exist or are not enforced, especially in sub-Saharan Africa, the consumption of food contaminated with AFs are directly linked to liver disease, tumor development, stunted development in children and other medical defects (reviewed in [3]). The AF compounds produced by A. flavus include the structurally similar forms B1, B2, G1 and G2. Most A. flavus strains produce only the B forms of AF, however other related species, such as A. parasiticus and A. nomius, produce both B and G AFs [4,5].

Mitigation strategies to minimize contamination of maize have involved the implementation of different approaches. These include the development of genetically modified (GM) maize crops [6,7], which are resistant to pests and minimize damage that provides entry points for A. flavus; classic breeding for resistance [8], and pre-harvest biocontrol by using naturally non-aflatoxigenic strains of A. flavus to outcompete the toxigenic strains [9]. The latter approach has resulted in several commercial biocontrol products currently being used in maize, groundnuts and cotton.

Environmental factors, including water activity (aw), temperature, light, as well as their interactions, have been demonstrated to have a significant influence on germination, growth and AFs production by strains of A. flavus [10,11,12,13]. Recent studies by Medina et al. (2015) demonstrated that water activity had a statistically significant effect on growth rates, whereas the effect of temperature and CO2 was negligible under the conditions tested. Further studies by Medina et al. (2017) examined the impact of aw and temperature interactions on global transcriptomic changes in A. flavus and how this relates to AFB1 production in maize-based matrices.

There has been interest in the resilience of mycotoxigenic fungi in relation to potential future climate change (CC) scenarios and whether interactions between aw,increased temperature andelevated CO2 might impact mycotoxin production at a fundamental genomic level and phenotypic toxin production. By the year 2100 temperatures could increase by 4 °C, and CO2 levels are anticipated to reach approximately 1000 ppm, depending on mitigation efforts [14]. Medina et al. (2015) showed that when A. flavus was exposed in vitro to these three-way interacting CC factors growth was unaffected, while expression of biosynthetic genes (aflD; aflR) and phenotypic AFB1 production were significantly increased. This is likely to result in economic impacts. For example, it has been suggested that CC could increase economic losses in the maize industry by $50 million to up to $1.7 billion if AF production increases under such scenarios [15]. A recent review suggested that food security, especially in LMCs, could be profoundly impacted on due to such increases in mycotoxin contamination of staple food commodities [16].

There has been interest in applying current technologies, including functional genomics, transcriptomics, and proteomics to measure the impact of interacting CC-related abiotic parameters on fungal growth, the regulation of the AF biosynthetic cluster, and toxin production. A recent study examined the impact of interactions between aw and temperature in the transcriptome of A. flavus when colonizing maize [17]. It was determined that aw (0.99, 0.91) and temperature (30, 37 °C) influenced the colonization of maize grain by A. flavus with a significant effect on AFB1 production at 0.91 aw and 37 °C. Both environmental factors affected biological processes and the numbers of up- and downregulated genes. The interacting environmental factors influenced the functioning of the secondary metabolite clusters for AFs and CPA such that an elevated number of genes were co-regulated by both aw and temperature. For example, an interaction effect for 4 of the 25 AFB1 genes, including regulatory and transcription activators occurred. For CPA, all 5 biosynthetic genes were affected by aw stress, regardless of temperature.

The objectives of the present study were to examine, for the first time, the impact of interactions between aw (0.99, 0.91), temperature (30, 37 °C) and CO2 (350, 650 and 1000 ppm) on A. flavus (NRRL 3357) colonization of stored maize grain to evaluate the impacts on (1) the whole genome, including the aflatoxin biosynthesis gene cluster using RNA-seq, (2) effects on key biosynthetic genes using q-PCR, and (3) effects on AFB1 production. We used functional genomics to better understand the transcriptomic response of the fungus to CC parameters of aw, temperature and CO2 levels, with a view towards predicting changes in fungal infection and toxin production associated with resilience to such climatic stress factors.

2. Results

2.1. The Interaction of Water, Temperature, and CO2 Impact AFB1 Production in Maize Grain

Overall, AFB1 production showed a positive correlation with water availability and CO2 levels. AFB1 was extracted and quantified after 10 days incubation at 30 °C or 37 °C, 0.91 or 0.99 aw, and 350, 650, or 1000 ppm CO2. At low aw (Figure 1A) the quantity of AF produced was 5- to 10-fold lower than at high aw levels (Figure 1B). However, the effect of CO2 was positively correlated with AFB1 production at both high and low aw levels. With an elevated CO2 level of 650 ppm, there was an interaction effect with temperature and aw. At 0.91 aw, AFB1 production was highest at 37 °C. Conversely, at 0.99 aw the temperature condition of 30 °C exhibited higher toxin levels. The trends observed indicate that CO2 affects toxin production. The maximum quantity of toxin was observed at elevated levels of CO2 (650 ppm) with low water activity and at 30 °C, however, no additional accumulation of AFB1 was observed at 1000 ppm.

Figure 1

Aflatoxin B1 (AFB1) production by A. flavus under different combinations of environmental conditions. (A) At 30 °C and under low water activity levels of 0.91 aw (“30/91”), the effects of CO2 are minimal; however, at 37 °C (“37/91”) there are increases in AFB1 production correlating with higher CO2 levels. (B) At high water activity levels the biosynthesis of AFB1 was significantly elevated at both 30 °C (“30/99”) and 37 °C (“37/99”). ANOVA: analysis of variance statistical analysis. * p ≤ 0.05.

2.2. Effect of Three-Way Interacting CC Conditions on Gene Expression

RNA-seq of A. flavus genes showed a large global effect with water and temperature stress, but a limited effect with elevated CO2 levels. Between 6.22 × 105 and 3.37 × 107 reads mapped to exogenic regions of A. flavus strain 3357 (Table 1). Principle component analysis (PCA) indicated that aw caused the largest variance observed (Figure 2A), accounting for 67% of the variance, however samples with similar water activity cluster together on PC1, indicating there little sample variation for this metric. The second principle component indicated that a change in temperature was the second most important factor, but only for samples stored with freely available water (0.99 aw). All samples stored under water stress (0.91 aw) clustered together and showed relatively little inter-sample variance. At high aw, a change in temperature resulted in a high variance in gene expression. The total number of genes affected under the various conditions are illustrated in Figure 2B. For each data point in Figure 2B, the unidentified variable (aw level, temperature or CO2 level) is the baseline condition (30 °C, 0.99 aw, or 350 ppm CO2) for those samples being compared. Of note, ~10-fold more genes responded to temperature and CO2 stress at high aw levels (Figure 2B, upper graph). Likewise, ~3x more genes were affected by temperature and ~25x more genes were affected by CO2 in the non-stressed temperature conditions (30 °C; Figure 2B, middle graph). High CO2 levels (1000 ppm) had a significant impact on the number of genes expressed, decreasing the number of genes affected by water and temperature by ~3 fold (Figure 2B bottom graph). In summary, elevated temperature, water stress, and increased CO2 decreased the number of differentially expressed genes when compared to the non-stressed conditions.

Table 1

Number of reads mapping to exons for all conditions tested.

		30 °C	37 °C
0.91 aw	350 ppm	3.37 × 107	2.13 × 107
	650 ppm	2.64 × 107	3.00 × 106
	1000 ppm	3.04 × 107	6.84 × 106
0.99 aw	350 ppm	1.04 × 107	8.01 × 105
	650 ppm	4.69 × 106	2.24 × 106
	1000 ppm	6.22 × 105	1.03 × 106

Figure 2

(A) Principle component analysis and (B) total gene counts of differentially expressed genes indicate that of the three environmental variables tested, water activity is the primary driver of transcriptional changes, with temperature also having a measurable impact. Carbon dioxide levels (1000 ppm vs. 350 ppm) has the largest impact at 0.99 aw (B, top) and at 30 °C (B, middle). High carbon dioxide levels (1000 ppm) reduces the level of genes affected by water activity and temperature by approximately 1/3 (B, bottom).

A significantly larger number of genes were affected by water stress, whereas approximately equal numbers were affected by temperature and CO2 (Figure 3, left and center). Most of the genes that were affected by changing temperature and CO2 levels were also affected by low aw, with 472 genes being upregulated by all three environmental conditions and 564 genes being downregulated. Figure 3 illustrates how changing all three environmental variables simultaneously resulted in a significantly larger number of genes (4853) being affected than changing variables independently. In other words, there is a cumulative effect of changing environmental conditions on global gene expression. The group labeled “Genes affected individually” refers to genes that are affected by just aw, temperature, or CO2 levels when examined separately. Sixty-nine genes were affected by each of the three conditions individually, but not affected when all three conditions are applied simultaneously.

Figure 3

Venn Diagrams illustrate the number of genes up regulated (left) and down regulated (middle) by individual conditions. For each comparison described, control conditions are assumed (30 °C, 0.99 aw and 350 ppm CO2). When all three conditions are changed simultaneously, 5820 genes are affected, however when only one environmental condition is changed, only 967 of these genes are affected, indicating a significant difference between the cumulative and individual effects of environmental changes (right).

2.3. Effect of Three-Way Interacting CC Conditions on Biological Processes

After 10 days A. flavus colonization of the maize kernels, quantitative PCR was conducted on two AF cluster genes: aflR, the regulating transcription factor, and aflD, a reductase (Figure 4A). Gene expression levels are shown relative to a control condition of 30 °C, 0.99 aw, and 350 ppm CO2. Four important observations were made: (1) at 30 °C/0.99 aw there were decreases in aflR and aflD expression, even at elevated CO2 levels (650 ppm; 1000 ppm) relative to the control (350 ppm). This is in contrast to what was observed at 37 °C, where aflD is activated at higher levels in the 650 and 1000 ppm CO2 treatments at 0.99 aw; (2) in two separate stressed conditions 30 °C/0.91 aw (water stressed conditions) and 37 °C/0.99 aw (high temperature), there are sharp differences in aflR expression levels between 650 and 1000 ppm, indicating that some CO2 level between the two may represent a threshold; (3) at low aw and high temperature, there were higher levels of activity than at 0.99 aw/350 ppm; (4) at 37 °C/0.99 aw, a dose response in relation to CO2 was observed.

Figure 4

Gene Expression of the aflatoxin gene cluster. (A) Quantitative PCR analysis shows the expression of aflR and aflD, an aflatoxin cluster transcription factor and structural gene, respectively. After 10 days of incubation on maize kernels gene levels at 30 °C generally decrease, however the effects of high CO2 (1000 ppm) levels at 0.91 aw indicate decrease values (left), possibly in response to elevated AFB1 levels. At 37 °C gene levels remain high, however, again at 1000 ppm CO2 the transcription factor aflR is decreased. (B) The heat map of regularized log transformed counts indicate hierarchal clustering associated with water activity levels. The clustering also indicates expression patterns that suggest early genes in the pathway may be responsive to high CO2 levels (See Results).

The heat map (Figure 4B) indicates relative expression levels for all genes in the AF gene cluster obtained from RNA-seq. Hierarchical clustering of sample conditions (top brackets) and genes (left side) are indicated. The AF cluster genes aflD and aflR are clustered relatively close together (genes 3 and 6, respectively), indicating similar expression profiles, which is in agreement with qPCR results where 9 of the 12 samples are similarly expressed. Hierarchical clustering of the samples indicated water availability as the primary determinant of expression patterns. After 10 days, the 30 °C/0.99 aw/350 ppm CO2 condition exhibited higher overall expression levels, distinct from any of the stressed conditions. Notably, expression at 30 °C/0.99 aw/350 ppm is higher than at 650 ppm and 1000 ppm CO2. Three of the genes in the AF biosynthesis pathway, fatty acid synthases aflA and aflB, and the polyketide synthase aflC, all of which are responsible for the synthesis of norsolorinic acid, were enriched in all samples. The reductases coded by aflD and aflF, responsible for the reduction of norsolorinic acid, were also enriched at relatively high levels; however, a third reductase, aflE, exhibited the highest enrichment level in the most stressed condition examined (37 °C/0.91 aw and 1000 ppm CO2).

Enrichment analysis was conducted to identify KEGG biological processes with an over-representation of genes being up- or downregulated. Table 2 lists the KEGG categories affected by changes in conditions, both individually and combined, from which several important observations can be made. Primarily, energy-related metabolic processes such as glycolysis/gluconeogenesis, purine/pyrimidine metabolism, and starch/sucrose metabolism are prominent among conditions analyzed. Most genes in this category that are essential for glycolysis, including triosephosphate isomerase (AFLA_094630), fructose bis phosphate aldolase (AFLA_030930), and pyruvate kinase (AFLA_087900) are all downregulated in all the conditions examined. The KEGG categories affected, that are unique to changes in CO2 concentration (30 °C/0.99 aw/350 ppm CO2 vs. 30 °C/0.99 aw/1000 ppm CO2), are cysteine and methionine metabolism (p = 0.009) and folate biosynthesis (p = 0.030). Three categories of KEGG were affected only when all three environmental conditions were changed simultaneously (37 °C/0.91 aw/1000 ppm CO2 vs. 30 °C/0.99 aw/350 ppm CO2). These were taurine and hypotaurine metabolism (p = 0.014), cyanoamino acid metabolism (p = 0.017), and ether lipid metabolism (p = 0.049).

Table 2

KEGG Categories with a statistically significant overrepresentation of genes that are differentially expressed and their associated p values.

Carbon Dioxide		Temperature
(30 °C/0.99 aw/350 ppm vs. 30 °C/0.99 aw/1000 ppm)		(37 °C/0.99 aw/350 ppm vs. 30 °C/0.99 aw/350 ppm)
KEGG Category	p-Value	KEGG Category	p-Value
Glycolysis/Gluconeogenesis	0.003	Glycolysis/Gluconeogenesis	0.000
Purine metabolism	0.008	Starch and sucrose metabolism	0.002
Cysteine and methionine metabolism	0.009	Methane metabolism	0.017
Fructose and mannose metabolism	0.011	Riboflavin metabolism	0.024
Pyrimidine metabolism	0.025	Glutathione metabolism	0.024
Folate biosynthesis	0.030	Glycosphingolipid biosynthesis - globo series	0.025
Carbon fixation in photosynthetic organisms	0.041	Fructose and mannose metabolism	0.029
Inositol phosphate metabolism	0.049	Pentose phosphate pathway	0.043
Water Activity		Combined
(30 °C/0.91 aw/350 ppm vs. 30 °C/0.99 aw/350 ppm)		(37 °C/0.91 aw/1000 ppm vs. 30 °C/0.99 aw/350 ppm)
KEGG Category	p-Value	KEGG Category	p-Value
Riboflavin metabolism	0.001	Glycolysis/Gluconeogenesis	0.000
Inositol phosphate metabolism	0.002	Methane metabolism	0.002
Glyoxylate and dicarboxylate metabolism	0.011	Riboflavin metabolism	0.005
Purine metabolism	0.012	Inositol phosphate metabolism	0.009
Methane metabolism	0.019	Taurine and hypotaurine metabolism	0.014
Starch and sucrose metabolism	0.020	Fructose and mannose metabolism	0.014
Glycolysis/Gluconeogenesis	0.023	Starch and sucrose metabolism	0.015
Cyanoamino acid metabolism	0.031	Cyanoamino acid metabolism	0.017
Fructose and mannose metabolism	0.039	Purine metabolism	0.019
		Pyrimidine metabolism	0.022
		Ether lipid metabolism	0.049

2.4. Effect of CO2 and Interactions with Other Abiotic Factors on Secondary Metabolite Gene Clusters

To identify secondary metabolic gene clusters that were affected by the interacting environmental conditions tested according to the RNA-seq results, a Secondary Metabolite Unique Regions Finder (SMURF) analysis of the A. flavus genome was conducted. Table 3 lists several SMURF-identified secondary metabolic gene clusters labeled according to Georgianna et al. [18] and their corresponding relative expression values (log2 fold change). The full list of secondary metabolite-associated genes identified by SMURF is provided in Supplementary Table S1.

Table 3

Known or putative gene clusters in Aspergillus flavus identified by the gene for their primary backbone enzyme. The fold change values indicate upregulation (positive number) or downregulation (negative number) according to RNA sequencing results. Clusters in bold indicate the gene is affected by both carbon dioxide levels tested.

				Effect of Carbon Dioxide		Effect of Water		Effect of Temperature
				At 30 °C/0.99 aw:		At 30 °C, 0.91 aw vs. 0.99 aw:		At 0.99 aw, 37 °C vs. 30 °C:
#	A. flavus 3357	Name	SM* Product	650 ppm/350 ppm	1000 ppm/350 ppm	350 ppm	1000 ppm	350 ppm	1000 ppm
5	AFLA_006170	polyketide synthetase (PksP)	naphthopyrone	-	3.37	5.01	1.96	1.86	-
10	AFLA_016140	scytalone dehydratase (Arp1) (conidial pigment biosynthesis)	conidial pigment 1,8-dihydroxynaphthalene-melanin	-	-	−3.54	−2.21	-	-
15	AFLA_045490	dimethylallyl tryptophan synthase, putative	aflatrem, ATM2	2.50	-	-	-	-	-
19	AFLA_060680	dimethylallyl tryptophan synthase	Unknown	−3.38	−4.09	−4.17	-	-	2.68
20	AFLA_062860	polyketide synthase (PkfA)	3-(2,4-dihydroxy-6-methylbenzyl)-O rsellinaldehyde	1.35	-	-	−1.87	2.42	-
21	AFLA_064240	nonribosomal peptide synthase (wykN)	WYK peptidase inhibitor	-	2.23	-	−2.14	1.31	−2.19
23	AFLA_066840	hybrid NRPS/PKS enzyme	Leporins	1.59	2.45	-	−3.27	-	−3.40
35	AFLA_101700	NRPS enzyme (lnaA)	piperazines	-	2.15	-	−1.94	2.66	-
36	AFLA_104210	PKS-like enzyme, putative	dihydrocurvularin	-	-	−3.04	-	−2.94	-
39	AFLA_108550	polyketide synthase	monodictylphenone	-	-	-	-	-	-
41	AFLA_114820	polyketide synthase (fluP) (pksL2)	6-MSAi	-	1.55	-	−2.17	-	-
44	AFLA_116890	polyketide synthase (PkiA)	6-hydroxy-7-methyl-3-nonylisoquinoline-5,8-dione	-	−2.52	−3.49	-	-	-
54	AFLA_139410	polyketide synthase (aflC/pksA/pksL1)	aflatoxin	-	1.20	0.95	-	-	−2.05
55	AFLA_139490	hybrid PKS/NRPS enzyme	cyclopiazonic acid	-	-	−2.75	−5.39	-	−3.09
	AFLA_125760	Class 2 Terpene Cyclase	Unknown	1.99	4.65	3.98	-	3.48	-

* Putative product based on identity to characterized enzymes.

Many of the genes listed in Table 3 have high sequence identity to previously characterized genes (in black). Of note, four of the genes (shown in red) are affected by both of the elevated CO2 levels (650 ppm and 1000 ppm). These consist of two dimethylallyl tryptophan synthases (DMTS), a terpene cyclase, and the hybrid NRPS/PKS shown to be responsible for leporin biosynthesis (see Discussion).

2.5. Identification of Gene Networks Affected by CO2 Levels

Gene co-expression networks were determined by analyzing all RNA sequencing results using WGNCA. The results were then visualized using Cytoscape. This analysis revealed two prominent gene networks with significant numbers of genes affected by CO2 levels (Figure 5).

Figure 5

Weighted Gene Network Co-expression Analysis followed by visualization in Cytoscape shows two unidentified networks heavily influenced by increased carbon dioxide levels (1000 ppm CO2 vs. 350 ppm CO2) at 30 °C and 0.99 aw. The color of the node indicates the log2 fold change values. (A) 268 out of 905 genes in the network are differentially expressed, with most showing increased levels. (B) 114 genes out of 415 genes in the network show altered expression levels. See Results for description of the genes in the network.

Genes with altered differential expression at 30 °C, 0.99 aw and 1000 ppm CO2, compared with 350 ppm, are colored and genes not differentially expressed are shown in grey. The network illustrated in Figure 5A shows 905 genes interacting in total, with 268 genes differentially expressed. The network has three distinct clusters on the left, right, and center, and the center cluster is co-expressing with genes in the outer clusters (indicated by black lines (edges)). The left and right clusters (95 and 67 genes, respectively) have edges connecting to the 45 genes in the center cluster. Gene ontology (GO) enrichment analysis of the individual clusters indicated the molecular function of genes on the left involve primarily ATP-binding protein kinases (e.g., GO:0005524, ATP binding; GO:0004672, protein kinase activity) (See Supplemental Table S2). The right cluster is enriched in structural and secretory-related genes (e.g., GO:0000166, nucleotide binding, Rheb, and Myo5; GO:0051056 regulation of small GTPase mediated signal transduction, Sar1). The center of the network contained no enriched GO categories, however it does contain several transcription factors and transcriptional regulatory elements (AFLA_029620, AFLA_114920, AFLA_003630). Figure 5B is a second network (network 2) identified by WGNCA analysis and Cytoscape rendering. It shows 415 genes in the network, 114 of which are significantly affected by changing CO2 levels. Network 2 consists of a significantly enriched number of ribosomal-related genes (GO:0003735, structural constituent of ribosome). Other non-ribosomal genes include and Hsp90-binding chaperone sba1 (AFLA_095590), a mycelial catalase cat1 (AFLA_090690), and the AF cluster gene aflT (AFLA_139420) (Supplementary Table S2).

3. Discussion

To date, the majority of research pertaining to environmental effects on AF production, fungal growth, and plant pathogenicity have focused on the effects of water availability and temperature. Work by Medina et al. [19,20] was among the first attempts to examine the effects of aw and temperature in the context of higher CO2 levels. They found that the interactive effect of water stress, high temperature and high CO2 increased both gene expression of the AF biosynthesis genes aflR and aflD, and AFB1 production in maize kernels. While the experimental procedure here permitted us to observe changes in gene expression and AF production after 10 days of growth, it is limited in that we cannot observe a temporal response of gene activation followed immediately by AF biosynthesis. However, we further expand on previous findings, allowing us to characterize putative secondary metabolic gene clusters, important developmental genes, and identify networks of co-expressed genes affected by CO2 levels.

Most evidence involving fungal carbon metabolism involves the release of CO2 via cellular respiration, however CO2 serves other functions. Hall et al. [21] reported that CO2 serves as an intra-colony signaling molecule important for colonization and pathogenesis in the fungus Candida albicans. It has also been shown that insects, in symbiosis with fungi, demonstrate preferences for certain ranges of elevated CO2 in which to conduct their fungal-rearing [22]. Furthermore, elevated CO2 present in the ambient environment of soil samples containing a mixed fungal population decreases respiration activity [23] and decreases AF accumulation in A. parasiticus [24]. The data here demonstrates that while marked effects in transcription (and by consequence growth and toxin production) are primarily a result of water availability, and secondarily temperature, changing CO2 levels altered fungal response to both water and temperature changes. This is made apparent by the total numbers of genes affected and by our observing that changes brought on by low water or high temperature stresses can vary depending on CO2 availability. Furthermore, evidence indicates that the combination of environmental stressors, including high CO2, have a compounding effect compared to that of individual stressors.

The production of secondary metabolites under different environmental conditions is of concern due to potential shifts in toxigenic potential. The biosynthetic gene cluster 19 in A. flavus contains a dimethylallyl tryptophan synthase (DMAT), which prenylates tryptophan or other aromatic substrates, however the secondary metabolite biosynthesized has yet to be identified. While the expression profile of individual cluster genes does not necessarily correlate with metabolite biosynthesis, DMATS genes are involved in the biosynthesis of several toxic metabolites produced by fungi that infect and contaminate crops such as cyclopiazonic acid [25] in A. flavus and ergotamine in Claviceps purpurea [26]. Toxic metabolites can serve offensive or defensive functions, increasing virulence of the fungus toward its host, or protecting the fungus from threats by other microorganisms or predation. The backbone gene in cluster 23, a hybrid Non-ribosomal Peptide Synthase/Polyketide Synthase (NRPS/PKS), is upregulated with increasing CO2 levels. This cluster is responsible for the biosynthesis of leporins [27]. Leporin A has anti-insectan properties [28], while leporin B [29] can chelate iron, forming a trimer with Fe3+ [27]. Siderophore biosynthesis and siderophore-mediated iron uptake have been found to increase under hypoxic conditions in A. fumigatus [30]. AFLA_125760, which was also upregulated at elevated CO2 levels, encodes for a terpene cyclase flanked by a putative P450 alkane hydroxylase (AFLA_125750) and a steroid alpha reductase (AFLA_125740). These types of cyclases convert the linear triterpene squalene into cyclized products under hypoxic conditions. Cyclized triterpenes, such as fungal ergosterol and bacterial hopenoids, provide cell membrane structural integrity and fluidity. Production of ergosterol, the most prevalent cyclized triterpene found in fungal cell membranes, is affected by mechanical and oxidative stress [30,31]. Other cyclic triterpenes may also be expressed under hypoxic and stressful conditions to improve cell membrane integrity.

Here we have demonstrated for the first time transcriptome-wide changes that occur in a pathogenic fungus while colonizing maize kernels under these interacting CC-related environmental conditions. AF biosynthesis and genes from the AF biosynthetic cluster responded to elevated CO2 levels, as did several other identified secondary-metabolic gene clusters. Further, there is a global change in the transcriptome in response to water and temperature stress under high CO2 conditions. This works lays a solid foundation for further research to establish which genes and gene networks should be targeted for fungal inhibition in future CC scenarios.

4. Materials and Methods

4.1. Fungal Strain

Aspergillus flavus strain NRRL 3357 (ATCC 200026, GenBank assembly accession: GCA_000006275.2) [32,33] was obtained from the Southern Regional Research Center, New Orleans, LA. This strain has been characterized extensively in several areas, including metabolite production [34,35], genome profiling, [36], for the development of biological control [37], molecular ecology [38] and many others. Spore stocks were stored at 4 °C or sub-cultured on Malt Extract agar (MEA; CM59, Oxoid Ltd., Basingstoke, UK) when needed.

4.2. Sample Preparation and Treatment

Undamaged French feed maize kernels were used in this study. Initially a standard curve was conducted whereby known amounts of water were added to multiple 10 g samples of kernels, incubated at 4 °C for 48 h, and the resulting aw of the kernels was determined using an Aqualab 4 TE water activity meter (Decagon Devices, Pullman, WA, USA).

Based on aw standard curve results, a volume of water was added to each 10 g sample of maize kernels to obtain aw levels of 0.99 and 0.91 (minus 200 μL for the later additions of A. flavus spore suspension). The maize kernels were then placed in glass culture vessels containing a microporous lid, which allows for moisture and air exchange (Magenta, Sigma Ltd., Castleford, UK). Subsequently, 200 μL of spore suspension (approx. 106 spores/mL) were added to make up the predetermined amounts of water required and thoroughly mixed. The inoculated vessels for each treatment were placed in enclosed environmental chambers. To control humidity levels, 2 glass jars (500 mL) containing glycerol-water solutions, appropriate to maintaining the equilibrium relative humidity at the target aw level, were placed in each chamber. The chambers were incubated at 30 and 37 °C for 10 days. To control CO2 and maintain humidity, specialty certified CO2 gas cylinders (British Oxygen Company, Jierfude, UK) containing either air, 650 ppm CO2 or 1000 ppm CO2 were used for flushing the environmental chambers. The gases were bubbled through a glycerol:water solution of the required aw level before flushing through the chambers and the valves closed. CO2 flushing was performed every day as described previously [13]. The glycerol-water solutions in the chambers were replaced with fresh solutions every 2 days during the incubation period. Three replicates per treatment were used in all cases. At the end of the incubation period, samples were snap frozen using liquid N2 and kept at −80 °C until a portion of each sample could be used for RNA extraction and purification, or dried for AFB1 extraction and clean-up prior to quantification using HPLC analysis.

4.3. Aflatoxin Analysis

AFB1 extraction was performed using AflaStar™®—Immunoaffinity Columns (IAC, Romer Labs Inc., St. Louis, MO, USA), following the manufacturer’s instructions. Briefly, 5 g of the sample were dried overnight at 80 °C and stored at room temperature. The samples were ground and 4 g was placed into a 50 mL Falcon tube, to which 16 mL of a methanol:water (60:40 v:v) solution was added. The samples were shaken for 1 h at room temperature, and then filtered through qualitative filter paper (QL 110, Fisher Scientific UK Ltd., Loughborough, UK). The extract (1 mL) was diluted in a 15 mL Falcon tube with 9 mL of PBS buffer (0.05 M/0.15 M NaCl, pH 7.4, Fisher Bioreagents®, Fisher Scientific UK Ltd., Loughborough, UK), and pH was checked with pH strips. The diluted extract was applied to the IAC, and allowed to drip through. After further cleaning, 3 mL of Methanol (HPLC grade) was used to elute the AFs. The eluent was dried and standards were prepared using 200 μL AF (R-Biopharm Rhône Ltd., Darmstadt, Germany) stock solution comprised of 1 ng/uL AFB1. The stock solution was pipetted into 2 mL Eppendorf tubes and left to evaporate to dryness overnight inside a fume hood. For quantification of AFs, 200 μL hexane was added to the residue followed by the addition of 50 μL triflouroacetic acid (TFA). The mixture was then vortexed for 30 s and then left for 5 min. Thereafter, a mixture of water:acetonitrile (9:1, v:v) was added and the entire contents of the tube were vortexed for 30 s, after which the mixture was left for 10 min to allow for thorough separation of layers. The hexane layer was discarded and the aqueous layer filtered through nylon syringe filters (13 mm × 0.22 μm; Jaytee Biosciences Ltd., Herne Bay, UK) directly into amber salinized 2 mL HPLC.

A reversed-phase Agilent 1200 series HPLC system with fluorescence detection was used to confirm the identity and quantify AFB1. This consisted of an in-line degasser, auto sampler, binary pump and a fluorescence detector (excitation and emission wavelength of 360 and 440 nm, respectively). Separation was achieved through the use of a C18 column (Agilent Zorbax Eclipse plus C18 4.6 mm × 150 mm, 3.5 μm particle size; Agilent, Berks, UK) preceded by a guard cartridge with the same packing material. Isocratic elution, with a mobile phase that included methanol:water:acetonitrile (30:60:10, v:v:v), was performed at a flow rate of 1.0 mL/min. The injection volume was 20 μL. A set of standards was injected (1 to 5 ng AFB1, B2, G1 and G2 per injection) and standard curves were generated by plotting the area underneath the peaks against the amounts of AFB1 standard injected. Linear regression was performed in order to establish a correlation relationship (correlation coefficient, R2 = 0.99).

4.4. Total RNA Extraction

For RNA-seq, tissue harvest was performed after 10 days using three replicates. This time frame was chosen because previous studies with both A. flavus and A. parasiticus suggested that gene expression of many of the biosynthetic genes was optimal after 8–10 days growth, although there does appear to be a sequential expression of groups of AF biosynthetic genes [39,40,41]. Studies on stored maize grain have also shown optimum AFB1 production at between 5–10 days at 0.98 aw and 10 days at 0.95 aw [42]. We have compromised and used 10 days in this study so that we can obtain molecular information and relevant toxin data. After 10 days of kernel infection (see above), 1 g of frozen milled maize was ground to powder using a mortar and pestle in the presence of liquid nitrogen, and then placed into a 2 mL extraction tube for isolation of total RNA. Total RNA was extracted using the RNAeasy Plant Mini kit (Qiagen, Hilden, Germany). One hundred milligram of the resulting powder was used for isolation of total RNA. The powder was resuspended in 1 mL lysis buffer supplemented with 10 μL β-mercaptoethanol in a 2 mL RNase free micro reaction tube. After vortexing, the tube was quickly frozen in liquid nitrogen. The sample was then thawed on ice. All further procedures were essentially the same as recommended by the manufacturer’s protocol.

4.5. RNA Sequencing

Library preparation for RNA sequencing was conducted using the NEB Ultra Directional RNA Library Prep Kit. Sequencing was performed on an Illumina HiSeq 2000 instrument. All samples had three biological replicates except the following samples: 30 °C/0.99 aw/350 ppm, 30 °C /0.99 aw/650 ppm and 37 °C/0.99 aw/650 ppm, which had two replicates, and 30 °C/0.99 aw/1000 ppm, 37 °C/0.99 aw/1000 ppm and 37 °C/0.99 aw/350 ppm which consisted of one replicate. Untrimmed sequencing reads were mapped to the A. flavus NRRL3357 (assembly JCVI-afl1-v2.0, http://www.ncbi.nlm.nih.gov/genome/360?genome_assembly_id=28730) reference sequence using GSNAP [43,44]. Reads (www.ncbi.nih.gov, accession ID: PRJNA380582) mapping to exons were counted using featureCounts [45] followed by differential expression testing with DESeq2 [46]. Genes were considered differentially expressed if they had an adjusted p-value < 0.05. Gene ontology and KEGG term enrichment was done using the GOSeq R Bioconductor package. The gene co-expression network was made using WGCNA (Weighted Gene Network Co-expression Analysis) with a signed network, the bi-weight mid-correlation method, and a soft-thresholding power of 10. Variance stabilized counts from DESeq2 were used as input to WGNCA. Genes that had an average read count of less than 10 across all samples were removed prior to WGNCA analysis. Principle Component Analysis plots and heat maps were generated using R. Venn Diagrams were created by BioVenn [47]. Gene ontology term enrichment analysis for WGCNA clusters was conducted using Fungifun2 [48]. Rendering of Gene Networks was performed using Cytoscape v3.4.0 [49]. KEGG annotation was produced through the BlastKOALA [50] service.

4.6. RNA Isolation, cDNA Synthesis and Quantitative PCR

RNA for qPCR was isolated using the RNeasy Plant Mini Kit (Qiagen) from infected maize kernels (see above). cDNA was synthesized according to the Omniscript RT kit protocol (Qiagen) using 500 ng of total RNA. The Bio-Rad CFX96 Real Time PCR Detection System (Bio-Rad, Watford, UK) was used for qPCR with TaqMan probes targeting two AF cluster genes, aflD and aflR, and the housekeeping gene β-tubulin gene [20,38]. Primer and probe sequences were obtained from Medina et al. [17]. Reactions from three biological replicates were prepared in 12.5 μL reaction mixtures in MicroAmp optical 96-well reaction plates and sealed with optical adhesive covers (Bio-Rad). The optimal thermal cycling conditions included an initial step of 10 min at 95 °C and all 45 cycles at 95 °C for 15 s, 55 °C for 20 s and 72 °C for 30 s. Quantification cycle (Cq) determinations were automatically performed by the instrument using default parameters. The expression ratio was calculated using the 2−ΔΔ method [51].