RNA-Seq Revealed Differences in Transcriptomes between 3ADON and 15ADON Populations of Fusarium graminearum In Vitro and In Planta

Fusarium graminearum is the major causal agent of Fusarium head blight (FHB) in barley and wheat in North America. The fungus not only causes yield loss of the crops but also produces harmful trichothecene mycotoxins [Deoxynivalenol (DON) and its derivatives-3-acetyldeoxynivalenol (3ADON) and 15-acetyldeoxynivalenol (15ADON), and nivalenol (NIV)] that contaminate grains. Previous studies showed a dramatic increase of 3ADON-producing isolates with higher aggressiveness and DON production than the 15ADON-producing isolates in North America. However, the genetic and molecular basis of differences between the two types of isolates is unclear. In this study, we compared transcriptomes of the 3ADON and 15ADON isolates in vitro (in culture media) and in planta (during infection on the susceptible wheat cultivar ‘Briggs’) using RNA-sequencing. The in vitro gene expression comparison identified 479 up-regulated and 801 down-regulated genes in the 3ADON isolates; the up-regulated genes were mainly involved in C-compound and carbohydrate metabolism (18.6%), polysaccharide metabolism (7.7%) or were of unknown functions (57.6%). The in planta gene expression analysis revealed that 185, 89, and 62 genes were up-regulated in the 3ADON population at 48, 96, and 144 hours after inoculation (HAI), respectively. The up-regulated genes were significantly enriched in functions for cellular import, C-compound and carbohydrate metabolism, allantoin and allantoate transport at 48 HAI, for detoxification and virulence at 96 HAI, and for metabolism of acetic acid derivatives, detoxification, and cellular import at 144 HAI. Comparative analyses of in planta versus in vitro gene expression further revealed 2,159, 1,981 and 2,095 genes up-regulated in the 3ADON isolates, and 2,415, 2,059 and 1,777 genes up-regulated in the 15ADON isolates at the three time points after inoculation. Collectively, our data provides a foundation for further understanding of molecular mechanisms involved in aggressiveness and DON production of the two chemotype isolates of F. graminearum.


Introduction
Fusarium graminearum is the major causal agent of Fusarium head blight (FHB) in North America and other regions of the world. The pathogen not only causes direct yield losses, but also produces various types of trichothecene mycotoxins [Deoxynivalenol (DON) and its acetylated forms (3-acetyl-4-deoxynivalenol = 3ADON and 15-acetyl-4-deoxvevalenol = 15ADON), nivalenol (NIV) and its acetylated form 4-acetylnevalenol] [1]. DON not only acts as a virulence factor during disease development [2,3] but also poses severe health hazards to both human and animals [4][5][6]. FHB and DON are primarily managed through an integrated approach that combines use of moderately resistant cultivars and a timely fungicide application [7]. However, sources of effective FHB resistance are limited and use of a single source of resistance such as Sumai 3 in commercial cultivars may create selection pressure on the pathogen and lead to an outbreak of more virulent/aggressive pathogen population [8,9].
In recent years, population genetics, global species structure and trichothecene chemotype diversity of the FHB pathogen complex have been extensively studied [10][11][12][13]. These studies indicate that one chemotype of the pathogen is dominant in specific geographic regions although other types may coexist in small fractions. Dominance of NIV-type isolates along with the low frequency of 15ADON-or 3ADON-type isolates is more common in Asian regions [14][15][16]. In North America, dominance of 15ADON-type isolates along with the presence of 3ADON-or NIV-type isolates was observed [17][18][19][20]. However, recent studies indicated that 3ADON-type isolates were significantly increased in Canada [10], North Dakota [18], and China [21] in recent years. The newly emerging 3ADON-type population appears to be more aggressive based on growth rate, disease severity on different cultivars with varied levels of resistance, and DON production in vitro [10,18]. Our data from two years of field experiments using susceptible and moderately resistant wheat cultivars indicated that the 3ADON producers accumulated a higher level of DON on grains irrespective of host resistance after spray inoculation (Ali and Zhong, unpublished data). The recovery of Fusarium isolates from artificially inoculated heads with a mixture of both 3ADON-and 15ADON-type isolates indicated that the recovery frequencies were similar for both types of isolates (Ali and Zhong, unpublished data), suggesting that the 3ADON-type isolates do not have the advantage of outcompeting 15ADON-type isolates during infection. However, it is still not known why the 3ADONtype population produces a higher level of DON than the 15ADON-type population.
Trichothecenes produced by F. graminearum can be broadly categorized into two groups (type A or type B trichothecenes) based on presence or absence of oxygen atoms at carbon atoms 7 (C-7) and 8 (C-8) [1,22]. Type A trichothecenes (T-2 toxin, HT-2 toxin, and 4, 15-diacetoxyscirpenol) lack the hydroxyl group at C-7 but have a hydroxyl group and ester group or no oxygen substitution at carbon atom C-8, whereas Type B trichothecenes (DON, NIV, and their derivatives) have a hydroxyl group at C-7 and a keto (carbonyl) group at C-8 [23][24][25]. DON and NIV have structural differences with the former lacking, and the latter possessing an oxygen atom at carbon C-4 [22,26].
Trichothecene biosynthesis in Fusarium involves a complex pathway consisting of oxygenation, isomerization and esterification steps [22]. The enzymes involved in these biosynthetic steps are encoded by 15 TRI genes located at three loci, including the TRI core cluster with 12 genes (Tri3, Tri4, Tri5, Tri6, Tri7, Tri8, Tri9, Tri10, Tri11, Tri12, Tri13, and Tri14) [27], the two-gene (Tri1-Tri16) locus [28][29][30], and a single gene locus (Tri101) [31]. Whether DON or NIV is produced depends on the function of two genes, Tri13 (cytochrome P450 monooxygenase) and Tri7 (acetyltransferase) [32]. NIV producers contain functional Tri13 and Tri7 while DON producers carry non-functional copies of these two genes due to multiple deletion or insertion events [32,33]. The genetic basis of 3ADON and 15ADON production and their biological significance are still not clear. Recently, Alexander et al. [34] indicated that Tri8 (trichothecene C-3 esterase) regulates the production of 3ADON or 15ADON, and is required to convert the diacetylated 3-and 15-ADON intermediate into 3ADON and 15ADON, respectively. However, comparative studies of transcriptomes of the two trichothecene types (3ADON and 15ADON) during host infection have not been conducted and the molecular mechanisms involved in higher DON production and more aggressive nature of 3ADON isolates are not known.
Gene expression profiles and their relationships with virulence or aggressiveness and DON accumulation during F. graminearum-wheat and -barley interactions have been studied using DNA microarrays [35][36][37][38][39]. With the development of next generation sequencing technologies, new tools such as RNA-sequencing provide more effective approaches to study the gene expression profile changes of organisms under different conditions [40,41]. The RNA-seq method is more sensitive than microarrays especially in detecting those transcripts that are rarely expressed [41]. Wang et al. [41] showed that 8 million reads were sufficient to reach RNA-Seq saturation for most samples with large genome sizes. Bashir et al. [42] have demonstrated that more than 90% of the transcripts in human samples can be adequately covered with just a million sequence reads. The coverage is more than sufficient to reach the saturation needed in RNA-Seq for small genomes such as those of filamentous fungi, while costs are comparable to the DNA microarray approach. Walkowiak et al. [40] used RNA-seq to compare gene expression profiles of one 3ADON strain and one 15ADON strain during interaction in cultures and identified genes expressed differentially between isolates and during their interaction. The overall goal of our study was to understand the molecular mechanisms that make 3ADON-and 15ADON-type populations different during infection on a susceptible cultivar with the following specific objectives: i) to compare the transcriptomes of the 3ADON-and 15ADON-type populations in vitro and in planta using the RNA-seq approach, and ii) to identify the expression differences of candidate genes related to aggressiveness and DON production between the two types of F. graminearum isolates.

Materials and Methods
Fungal isolates and wheat materials F. graminearum isolates were collected during 2008 to 2010 from North Dakota (45°56 0 -49°0 0 0 N, 96°33 0 -104°03 0 W) and genotyped with trichothecene type specific primers [18]. Ten 3ADON-type isolates and ten 15ADON-type isolates (Table 1) were randomly selected and further characterized for disease aggressiveness and DON accumulation on three spring wheat genotypes [Grandin (susceptible), Steele-ND (moderately resistant), and ND 2710 (resistant)] in greenhouse experiments using the point inoculation method [43]. The plant growth conditions, inoculum preparation, inoculation, and disease scoring were the same as previously described [18]. For DON content measurement, grains from inoculated heads of each wheat genotype were harvested at maturity, manually threshed, ground to a fine powder using a coffee grinder, and then sent to the Veterinary Diagnostic Laboratory, North Dakota State University, for mycotoxin analyses.

Sample collection for RNA extraction
To collect the in vitro samples, spores from each of the ten 3ADON-or 15ADON-type isolates ( The in planta samples were collected from the FHB susceptible wheat cultivar 'Briggs' [44] at three different time points after inoculation with the 3ADON-and 15ADON-type populations, respectively. Briefly, seven seeds of Briggs were planted in each of the 15-cm plastic pots filled with Sunshine pot mix (Sun Gro Horticulture Canada Ltd.) and maintained in a greenhouse with 16 h supplemental lights at 23±1°C. Plants were fertilized with slow releasing Osmocote + (15:9:12) (Everris NA, Inc, Marysville OH) and Plantex 20-20-20 (Plant Products Co. Ltd, Ontario, Canada) at a two-week interval. At anthesis, eight to ten heads/pot were tagged and inoculated with a mixture of spores from ten isolates of the same type (Table 1) using the point inoculation method with 1000 spores per spikelet. The spikelets nearby to the inoculated one for each spike were marked with sharpie markers. After incubation in humidity chambers (misting run for 30 s on every 8 m to maintain 100% humidity) for 48 h at 26-27°C and with 18 h lights, the plants were returned to the greenhouse under normal conditions. The inoculated spikelets were collected from 8-10 heads per replicate at 48, 96 and 144 hours after inoculation (HAI), respectively, frozen immediately in liquid nitrogen, and stored at -80°C until use. For each type of isolates, two replicates were used. Thus, a total of 12 in planta samples [two fungal populations (15ADON-type and 3ADON-type), two replicates, and three time points (48,96, and 144 HAI)] were obtained for RNA extraction.

RNA extraction, library preparation, and sequencing
Total RNA was extracted from approximately 30 mg of each in vitro sample (mycelia) or in planta sample (plant and fungal tissues) (ground to a fine powder) using the SV total RNA Mapping sequence reads to the reference genome and identification of differentially expressed genes (DEGs) Mapping of sequence reads to reference genome, and analyses of transcript abundance and differential gene expression were performed as described by Trapnell et al. [45]. All sequence reads were trimmed to remove the low-quality sequences (the first 13 bases). The trimmed reads (37 bases) were then aligned to the F. graminearum reference genome downloaded from the Broad Institute (http://broadinstitute.org/annotation/genome/fusarium_group/ MultiDownloads.html) using Bowtie v0.12.5 (http://bowtie-bio.sourceforge.net/index.shtml) [46] and TopHat v2.0.0 (http://tophat.cbcb.umd.edu/) [45,47] with default settings. Cufflinks v0.9.3 (http://cufflinks.cbcb.umd.edu/; [48]) was used to calculate transcript abundance based on fragments per kilobase of transcript per million fragments mapped (FPKM) using all parameters on default settings. The transcript was considered as expressed when the FPKM value was greater than 0.1 and the lower boundary for FPKM value was greater than zero at 95% confidence interval. Once the transcript abundance was calculated for individual sample files using Cufflinks, the output files were further merged pairwise for each comparison (in vitro comparison between two populations, in planta comparison between two populations and in planta versus in vitro for each population) using Cufflinks utility program-Cuffmerge [45]. The pairwise comparisons of gene expression profiles between the two populations were done using the Cuffdiff program of the Cufflinks version 1.3.0 [48]. The genes were considered significantly differentially expressed if Log2 FPKM (fold change) was !1.0 and false discovery rate (FDR, the adjusted P value) was <0.01. The q-value which is a positive FDR analogue of the p-value was set to <0.01 [49].
In order to visualize the expression data from all samples into two dimensions, principal component analysis (PCA) was performed using JMP Genomics v 6.0 (SAS Institute Inc., Cary, NC) for all genes (except novel transcripts). The expression data were transformed using mean normalization prior to PCA. The expression data of individual conditions were divided by their mean values across all treatment conditions to neutralize influence of hidden factors.

Functional categorization of differentially expressed genes (DEGs)
The differentially expressed genes (DEGs) were functionally categorized online for all pairwise comparisons according to the Munich Information Center for Protein Sequences (MIPS) functional catalogue [50]. The functional categories and subcategories were regarded as enriched in the genome if an enrichment P-and FDR-value was below <0.05. The KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analyses were performed using interface on Blast2GO (Blast2GO v2.6.0, http://www.blast2go.com/b2ghome) for all DEGs to identify gene enrichment on a specific pathway.
Quantitative real time-PCR (RT-PCR) analysis. RT-PCR analysis was done according to the method described by Leng et al. [51]. Briefly, total RNA was extracted from each of the samples as collected for RNA-seq using the SV Total RNA Isolation Kit (Promega, Madison, WI) and purified by treatment with DNase I (NEB, Ipswich, MA) according to the manufacturers' manuals. The reverse transcription reaction was performed on 2μg of total RNA using the SMART MMLV Reverse Transcriptase (Takara, Mountain View, CA). cDNA was diluted 20 times and used as template for quantitative RT-PCR, which was performed with the CFX96 real time PCR system (Bio-Rad, Hercules, California). Primers used for RT-PCR were listed in S5 Table. For each cDNA sample, three replications were performed. Each reaction mixture (20 μl) contained 5μl of cDNA template, 10μl of SYBR 1 Green PCR Master Mix (Applied Biosystems, Foster, CA) and 0.3μl of each primer (10 μM). Relative expression levels of genes were normalized using the beta-tubulin gene as internal control, and were calculated as the fold change by comparison between in planta and in vitro (axenic culture) samples.

Statistical analyses
Significance of differences between 3ADON and 15ADON populations in DON accumulation on each wheat genotype, number of total reads, and total expressed genes were analyzed using proc test in SAS (version 9.4; SAS Institute, Cary, NC). The standard deviation and standard error were calculated using Microsoft Excel 2013. The area-proportional Venn diagrams were created using the online tool (http://bioinforx.com/free/bxarrays/overlap.php) and Microsoft PowerPoint 2013.

Sequence data accessibility
The sequence data have been uploaded to the NCBI Short Read Archive [GEO# GSE83735]

Results
Aggressiveness and DON accumulation of the two fungal populations on spring wheat genotypes No significant difference in aggressiveness was observed between the 3ADON and 15ADON populations on the three wheat genotypes based on disease severity (DS). The average DS on the resistant genotype (ND 2710) was 29.6±7.2% and 27.9±10.1% for the 3ADON and 15ADON populations, respectively. On the susceptible cultivar (Grandin), the average DS was 65.3±11.4% for the 3ADON population and 56.8±16.3% for the 15ADON population. On the moderately resistant cultivar (Steele-ND), the average DS was 52.0±9.5% for the 3ADON population and 54.2±7.1% for the 15ADON population. However, the 3ADON population accumulated a significant higher DON content on Grandin (t = 8.1, p< 0.0001) and ND 2710 (t = 3.4, p = 0.0034) than the 15ADON population, although the difference was not significant on Steele-ND (t = 0.8, p = 0.44) (Fig 1).

General features of in vitro and in planta transcriptomes of 3ADON and 15ADON populations
A total of 559,577,636 sequence reads totaling 27.97 gigabase-pair (Gbp) were generated from the 4 in vitro and 12 in planta RNA samples ( Table 2). The numbers of sequence reads from the in-vitro RNA samples ranged from 34.8 to 40.9 million for the 3ADON population, and 32.2 to 34.5 million for the 15ADON population. The sequence reads from the in planta samples ranged from 26.4 to 49.5 million for the 3ADON population, and 27.8 to 39.1 million for the 15ADON population. No significant differences in a total number of sequence reads were observed between the two populations (t = 0.99, p = 0.34).
Of the sequence reads from in vitro samples, 80.1-82.4% were aligned to the reference genome of F. graminearum (PH-1), and 79.1-81.7% of the reads had a single match to the reference genome while 1.07% of the reads aligned to multiple genomic locations within the reference genome. Of the sequence reads from in planta samples, only 5.3-13.3% mapped to the reference genome, with unique matches accounting for 5.1-12.5% of total reads ( Table 2). The number of expressed transcripts ranged from 14,242 to 14,564 in vitro, and 12,163 to 13,586 in planta (Fig 2).
The principal component analysis (PCA) of transcript abundance, measured as the FPKM (fragments per kilobase of transcript per million fragments mapped) values, indicated that the RNA-seq data were very similar between the two biological replicates of each treatment (Fig 3). The principal component 1 (PC1) accounted for 21.6% variation, and clearly distinguished the

Comparison of transcriptomes between 3ADON and 15ADON populations under in vitro conditions
Under in vitro growth conditions, 479 and 801 genes were up-and down-regulated, respectively, in the 3ADON population as compared to the 15ADON population. Among them, 454 up-regulated and 723 down-regulated genes were identified in Munich Information Center for Protein Sequences (MIPS) database using functional catalogue (FunCat) analysis (Table 3, S1  Table). Among the 454 up-regulated genes found in FunCat, 276 were under 'unclassified proteins' and were non-significant (p = 0.977, S1 Table). The significantly enriched categories were those for metabolism (p = 1.8×10 −09 ), cellular transport, transport facilities and transport  3ADON and 15ADON populations in vitro and in planta. Transcript-fragment reads from each sample were mapped to the reference genome (F. graminearum PH-1) using Bowtie version 0.12.5 [46] and TopHat version 2.0.0 [47]. Fragments per kilobase of transcript per million fragments mapped (FPKM) were calculated by Cufflinks version 0.9.3 [48]. HAI: hours after inoculation, 3ADON: the 3ADON population producing 3-acetyl-deoxynivalenol and DON, 15ADON: the 15ADON population producing 15-acetyl-deoxynivalenol and DON, vertical bar represents standard error of means between replications, MBA: mung bean agar (in vitro).
KEGG pathway analysis of differentially expressed in vitro genes indicated that majority of up-regulated genes in the 3ADON population were involved in metabolism of starch, sucrose, methane, drugs, and in inter-conversions of pentose and glucuronate, while the down-regulated genes were those for purine, and thiamine metabolism (Fig 4, S1 Table).  Table 3). More unique genes were expressed during early infection establishment, but the number of common genes was higher during disease spread even for the 3ADON population ( Fig 5D, S3 Table).
In order to identify uniquely up-regulated genes and their functional categorization in each population, the genes specifically up-regulated in the 3ADON and 15ADON populations at each time point (633+526+688 in 3ADON; 889+604+350 in 15ADON) were summed up and duplicate genes among three time points were discarded. This gave a total of 1396 unique genes in the 3ADON population, and 1398 unique genes in the 15ADON population, respectively. Of these in planta up-regulated unique genes, 1257 from 3ADON population and 1278 from 15ADON population were identified in Munich Information Centre for Protein Sequences (MIPS) functional catalogue database [50]. The highest number of genes from both 3ADON (552/1257, 43.9%) and 15ADON (516/1278, 40.3%) populations belonged to 'unclassified proteins' category ( Fig 6, S4 Table). In the 3ADON population, the up-regulated genes with functional annotations were significantly enriched in the following categories: 'metabolism' (400 genes, p = 0), 'protein synthesis' (134 genes, p = 1.2×10 −38 ), and 'protein with binding function or cofactor requirements' (structural or catalytic) (339 genes, p = 1.7×10 −5 ). In the 15ADON population, the major enriched categories were 'metabolism' (536 genes, p = 0), 'energy' (77, p = 0.0017), and 'cell rescue, defense and virulence' (162 genes, p = 0.0014) (Fig 6, S4 Table). Also significantly enriched were the major functional sub-categories involved in nitrogen, sulfur and selenium metabolism (p = 0.0199), pentose-phosphate pathway  Table).
Further analyses of the same gene sets using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway identified 205 genes involved in 65 metabolic pathways for the 3ADON population, and 178 genes involved in 72 metabolic pathways for the 15ADON population (Table 4, S4 Table). In the 3ADON population, more genes involved in pathways for Values with log2 fold change >1 and false discovery rate (<0.01) were considered as differentially expressed. Only pathways having at least two genes up-regulated on either of the population are shown. Detailed information on KEGG pathway analysis is given in S1 Table. doi:10.1371/journal.pone.0163803.g004 metabolism of purine, arginine, proline and pyrimidine, citrate cycle (TCA cycle), valine, leucine and isoleucine biosynthesis, pentose phosphate pathway, carbon fixation pathways in prokaryotes or in photosynthetic organisms were up-regulated (Table 4). In contrast, in the 15ADON population, more genes involved in nitrogen metabolism, pentose and glucuronate inter-conversions, starch and sucrose metabolism, drug metabolism-cytochrome P450, N-Glycan biosynthesis, various types of N-glycan biosynthesis, benzoate degradation, chloroalkane and chloroalkene degradation, metabolism of xenobiotics by cytochrome P450 were up-regulated (Table 4, S4 Table).
All TRI genes (Tri1, Tri3, Tri4, Tri5, Tri6, Tri8, Tri9, Tri10, Tri11, Tri12, Tri14, and Tri15) involved in the biosynthesis or regulation of trichothecene production were highly up-regulated in planta in both populations as compared to the in vitro conditions (Fig 7, the values are shown in Log2 scale). The highest fold increase was observed for the Tri3 gene in both populations followed by Tri4, Tri5, Tri1 and others (Fig 7A and 7B). Five Tri genes (Tri1, Tri6, Tri10, Tri11, and Tri12) showed higher expression in the 3ADON population as compared to the 15ADON population. The result is consistent with quantitative RT-PCR analysis of mRNA expression for two of the Tri genes (Tri6 and Tri10) (Fig 7C and 7D).
In general, the expression levels of the genes involved in biosynthesis of trichothecenes and secondary metabolites (polyketides and non-ribosomal peptides) varied in the in vitro and in planta samples collected at the three time points. While, as stated above, most of the trichothecene biosynthesis genes were up-regulated in planta compared to in vitro, some other secondary metabolite genes (NPS4, NPS7, NPS11, PKS6, PKS11, and PKS12) had a higher expression in vitro than in planta in both populations (Fig 8). However, there was an apparent difference in Tri gene expression pattern between the two chemotypes. Transcripts of some Tri gene were high at 48 HAI, but declined at 96 HAI and recovered at 144 HAI in 3ADON population while transcripts of most Tri genes in 15ADON population gradually diminished during the plant  Table. doi:10.1371/journal.pone.0163803.g006 Table 4

Comparison of transcriptomes between 3ADON and 15ADON populations in planta
A total of 185, 89, and 62 genes were up-regulated in the 3ADON population while 292, 361, and 241 genes were down-regulated at 48, 96 and 144 HAI, respectively, as compared to the 15ADON population (Table 3). Among them, 167, 63 and 44 genes were uniquely up-regulated at 48, 96 and 144 HAI, respectively. Only four genes (FGSG_04621, FGSG_04694, FGSG_06540, and FGSG_10632) were up-regulated at all three infection time points (Fig 9, S2  Table).   Functional catalogue (FunCat) analysis [50] indicated that more than 50% of differentially expressed genes identified in planta were for un-classified proteins, and 20.9% (37/177) and 31.2% (89/285) of the up-and down-regulated genes at 48 HAI were in the functional category 'metabolism' (Table 5, S2 Table). The functional category 'cell rescue, defense and virulence' was only significantly enriched among up-regulated genes at 96 HAI and among down-regulated genes at 48 HAI, respectively (S2 Table). Among the up-regulated genes, the functional sub-categories in glutamine degradation (p = 0.0055), C-compound and carbohydrate metabolism (p = 0.011), glutamine metabolism (p = 0.0153), arginine biosynthesis (0.0241), metabolism of urea cycle, creatine and polyamines metabolism (0.036), lipid, fatty acid and isoprenoid metabolism (0.046), assimilation of ammonia, and metabolism of the glutamate group (0.0499) were significantly enriched at 48 HAI. In addition, the genes for non-vesicular and cellular import, transport of compounds such as allantoin and allantoate, cofactor and sugar, and for post-transcriptional control were also significantly enriched (p<0.05) during the early infection stage (48 HAI). At 96 HAI, the genes for detoxification either by degradation (p = 0.0028) or by modification (p = 0.0062), and for defense related proteins (p = 0.0098) were highly enriched. At 144 HAI, the genes for secondary and polyketides metabolism, acetic acid derivatives metabolism, detoxification by modification or by degradation, and for degradation of ester compounds were abundant (P<0.05) (S2 Table).  Table. doi:10.1371/journal.pone.0163803.g009 A total of 50 genes in the 3ADON population showed at least five-fold expression differences at one or more time points in planta as compared to the 15ADON population, with 36 of them exhibiting higher expression and 24 having lower expression (Table 6). Among these 50 differentially expressed genes, 19 and 12 genes were up-regulated in the 3ADON and 15ADON populations, respectively, at least at one time point in planta in comparison with in vitro expression (Table 6).

Discussion
In this study, we compared relative aggressiveness and amount of DON accumulation in wheat grains between the 3ADON and 15ADON populations of F. graminearum, and also examined differences in transcriptomes between the two populations under in vitro and in planta conditions. The 3ADON population caused a higher level of disease than the 15ADON population although the difference was not significant. This result is in agreement with that of Gilbert et al [54], who compared aggressiveness of 3ADON and 15ADON isolates on the susceptible cultivar 'Robin' and the moderately resistant line '5602 HR' , and failed to find significant differences in aggressiveness. However, we found that the 3ADON population produced a significantly higher level of DON on ND 2710 and Grandin, but not on Steele ND (Fig 1). This result is consistent with previous studies [18,55], which showed that the 3ADON population accumulated a significantly higher level of DON than the 15ADON population on both susceptible and resistant cultivars in greenhouse and field experiments.  Using microarrays, Lysøe et al. [35] studied the global gene expression pattern of F. graminearum (a 15ADON producer) during infection of the susceptible wheat cultivar 'Bobwhite' , and showed the number of expressed genes increased from 48 hai (>4000 genes) to 96 hai (>8000), but declined after 144 hai. However, our data showed the number of expressed genes was relatively consistent for the 3ADON population at the three time points although the number of expressed genes decreased as the infection progressed for the 15ADON population ( Fig  2). The discrepancy between our study and Lysøe et al. [35] could be due to different methods used for RNA sampling and expression measurements. We collected only inoculated spikelets at each infection time point for RNA extraction, while Lysøe et al. [35] harvested whole spikes for RNA isolation. The higher number of expressed genes found at a later stage for the 3ADON population may reflect its more vigorous aggressiveness and higher DON production than the 15ADON population during infection. In contrast to microarrays, the RNA-seq method potentially quantifies all of the expressed RNA from both the pathogen and the plant. We detected more than 12,100 gene transcripts at all three infection stages studied while Lysøe et al. [35] detected less. The significantly higher number of transcripts detected in our study than in Lysøe et al. [35] may be due to the more sensitive nature of RNA-seq than the microarray technique for detecting rarely expressed transcripts [41,56].
The comparative analyses of in vitro transcriptomes between the 3ADON and 15ADON populations identified differentially expressed genes associated with several physiological and cellular metabolic processes. Remarkably, the genes annotated for cellular transport and transport facilities and those involved in active metabolism of internal metabolites and on nutrients uptake were significantly enriched in the 3ADON population. Seong et al. [57] identified a considerably higher number of genes (216) annotated for permeases or transporters during conidial germination (0-24h), and in fresh spores and hyphae under nutrient limiting conditions. Hallen et al. [58] analyzed changes in gene expression during perithecial development. Among 162 predicted ions transporter genes, 44 were upregulated at least two fold during sexual development. In many fungi, nutrient deficiency is prominent phenomenon during sporulation and early infection process that requires transport of various nutrients, nitrogen and carbon sources [59]. Thus, expression and enrichment of transporter genes in 3ADON population might be essential for metabolism and uptake of various carbon and nitrogen compounds during the nutrient starvation period to facilitate higher spore yield. FGSG_08403 is one of 12 highly upregulated TFs in 3ADON population. This gene is required for perithecial development and ascospore production in F. graminearum [60]. It is not known if the higher expression of this gene in 3ADON population would be an advantage over the 15ADON population. Further investigation is needed to answer this question.
We found none of the Tri genes except Tri15 and Tri8 showed differential expression between the two populations under in vitro growth conditions. Tri15 (FGSG_11025) showed a 2.8 fold increased expression while Tri8 (FGSG_03532) had a 7.7 fold decreased expression in 3ADON population as compared to the 15ADON population (S1 Table). This result is in contrast with Walkowiak et al. [40] who found all of Tri genes except Tri10 were upregulated in the 15ADON isolate (FG2) compared to the 3ADON isolate (FG1) in vitro. The discrepancy might be due to use of different growth conditions and F. graminearum isolates in the two studies. We used five days old mycelia grown on MBA plate while Walkowiak et al. [40] harvested mycelia from liquid culture after six hours of incubation for RNA extraction. Additionally, we used two groups of isolates for gene expression comparisons while only one pair of isolates was compared in the study of Walkowiak et al. [40].
We identified a large number of genes up-regulated in planta compared to in vitro in 3ADON and 15ADON populations (Fig 5A-5C, S2 Table). In general, all known Tri genes required for trichothecene biosynthesis and genes involved in secondary metabolite production were differentially expressed in vitro versus in planta in both populations (Figs 7 and 8). However, the dynamic expression changes of Tri genes during the three infection time points were different for the two types of isolates. Most of Tri genes in the 3ADON isolates showed an up-, down-, and up-regulated expression pattern at 48, 96, and 144 HAI compared to the in vitro conditions, whereas expression of most Tri genes in the 15ADON isolates was gradually decreased as the infection progressed (Figs 7 and 8). Among these genes, Tri5 encodes trichodiene synthase involved in the first step in the trichothecene biosynthetic pathway and is required for DON synthesis [61], while Tri8 encodes C-15 esterase or C-3 esterase, which is required for conversion of intermediate precursor into 3ADON or 15ADON product, respectively [34]. Also, four Tri genes, namely Tri1, Tri6, Tri10 and Tri11, had at least 1.4 times greater expression in 3ADON population compared to 15ADON population across all time points (Fig 7). Tri1, which is located outside the Tri-cluster, encodes P450 oxygenase and negatively regulates the production of calonectrin [62], the intermediate precursors for 3ADON biosynthesis in F. culmorum [63]. Tri6 is a pathway-specific transcriptional regulator in trichothecene biosynthesis pathway [30,64], and further regulates expression of six other Tri genes (Tri1, Tri3, Tri6, Tri7, Tri12 and Tri14) within the Tri-cluster and additional 192 potential genes in F. graminearum [65]. The Tri10 gene, another regulatory factor, is required for trichothecene biosynthesis and regulation of expression of six Tri genes (Tri3, Tri7, Tri8, Tri9, Tri11 and Tri12) in F. sporotrichioides [30]. Tri11 encodes cytochrome P-450 monooxygenase required for hydroxylation in trichothecene biosynthesis [62]. The higher expression of these global regulating genes in 3ADON population might explain the difference in DON production between 3ADON and 15ADON populations.
Functional analysis of in planta up-regulated genes compared to in vitro between the two populations further revealed a chemotype specific gene expression pattern, and identified genes involved in distinct metabolic and molecular mechanisms. Major enriched functional categories identified were consistent with previous study [35], although some of the categories were specific to either 3ADON or 15ADON population. Genes under functional categoriesnitrogen, sulfur and selenium metabolism; protein synthesis; protein with binding function or cofactor requirements and sub-categories under metabolism were highly specific and enriched in 3ADON population. Within the metabolism sub-category, the genes involved in the assimilation of ammonia, metabolism of the glutamate group, degradation, and biosynthesis of amino acids were highly enriched. It is a well-established concept that nutrients availability and their acquisition by pathogens are pre-requisites for successful colonization and fungal establishment [59]. Various nitrogen and carbon sources are required for trichothecene biosynthesis, secondary metabolite production, and virulence in many fungi including F. graminearum [66][67][68]. The nitrogenous compounds such as ammonia, glutamine, glutamate, and asparagine are the primary nutrient sources for many fungi including S. cerevisiae, A. nidulans, N. crassa and others. However, in case of lack or very low concentration of primary sources, fungi utilize many alternative nitrogen sources such as nitrate, nitrite, purines, amides, most amino acids and proteins after de novo secretion of pathway-specific catabolic enzymes and permeases [69,70]. Thus, the expression of genes to utilize these nutrient sources during infection might have fitness advantage to the 3ADON population over the 15ADON population in survival, sporulation, aggressiveness, and higher DON accumulation.
Members of the C2H2 (Cys-Cys-His-His) zinc finger transcription factor (TF) family were expressed more abundantly in the 3ADON population than in the 15ADON population (S4 Table). Transcription factors have a diverse role in signal transduction, respiration, nitrogen utilization, peroxisome proliferation, stress tolerance, drug resistance, gluconeogenesis, sugar and amino acid metabolism and so on [71]. Among the 76 TFs with C2H2 zinc finger domains identified in F. graminerum genome [52], eight (FGSG_00764, FGSG_01298, FGSG_01350, FGSG_04288, FGSG_06701, FGSG_10350, FGSG_10470, FGSG_13964) were exclusively upregulated in the 3ADON population (S4 Table). The function of the C2H2 zinc finger proteins in F. graminearum is unknown, but their roles have been studied in other fungi, such as calcium signaling in Aspergillus nidulans (CrzA) [72], regulation of biological processes (sexual development) (SteA) in Aspergillus [73], ustilagic acid biosynthesis (Rua1) in Ustilago maydis [74]. Of the TF encoding genes characterized by Son et al. [60] using deletion mutation, five genes (FGSG_00764 and FGSG_01298 with C2H2 zinc finger domain, FGSG_09286 and FGSG_10142 with bZIP domain, and FGSG_09871 with bromo domain) up-regulated in the 3ADON population were involved in either virulence or DON biosynthesis or both. Thus, the 3ADON population may have unique regulation of genes making it different from the 15ADON in aggressiveness and DON production as indicated by Puri and Zhong [18] and Ward et al. [10].
Our study further revealed a set of transporter-encoding genes uniquely up-regulated and enriched in the 3ADON population (S2 Table, Fig 6, Table 5). The genes involved in the transport of carbohydrate, sugars, allantoin and allantoate, vitamine/cofactor etc. and those required for host invasion and utilization of nutrient sources such as carbohydrates, proteins, lipids, and vitamins were enriched during early infection (48HAI). The genes for carbohydrate transport and cellular import were also enriched at 144 HAI. The allantoin and allantoate transport category, which is required to utilize uric acid, a host induced catabolic compound in response to pathogen infection, was specific to wheat infection under nutrient limiting condition [35]. The enrichment of these transporter genes at 48 HAI suggests their importance during early infection. The genes for cell rescue, defense and virulence were up-regulated in 3ADON population at 96 HAI. The genes for biosynthesis of secondary metabolites and detoxification of antimicrobial plant metabolites were also significantly up-regulated in 3ADON population at 144 HAI. One of them is the cytochrome P450s, which is required for oxygenation during secondary metabolite production and contributes to fungal virulence via detoxification of antimicrobial plant metabolites [75].
Taken together, our results revealed a set of genes that were differentially expressed between 3ADON and 15ADON populations during in vitro, in planta or both conditions. Future functional analysis of these genes will provide insights into the mechanisms involved in higher DON production and aggressive behavior of the newly emerging 3ADON population.
Supporting Information S1