Integrated Metabolo-Transcriptomics Reveals Fusarium Head Blight Candidate Resistance Genes in Wheat QTL-Fhb2

Background Fusarium head blight (FHB) caused by Fusarium graminearum not only causes severe losses in yield, but also reduces quality of wheat grain by accumulating mycotoxins. Breeding for host plant resistance is considered as the best strategy to manage FHB. Resistance in wheat to FHB is quantitative in nature, involving cumulative effects of many genes governing resistance. The poor understanding of genetics and lack of precise phenotyping has hindered the development of FHB resistant cultivars. Though more than 100 QTLs imparting FHB resistance have been reported, none discovered the specific genes localized within the QTL region, nor the underlying mechanisms of resistance. Findings In our study recombinant inbred lines (RILs) carrying resistant (R-RIL) and susceptible (S-RIL) alleles of QTL-Fhb2 were subjected to metabolome and transcriptome profiling to discover the candidate genes. Metabolome profiling detected a higher abundance of metabolites belonging to phenylpropanoid, lignin, glycerophospholipid, flavonoid, fatty acid, and terpenoid biosynthetic pathways in R-RIL than in S-RIL. Transcriptome analysis revealed up-regulation of several receptor kinases, transcription factors, signaling, mycotoxin detoxification and resistance related genes. The dissection of QTL-Fhb2 using flanking marker sequences, integrating metabolomic and transcriptomic datasets, identified 4-Coumarate: CoA ligase (4CL), callose synthase (CS), basic Helix Loop Helix (bHLH041) transcription factor, glutathione S-transferase (GST), ABC transporter-4 (ABC4) and cinnamyl alcohol dehydrogenase (CAD) as putative resistance genes localized within the QTL-Fhb2 region. Conclusion Some of the identified genes within the QTL region are associated with structural resistance through cell wall reinforcement, reducing the spread of pathogen through rachis within a spike and few other genes that detoxify DON, the virulence factor, thus eventually reducing disease severity. In conclusion, we report that the wheat resistance QTL-Fhb2 is associated with high rachis resistance through additive resistance effects of genes, based on cell wall enforcement and detoxification of DON. Following further functional characterization and validation, these resistance genes can be used to replace the genes in susceptible commercial cultivars, if nonfunctional, based on genome editing to improve FHB resistance.

caused by Fusarium graminearum Schwabe [telomorph: Gibberella zeae Schw. (Petch)] is one of the most devastating and alarming diseases of wheat ruining harvests, in many wheat producing regions of the world, including Canada [15,16]. The accumulation of mycotoxins such as deoxynivalenol (DON) and nivalenol (NIV) is of major concern due to their detrimental effects on humans and animals [17]. The development of FHB resistant cultivars is considered to be the best way to manage this disease and the accumulation of mycotoxins in grains, as it is the most efficient, economic, and ecofriendly approach to manage FHB [15].
FHB resistance is quantitative in nature, involving several genes, each with small or large effects, and the phenotype is the result of their additive effects. Three different types of FHB resistance in wheat have been recognized and used in breeding: (i) resistance to initial infection or spikelet resistance (type-I), (ii) resistance to spread within the spike or rachis resistance (type-II), and (iii) resistance to mycotoxin accumulation in grains (type-III) [18]. The development of resistant cultivars is very challenging because of limited understanding of genetics of resistance and lack of cost-effective means to phenotype [19]. The screening is generally done based on spray inoculation which leads to high experimental errors, leading to inconsistent ranking of genotypes over years. Inoculation under controlled conditions can significantly reduce the experimental error and can enable quantification of both spikelet and rachis resistance [20].
The QTL-Fhb2 localized on the short arm of chromosome 6B is the second major QTL conferring rachis resistance [23]. The QTL-Fhb2 has been mapped as a Mendelian factor, spanning a region of 4.2 cM flanked by two simple sequence repeat (SSR) markers, GWM-133 and GWM-644, using recombinant inbred (RIL) population derived from a cross between BW278 (resistant parent) and AC foremost (susceptible parent), and the QTL explained 24.1% of resistance to FHB [23]. Resistant parent BW278 is a descendant of Sumai-3, a Chinese bread wheat cultivar which possesses high level of rachis resistance [22]. Similarly, the QTL-Fhb2 was mapped on 6BS using double haploid population, which conferred 21% resistance to FHB [24]. Both the studies reported high levels of rachis resistance, but no study has reported the specific genes localized in the QTL-Fhb2 region conferring resistance, nor the mechanisms of resistance. As the QTL regions contain several genes, the dissection and functional characterization of each gene localized in the QTL region on genomic scale is a very challenging task, especially in wheat which possesses a highly complex genome and lacks a complete genome sequence [25]. Therefore, genes in the QTL regions and the resistance mechanisms governed by them are not studied in detail. Even though significant attempts have been made to identify candidate genes localized within the QTL-Fhb1 [6,26], QTL-Fhb5 [27], so far no FHB resistance gene has been identified and validated.
Technological advancements in genome sequencing and integration of omics platforms have offered novel insights to explore the regulation of metabolic pathways and their biosynthetic genes underlying disease resistance mechanisms [6,[28][29][30]. Metabolomics is a potential post-genomics tool to elucidate the host biochemical responses under biotic stress, to identify the candidate R genes, and to validate gene functions [6,[31][32][33]. Non-targeted metabolomics has been applied to reveal the host biochemical mechanisms of quantitative resistance in crop plants such as wheat [6,34], and barley against F. graminearum [35][36][37][38], and potato against Phytophthora infestans [28][29][30]33]. Non-targeted metabolic profiling of wheat near isogenic lines (NILs) with FHB resistant QTL-Fhb1 revealed deposition of HCCAs in the cell wall that reduced further spread of pathogen within rachis, thus imparting resistance [6]. Higher abundance of several resistance metabolites belonging to phenylpropanoids, flavonoids, fatty acids, and terpenoids that inhibited the growth of F. graminearum and trichothecene biosynthesis was identified in barley against F. graminearum [35,38]. Metabolic profiling of resistant and susceptible potato cultivars against late blight identified phenylpropanoids and their biosynthetic genes regulated by StWRKY1 [28].
The transcriptome is highly active and instantly changes with response to cellular perturbations. The study of wheat transcriptome under Fusarium stress revealed the expressed genes following FHB infection [39,40]. QTL-Fhb1 specific RNA-seq of Wangshuibai and its mutant NAUH117 (lacking a chromosome region including QTL-Fhb1 segment), revealed association of PR5, PR14, ABC transporter and jasmonic acid pathway genes in FHB resistance in wheat [39]. Transcriptome analysis using RNA-seq in maize upon Fusarium verticillioides inoculation revealed the involvement of shikimate, flavonoid, lignin and terpenoid biosynthetic pathways in imparting FHB resistance [40]. A lipid transfer protein (LTP) was found to be constitutively more abundant in NIL carrying QTL-Fhb5, based on microarray [27]. QTL-specific microarray analysis of Sumai-3 and two susceptible NILs showed up-regulation of 25 genes and the genes encoding PR proteins, like β-1-3 glucanase (PR-2), thaumatin like proteins (PR-5) and wheatwins (PR-4) were significantly over-expressed in genotypes containing Sumai-3 3BS region [41]. Microarray analysis of near isogenic lines carrying QTL-3BS, showed the up-regulation of genes involved in cell wall biogenesis upon fusarium infection [42]. A gene UDP-glycosyltransferase was reported to be highly over-expressed in NILs harboring two QTL-Fhb1 and QTL-Fhb5, based on microarray analysis, which has a major role in the detoxification of deoxynivalenol [42].
Considering this background, RILs carrying resistant and susceptible alleles of genes in QTL-Fhb2 were subjected to metabolome and transcriptome profiling upon F. graminearum inoculation. Our study is the first report that revealed six putative resistance genes localized within the QTL-Fhb2 region and also the plausible association of cell wall enforcing metabolites, explaining the underlying mechanisms of resistance, based on the integration of metabolomics and transcriptomics. The application of these genes, following validation, is discussed.

Disease severity of RILs
Spikelets showing necrotic lesions and bleaching symptoms were considered as diseased. The diseased symptoms were visible at 3 dpi as small, tiny necrotic spots on the inoculated pair of spikelets. The numbers of spikelets diseased in R-RIL were very few; the fungus was able to colonize only spikelets adjacent to the inoculated pair of spikelet, not further, even at 21 dpi. On other hand, in S-RIL the fungus was able to spread through rachis making the whole spike diseased at 21 dpi (Fig 1A). The rachis was lush green in R-RIL without any disease symptoms in rest of the spikelets, unlike in S-RIL where both rachis and spikelets were diseased showing blackish brown discoloration with necrotic lesions. This clearly demonstrated that the QTL imparts high rachis resistance as the fungus was unable to spread through rachis in R-RIL. The PSD in S-RIL was highly significant than in R-RIL ( Fig 1B). The AUDPC calculated from PSD was significantly higher (P < 0.001) in S-RIL (AUDPC = 5.95) than in R-RIL (AUDPC = 1.34) (Fig 1C).

Quantity of fungal biomass in RILs
The quantity of fungal biomass, quantified as relative copy number of Tri6 gene, was highly significant in rachis of S-RIL as compared to R-RIL ( Fig 1D). The reduced amount of fungal biomass in rachis of R-RIL, due to reduced ability of the pathogen to colonize the rachis, clearly demonstrates the role of QTL-Fhb2 in the R-RIL. Whereas, in the S-RIL the Tri6 gene copy number was high, as the fungus was able to move through the rachis to the adjacent spikelets, meaning absence of rachis resistance. The relative gene copy number of Tri6 was also calculated for the spikelets collected from both R-RIL and S-RIL, and it was found that the Tri6 gene copy number was also significantly higher in S-RIL than in R-RIL, but the fold change was . Spike and rachis of R-RIL and S-RIL, 21 dpi with F. graminearum spore suspension. A single alternate pair of spikelets in a spike was inoculated; black arrows indicate the site of inoculation. The spike and rachis in R-RIL shows only necrotic spots or diseased symptoms limited to the inoculated spikelet, while in S-RIL both spikelet and rachis are entirely diseased. (B). Proportion of Spikelets Diseased (PSD). A single pair of spikelets of a spike was inoculated in both the RILs and the proportion of spikelets diseased was recorded at 3 days intervals until 21 days, from which the PSD was calculated. The bar graph shows high PSD in S-RIL as compared to R-RIL. (C). The bar graph shows area under disease progress curve (AUDPC) calculated from PSD, significantly higher in S-RIL. (D) and (E). Fungal Biomass quantification in RILs. Three alternate pair of spikelets were inoculated with F. graminearum spore suspension and samples were collected at 7 dpi. The total genomic DNA was extracted and the relative gene copy number of Tri6 was estimated using 2 −ΔΔC T method. (D). Shows the relative gene copy number of Tri6 in rachis tissues; (E). Shows the relative gene copy number of Tri6 in spikelet tissues. In both graphs, the gene copy number of Tri6 is significantly higher in S-RIL as compared to R-RIL.
comparatively lower in spikelets than in rachis ( Fig 1E). These results clearly demonstrate that the resistant alleles of genes in QTL-Fhb2 impart high rachis resistance.

Resistance related constitutive (RRC) metabolites associated with QTL-Fhb2
A total of 550 RRC metabolites were differentially accumulated, of which 57 were putatively identified (S1 Table). These metabolites belonged to different chemical groups (Table 1;  Comparative transcriptome of RILs based on RNA-seq Whole transcriptome analysis (RNA-seq) was done with two genotypes, viz, R-RIL and S-RIL, across four treatments (RP, RM, SP, SM) with three biological replicates for each, collected at 48 hpi. The assembly of sequences, from the 12 sequenced samples, were compared and annotated. Transcriptome analysis using RNA-seq generated 155012 and 165036 transcripts in R-RIL and S-RIL, respectively. The number of transcripts up-regulated and down-regulated in R-RIL and S-RIL are depicted in Fig 3A. The gene ontology analysis classified the transcripts up-regulated in R-RIL and S-RIL based on their involvement in biological processes, molecular functions and cellular component in which they are localized (S1 Fig). The majority of the transcripts depicted their involvement in biological processes such as translation, transcription, response to oxidative stress, lipid metabolism, photosynthesis, DNA repair, protein folding, carbohydrate metabolic processes, trans-membrane transport, suggesting the involvement of various pathways and/or interactions in disease resistance. Differential gene expression (DGEs) analysis classified transcripts as differentially expressed with Log 2 FC values, or   (Table 2; S2 Table). The differentially expressed transcripts (highly upregulated and highly down-regulated) in R-RIL and S-RIL are shown in the form of heat maps with the respective gene IDs (Fig 3B). Pathway analysis of transcripts showed that majority of transcripts belonged to phenylpropanoid, flavonoid, fatty acid, oxylipin/jasmonic acid, and phospholipid pathways.
Genotype-specific transcriptional changes in response to F. graminearum The genotype-specific transcriptional modulations are clear indications from the gene ontology classification (S1 Fig), showing the differences in the functional categories of the transcripts in each RIL. The R-RIL dataset, in particular, upon pathogen treatment are considered best to identify the transcripts that would possibly be involved in imparting FHB resistance. Nevertheless, in S-RIL pathogen inoculated dataset, the transcripts down-regulated were not overlooked. The transcripts up-regulated in R-RIL were further classified according to their involvement in various biosynthetic pathways or regulators of gene expression (transcription factors, protein kinases, and secondary messengers) or transcripts belonging to RR-proteins involved in DON detoxification ( Table 2). The transcription factor bHLH041 (FPKM = 0.42) was detected only in RP, suggesting its role in FHB defense. Apart from bHLH, transcription regulatory genes like WRKY, R2R3 MYB and MYB-4 were up-regulated in R-RIL ( Table 2). The transcripts belonging to phenylpropanoid pathway genes such as agmatine coumaroyltransferase   , and chitinases were also up-regulated in RP. Transcripts involved in the detoxification were highly up-regulated in RP such as, UDP-glycosyltransferases, multidrug resistance proteins, pleiotropic drug resistance proteins, ABC transporters and glutathione S-transferases (Table 2; S2 Table). All the transcripts up-regulated in resistant genotype were reported to be involved in FHB resistance, thus implicating the involvement of a hierarchy of genes and/or interactions in imparting FHB resistance in wheat.

Genetic controls underlying QTL-Fhb2
The markers flanking the QTL were sequenced and the region (sequence) within the two flanking markers was considered as QTL-Fhb2 region using wheat survey sequence available (http:// wheat-urgi.versailles.inra.fr/). The transcripts aligning to the 6BS reference genome were pulled out separately and furthermore, the transcripts aligning to the sequence within the flanking marker co-ordinates were considered as genes localized within QTL-Fhb2 region. Based on high FC metabolites and transcripts in R-RIL we were able to localize six putative candidate genes within the QTL-Fhb2 region that were associated with biotic stress resistance functions (Fig 4). The putative candidate genes localized within the region were: 4CL, bHLH041 TF, GST, ABC-4, CS, and CAD. The expression values for these genes localized within the QTL-Fhb2 region are presented in Table 2 (transcripts marked with asterisk ( Ã )). The list of all the genes localized within the QTL-Fhb2 region is provided in S3 Table.

Confirmation of gene expression based on qRT-PCR
To validate the RNA-seq data a qRT-PCR analysis was carried out for a few selected genes, such as, 4CL, bHLH041, ABC-4, GST, chitinase1, CHS, PAL and MYB4. The expression values obtained from qRT-PCR analysis for the selected genes were similar to RNA-seq confirming the reproducibility of transcriptome data (Fig 5).

Discussion
Resistance in wheat against FHB is quantitative and more than 100 QTLs with major or minor effects have been identified and mapped on all wheat chromosomes, expect on 7D. However, the genetic controls underlying these resistant QTLs are yet to be revealed. These QTL regions, even when they are fine mapped, contain several genes, including those conferring resistance. Identification of candidate genes and their resistance functions are crucial to map the hierarchical network of genes involved in the biosynthesis of a given or set of RR metabolite(s), which directly suppresses the pathogen progress. The R genes biosynthesizing these RR metabolites are regulated by other R genes; thus the functional hierarchy must be confirmed before an individual gene can be transferred to a susceptible genotype [12]. Also these specific candidate genes can be transferred without linkage drag effects based on genome editing tools. In view of this, an attempt was made to dissect one of the major FHB resistant QTL-Fhb2 to identify the underlying genetic controls. Integrating two systems biology disciplines, metabolomics and transcriptomics, we were able to identify putative genes localized within the QTL-Fhb2 region and their plausible mechanisms of resistance.

Metabolic profiling identified potential FHB resistance biomarkers
Non-targeted metabolic profiling of genotypes with contrasting levels of FHB resistance identified several RR metabolites, highly accumulated, in different metabolic pathways such as glycerophospholipid, phenylpropanoid (HCAAs), flavonoid, fatty acid and terpenoid. Glycerophospholipids like phosphatidic acid (PA) (FC = 17.72), phosphatidylcholine (FC = 6.56), phosphatidylinositol (FC = 5.77) are known to be deposited to enforce cell walls. Furthermore, these are also known to perceive and transmit signals activating downstream genes eventually regulating R genes to biosynthesize RR metabolites and RR proteins [43]. These compounds are either converted into bioactive lipids (components of lipid bilayer of cell membrane) or stay as soluble molecules (messengers), further binding to the downstream mitogen activated protein kinases (MAPK), further affecting the enzymatic activities, vesicle trafficking and ion fluxes [43]. PA acid is an important secondary messenger in plant stress signaling [44]. These stresses involve pathogen attack, salinity, cold, drought, heat and wounding. In regard to pathogen response, PA acid has been shown to accumulate in response to various elicitors such as xylanase, flagellin, and chitosan [45]. Interestingly, a gene diacylglycerol kinase (DGK) that catalyzes the conversion of structural lipids (PC, PE, PS) into PA was upregulated in our study. This clearly illustrates that the early membrane modifications and their involvement in further activating defense responses are very crucial. Rachis resistance in wheat is mainly due to the deposition of HCCAs [6]. In our study, we detected a higher abundance of two HCCAs in particular, N-Caffeoylputrescine (FC = 5.03) and Feruloyl-2-hydroxyputrescine (FC = 3.31) upon F. graminearum inoculation. Suberin, a complex, intractable biopolymer, with polyaromatic domain of HCCAs, is deposited apoplastically between the primary cell wall and plasma membrane to enforce cell walls [46]. Once, the cell walls are thickened the pathogen won't be able to spread within the spike, through the rachis, thus imparting high levels of rachis resistance. Resistance in NIL carrying resistant alleles of QTL-Fhb1 was reported to be due to the deposition coumaroylputrescine, feruloylputrescine, coumaroylagmatine, cinnamoylserotonin, feruloylagmatine, p-coumaroylserotonin [6]. HCCAs such as feruloylputrescine, p-coumaroyltyramine, N-feruloyltyramine, 4-coumaroyl-3-hydroxyagmatine, feruloylagmatine, 4-coumaroylagmatine, terrestriamide, and feruloylserotonin were reported to impart late blight resistance in potato [30,33]. HCCAs imparted resistance in tomato against Pseudomonas syringae [47], in maize against F. graminearum [48], in onion against Botrytis allii [49]. These HCCAs are not only involved in cell wall thickening, but also they possess antioxidant, antiviral, antibacterial and antifungal activities [47]. Silencing of St-WRKY1 reduced not only HCAAs, but also increased susceptibility of potato to late blight, thus confirming the role of HCAAs in resistance. These HCCAs can be used as markers to screen FHB resistance genotypes, because RR metabolites being the end products of R gene function, guarantee resistance [14].
Apart from HCCAs, we also have identified several flavonoids with high abundance, as RRI and RRC metabolites, such as,  methylquercetin. The growth of Fusarium and macroconidia formation is completely inhibited by dihydroxyquercetin [50], and of the fungus Neurospora crassa by quercetin 3-methyl ether and its conjugated glucosides [51]. The deposition of flavonoid conjugates (of glucoside and methoxy) was higher in rachis tissues of NIL carrying resistance alleles of genes in QTL-Fhb1 [6]. The involvement of both preformed and induced flavonoids in plant defense against pathogens, herbivores, and environmental stress is well documented [52]. In resistant barley several flavonoids were accumulated in high abundance upon F. graminearum inoculation [35,38].
We detected high fold accumulation of preformed and induced free fatty acids such as dihydroxyoctadecenoic acid, 2,3-bis(trimethylsilyl)oxy-butanedioic acid bis (trimethylsilyl) ester, cucurbic acid, 2,4-dimethyl-2-eicosenoic acid, N-palmitoyl phenylalanine, and N-arachidonoyl alanine. Fatty acids are not only part of structural constituents, but also act as secondary messengers and regulators of signal transducing molecules or transcription factors [53]. Arachidonic acid acts as an elicitor in plant defense response to phytopathogens [54]. Several free fatty acids were accumulated in barley upon F. graminearum invasion [38]. The antifungal capabilities of octadecenoic, tetradecanoic, docosanoic, butenoic acid have been reported [55]. Hence, these fatty acids may be crucial components of FHB resistance in wheat, not only acting as physical barriers, but also as antimicrobials.

Transcriptome changes provided key insights to genetic reprogramming upon pathogen invasion
Transcripts belonging to phenylpropanoid and flavonoid pathways, including receptor kinases, transcription factors, detoxification, and signaling genes were highly regulated, following pathogen inoculation ( Table 2). The elicitors produced by pathogen are perceived by plant membrane receptors. In our study, we found higher transcript abundances of lectin receptor kinase (LRK), serine threonine-protein kinase (STPK), proline-rich receptor-like protein kinase (PERK1), and wall-associated receptor kinase 3 (WAK3). The role of LRKs [56], STPKs [57], and WAKs [58] in plant defense is well documented. The overexpression of WAK1 in Arabidopsis thaliana conferred higher resistance to Botrytis cinerea [58]. The resistance to fungal pathogen Magnaporthe grisea was due to a G-type lectin receptor kinase (Pi-d2) in rice [59]. These receptor kinases transduce signals downstream, activating several groups of TFs. The TFs belonging to bHLH, WRKY and MYB groups were up-regulated, depicting their involvement in regulating downstream R genes that biosynthesize RR metabolites and proteins. In our study, the phenylpropanoid pathway genes, such as, agmatine coumaroyltransferase (ACT), caffeic acid 3-o-methyltransferase (CoMT), laccase, phenylalanine ammonia-lyase (PAL), 4-coumarate CoA ligase (4CL), and cinnamyl alcohol dehydrogenase (CAD) were highly expressed in R-RIL. PAL, a hub gene, that biosynthesizes precursor for phenylpropanoid biosynthetic pathway, was highly up-regulated in Sumai-3 upon F. graminearum invasion [41]. ACT which is localized within wheat FHB resistant QTL-2DL imparts high rachis resistance by cell wall thickening [60]. 4CL is an important R gene for both lignin and flavonoid biosynthesis, and was induced in cucumber against powdery mildew [61], cotton against wilt fungus Verticillium dahlia [62], and potato against Phytophthora infestans [30]. Laccase and peroxidase (POD) involved in lignin biosynthesis were up-regulated in our study, emphasizing an increased lignification of cell walls around infection site in R-RIL. In our study, the peroxidase was highly expressed in RP (FPKM = 45.08). The involvement of POD in the defense responses of wheat to Fg infection has been reported [39]. Genes involved in flavonoid biosynthesis like chalcone synthase8 (CHS8), cinnamoyl reductase (CR) and dihydroflavonol 4-reductase (DHFR) were detected only in R-RIL. The disruption of flavonoid pathway significantly reduced flavonoid metabolites [61]. Resistance in wheat to the hemibiotrophic pathogen, Septoria tritici was due to a higher accumulation of CHS transcripts [63]. RR proteins restrict the spread of pathogens. In our study, we detected several PR proteins such as chitinase, peroxidase, and PR1 with higher expressions. Chitinases are very important in plant defense against many fungal pathogens, as they degrade the fungal cell walls, which are primarily made up of chitin. Expression of barley class-II chitinase gene in wheat conferred high level of resistance against F. graminearum under greenhouse and field conditions [64]. Expression of rice chitinase enhanced resistance against Magnaporthe grisea in rice [65], Uncinula necator in Italian ryegrass [66], and Puccinia coronata in grapevine [67].
Mycotoxins produced by Fg, such as trichothecenes, play a major role in pathogenesis, especially DON, a well-known virulence factor. Therefore, the resistance to DON is crucial to confer enhanced FHB resistance [68]. Mutant Fg strains (unable to produce DON) showed reduced FHB severity [69]. In previous studies, it has been reported that DON is converted into less toxic DON-3-O-glucoside (D3G) [70,71]. However, several genes are involved in DON reduction in plants, such as multidrug-resistant protein, multidrug resistance-associated protein, UDP-glycosyltransferase and ABC transporters were detected with higher transcript abundances upon Fg inoculation in wheat cv. Nobeokabouzu-komugi [17]. The accumulation of glutathione S-transferase (TaGSTF5) in wheat resistant cv. Ning7840 upon Fg invasion was significantly higher [72]. Similarly, in our study we found higher expressions of ABC transporter b family member 4, UDP-glycosyltransferase 85a2, UDP-glycosyltransferase 74e1, pleiotropic drug resistance protein 4 and glutathione S-transferase in R-RIL upon pathogen inoculation. Transgenic Arabidopsis and wheat expressing a barley UDP-glucosyltransferase (HvUGT13248) detoxifies deoxynivalenol and provides high levels of resistance to F. graminearum [68,73]. ABC transporter proteins (yeast PDR5 like) confined to plasma membrane confers partial resistance against trichothecenes in wheat by serving as drug efflux pumps [74]. The higher transcript abundances of detoxification genes, clearly explain the reduced levels of DON accumulation in R-RIL, thus contributing to FHB resistance.

QTL-Fhb2 imparts resistance through additive effects of cell wall reinforcement and DON detoxification
The markers flanking the QTL locus were sequenced and the sequence (http://wheat-urgi. versailles.inra.fr/) within the two flanking markers was retrieved and the potential R genes in that region were identified (Fig 4). QTL-Fhb2 was consistently mapped on chromosome 6BS, conferring high rachis resistance [23,24]. Based on high FC metabolites and transcripts, we identified 4CL, CS, bHLH041, GST, ABC4, and CAD as putative candidate R genes localized within the QTL-Fhb2 region. Based on our study, we propose a hypothetical model for FHB resistance in wheat line carrying resistant alleles of genes in QTL-Fhb2 (Fig 6). The importance of each candidate gene in the model on plant defense is discussed.
4CL, CS and CAD confer resistance through cell wall reinforcement 4-Coumarate CoA: Ligase (4CL). We detected a higher transcript abundance (FC = 1.23) of 4CL in R-RIL. 4CL is an important enzyme that catalyzes the conversion of coumaric, ferulic, caffeic, and sinapic acids into hydroxycinnamoyl-CoA thiol esters, which further serve as precursors in lignin, flavonoid, polyphenol, coumarin and suberin biosynthesis [75]. Two HCCAs, N-Caffeoylputrescine (FC = 5.03) and feruloyl-2-hydroxyputrescine (FC = 3.31) were in high abundance, showing a higher expression of 4CL to biosynthesize them in R-RIL than in S-NIL. Increased expression of 4CL increased HCCAs accumulation in potato against P. infestans, imparting resistance through thickening of cell walls [28][29][30]33]. Interestingly, we also detected higher abundances of metabolites and transcripts from lignin and flavonoid pathways, implicating the role of 4CL in the biosynthesis of precursors. Silencing of 4CL reduced lignin content in switch grass [76], and poplar [77]. As 4CL is a hub R gene in the phenylpropanoid pathway leading to the biosynthesis of lignins and flavonoids, we consider this to be a potential candidate R gene conferring high rachis resistance to FHB.
Callose synthase (CS). Plants restrict the spread of pathogens through deposition of an RR metabolite, such as callose (β-1,3-glucan) to form cell wall appositions or papillae [78]. Papillae are complex structures formed around the invading hyphae, which are produced inbetween plasma membrane and the cell wall. The papilla is composed of different classes of compounds such as, phenolics, reactive oxygen species, cell wall proteins and glucans [79]. In our study, we also detected a higher transcript abundance of callose synthase. Callose synthase 5 (CalS5) in Arabidopsis thaliana plays a predominant role in the synthesis of the callose wall and callose plugs, and containment of powdery mildew hyphae in Arabidopsis [80]. Arabidopsis thaliana callose synthase PMR4 expression in barley increased penetration resistance to powdery mildew [81]. Hence, we consider callose synthase as one of the integral candidates in FHB defense.
Cinnamyl alcohol dehydrogenase (CAD). CAD is a key enzyme in lignin biosynthesis that catalyzes reduction of cinnamaldehydes into cinnamyl alcohols, the last step of monolignol biosynthesis, before oxidative polymerization in the cell wall [82]. Lignin is a complex phenolic polymer which is deposited in the cell walls of many plants. Deposition of lignin (lignification) is known to confer resistance against invading pathogens. Lignification thickens the cell wall and makes it difficult for fungal appressoria to penetrate into the cell [83]. It also makes cell walls more water resistant, and in turn, less accessible to cell wall degrading enzymes [83]. Gene expression profiling and silencing showed monolignol biosynthesis is very important in penetration defense in wheat against powdery mildew invasion [84]. Cinnamyl alcohol dehydrogenase-C and D play a crucial role conferring resistance in Arabidopsis against bacterial pathogen Pseudomonas syringae pv. tomato [85]. AtCAD1 is involved in lignification of elongating stems in Arabidopsis thaliana [86]. CAD was upregulated in NILs containing QTL-Fhb1, upon pathogen invasion [6].
ABC transporter and GST aiding resistance through DON reduction ABC transporter-4 (ABC-4). DON inhibits eukaryotic protein synthesis and increases the virulence of F. graminearum by suppressing RR protein and metabolite biosynthesis in plants. A wheat ABC transporter (TaABCC3.1) imparts DON tolerance [87]. A TaABCC (ABC transporter C family) gene within FHB resistant QTL-2DS conferred resistance by reducing DON accumulation [88]. TaABCC gene underlying wheat resistance QTL-Fhb1 imparts FHB resistance [89]. Similarly, in our study we identified higher transcript abundance of ABC transporter b-family member 4 and consider this to play a significant role in rachis resistance, by reducing DON for virulence.
Glutathione S-transferase (GST). GSTs play an important role in plant resistance against biotic and abiotic stresses [90]. These are dimeric enzymes that catalyze the conjugation of electrophilic molecules to glutathione (GSH) [91]. DON the Fusarium virulence factor upon conjugation with GST is presumably detoxified, thus decreasing the pathogenicity of the fungus [27,92]. GSTs from barley are reported to detoxify DON [70]. Glutathione S-transferase genes were induced in Nicotiana benthamiana upon Colletotrichum destructivum and C. orbiculare infection and a GST gene was implicated in resistance [91]. GSTs were induced in potato upon Phytophthora infestans, wheat upon Erysiphe graminis and Arabidopsis upon Peronospora parasitica infection [93]. The expression of GST gene from Lilium regale conferred resistance to Fusarium oxysporum in tobacco [94]. Similarly, we presume that the GST identified in our study, reduces the virulence of pathogen by detoxification of DON, and in turn, aid in enhancing FHB resistance.
bHLH041: A novel candidate for FHB resistance. The bHLH group of transcription factors shares a basic helix loop helix protein structure. The size of these transcription regulators varies anywhere from 60-100 amino acids and consists of two highly conserved domains [95]. The N-terminal basic domain functions as a DNA-binding domain, and the second basic domain which is separated by a loop, determines the dimerization capacity of a protein [95,96]. These transcription regulators usually bind to a consensus sequence known an Ebox (CANNTG) [97]. The role of bHLH transcription regulators in plant response to pathogen attack has been well documented [96,98]. A transcription regulator TabHLH060 was highly expressed in wheat leaves upon invasion by an obligate pathogen Puccinia striiformis f. sp. tritici [96]. In our study, the transcription regulator TabHLH041 was detected only in RP (FPKM = 0.42). Further functional characterization and the identification of its potential downstream targets should increase our understanding about its role in FHB defense reactions.
Our study identified several important R genes localized in the QTL-Fhb2. Even though this QTL was fine mapped it contains several genes with different mechanisms of resistance, but acting cumulatively to impart high level of rachis resistance. These genes should be sequenced to verify if they are functional in R-RIL, but nonfunctional in S-RIL. These polymorphic candidate genes identified here can be used in breeding, following validation of gene resistance functions based on silencing these genes in resistant genotype. The functional alleles of these genes can be used to replace the alleles in susceptible commercial cultivars, if nonfunctional, based on genome editing.

Conclusion
FHB resistance is polygenic in nature; many R genes with major and/or minor effects contribute cumulatively conferring resistance. The mapped FHB resistant QTLs localize several genes governing resistance. We dissected FHB resistant QTL-Fhb2, based on integrated metabolotranscriptomics approach and putatively identified six resistance genes. These genes confer structural resistance through cell wall reinforcement, reducing the spread of pathogen through rachis within the spike. Several other genes detoxify DON, the virulence factor, and thus reducing the severity of the disease. In conclusion, we report that the wheat resistant QTL-Fhb2 confers high rachis resistance through combined effects of cell wall enforcement and reduction of DON. The candidate genes identified here, upon functional characterization, can pave the way to develop highly FHB resistant genotypes.

Plant and pathogen production, and inoculation
All experiments were conducted in greenhouse as a randomized complete block design, with two RILs (R-RIL and S-RIL), two inoculations (pathogen and mock as control) and three to five replications over time, depending upon the nature of the experiment. In each pot, 4 seeds were sown at 5 d intervals. Plants were grown at 25 ± 3°C, 70 ± 10% relative humidity and 16 hours of photoperiod throughout the growing period. Plants were irrigated every day and fertilized at 15 d intervals with 20-20-20 NPK, and the trace elements according to the growth stage of the plants [6]. The F. graminearum, isolate Z-3639 was initially cultured on potato dextrose agar (PDA, DIFCO Laboratories Detroit, Michigan, USA) medium to produce mycelia and then in rye-B agar medium with incubation under UV light for the production of macroconidia [99]. Cultures grown for a week were used to prepare macroconidial suspension for inoculations. The final concentration of spore in the suspension was adjusted to 10 5 ml -1 using sterilized distilled water and 10 μl suspension was inoculated per spikelet [20]. Plants were inoculated at about 50% anthesis (GS = 60-65) and were covered with plastic bags for 48 hours to maintain high humidity in order to facilitate initial infection.

Phenotyping of RILs
Disease severity analysis and fungal biomass quantification were carried out to phenotype RILs, for FHB resistance. For disease severity analysis, a single alternate pair of spikelets, in the middle of a spike of each RIL was inoculated with F. graminearum macroconidial suspension. The total number of spikes inoculated per replication were ten and the total number of replications were five. The number of spikelets diseased were recorded every three days until 21 days post inoculation (dpi) and the proportion of spikelets diseased (PSD) and the area under disease progress curve (AUDPC) was calculated. For fungal biomass quantification, three alternate pairs of spikelets per spike per RIL were inoculated with Fg spore suspension and distilled water, separately. The number of spikes inoculated per replication and the number of total replications were three. 7dpi, six successive pairs of spikelets (three alternate pairs inoculated and three uninoculated) were harvested; the rachis region underlying these six spikelets was collected. The total DNA was extracted from rachis tissues using DNeasy Plant Mini Kit (Qiagen, Germany). The concentration of DNA was normalized to 20 ng/μl and the fungal biomass was quantified using F. graminearum specific Tri6 as target gene and wheat specific Actin as a housekeeping gene [20]. The relative gene copy number of Tri6 was calculated following 2 −ΔΔC T method [100].

Sample collection, metabolite extraction, metabolic profiling, and data analysis
For metabolic profiling, the experiment was conducted in a randomized complete block, with two genotypes (R-RIL and S-RIL), two inoculations (pathogen = P and mock = M) and five replications per treatment. Three alternate pairs of spikelets of a spike were inoculated with F. graminearum spore suspension (P) or distilled water (M). 7dpi, six successive pairs of spikelets (three alternate pairs inoculated and three uninoculated) were harvested and the rachis region underlying these six spikelets was separately collected, ground in liquid nitrogen using prechilled mortar and pestle. 100 mg of tissue samples were weighed in 2 ml sterilized micro-centrifuge tubes and used for metabolite extraction. Metabolites were extracted using absolute methanol followed by 60% methanol in order to extract polar, semi-polar and non-polar metabolites [101]. The metabolites were analyzed using LC-ESI-LTQ-Orbitrap-MS. Randomization of samples was done to avoid any structural errors associated with the liquid chromatography and high resolution mass spectrometry (LC-HRMS = LC-LTQ-Orbitrap) analysis. The output data files obtained from LC-MS analysis were first converted into mzXML/.cdf and were exported to MZmine2 software for peak deconvolution, peak detection, spectral filtering and normalization of peaks [102]. The abundance of peaks were subjected to paired t-test (comparison of two treatments at a time) to identify treatment significant metabolites. Treatment significant metabolites with P < 0.05 were retained for further analysis. Metabolites with higher abundance in resistant genotype (R-RIL) than in susceptible genotype (S-RIL) were considered as resistance related (RR) metabolites. A RR metabolite based on inoculation (M) was considered as constitutive resistance related constitutive (RRC = RM>SM) metabolite. A metabolite with significantly higher abundance in the pathogen inoculated treatments than in mock inoculated treatments was considered as a pathogenesis related (PR) metabolite, in resistant (PRr = RP>RM) or susceptible (PRs = SP>SM) genotypes. A PRr metabolite in a resistant genotype with abundance greater than that in susceptible pathogen (PRs) inoculated was considered as a resistance related induced metabolite (RRI = (RP>RM) > (SP>SM)) metabolite [6,10]. The resistance metabolites were identified with putative compound names using different databases PlantCyc (http://www.plantcyc.org/), METabolite LINk (METLIN) (https://metlin.scripps.edu), KNAp-SAcK (http://kanaya.naist.jp/KNApSAcK/) and Kyoto encyclopedia genes and genomes (KEGG) (http://www.genome.jp/kegg/). The putatively identified metabolites were further confirmed based on: i) accurate mass match (with accurate mass error < 5 ppm) [6,35]; ii) fragmentation pattern match [35,38].
Sample collection, RNA extraction, library preparation, Illumina sequencing, and data analysis For transcriptome analysis, three alternate pairs of spikelets of three spikes per RIL were inoculated with F. graminearum spore suspension and water (control). Each treatment (resistant pathogen = RP, resistant mock = RM, susceptible pathogen = SP, susceptible mock = SM) consisted of three biological replicates. At 48 hpi, the three inoculated and three un-inoculated spikes were harvested, the rachis in inoculated region was harvested, and ground in liquid nitrogen using pre-chilled mortar and pestle. 100 mg of tissue samples were weighed in 2 ml sterilized micro-centrifuge tubes and the total RNA was extracted using Qiagen RNeasy Mini Kit (Qiagen, Germany). Total RNA was quantified using a NanoDrop Spectrophotometer ND-1000 (NanoDrop Technologies, Inc.) and its integrity was assessed using a 2100 Bioanalyzer (Agilent Technologies). Libraries were prepared from 250 ng of total RNA using the TruSeq stranded mRNA sample preparation kit (http://www.illumina.com/products/truseq_stranded_ mrna_sample_prep_kit.html), as per the manufacturer's recommendation. Using the Poly-A selection, mRNA molecules were separated and fragmented, followed by cDNA synthesis, ligation of adapters and cDNA fragments enrichment (PCR) (http://www.illumina.com/ applications/sequencing/rna/mrna-seq.html). Libraries were quantified using the Quant-iT™ PicoGreen 1 dsDNA Assay Kit (Life Technologies) and the Kapa Illumina GA with Revised Primers-SYBR Fast Universal kit (D-Mark). Average size fragment was determined using a 2100 Bioanalyzer (Agilent Technologies). Libraries were sequenced with an Illumina Hiseq 2000 sequencer, with 100 bp paired-end reads (http://support.illumina.com/sequencing/ sequencing_instruments/hiseq_2000.html). The Illumina Hiseq paired end raw reads were quality checked using FastQC. The raw reads were processed using ABLT in-house program to filter adapters and low quality bases towards 3'-end. The raw reads obtained were aligned to the reference genome using TopHat [103]. Cufflinks package was used to assemble transcript, estimates their abundances, and tests for differential expression and regulation [104]. The Cuffdiff program in Cufflinks was used to identify differentially expressed transcripts. The transcripts that were not aligned to reference genome were assembled using de novo assembly [39]. The transcripts were differentially classified as up-regulated, down-regulated and neutrally regulated. The transcripts with FC = ±1 were considered as neutrally regulated transcripts. Two treatments were compared at a time: RP vs RM, and SP vs SM, respectively, and further the fold change (FC) was calculated taking into consideration FC coming from resistant and susceptible genotypes to identify up, down, neutrally regulated expression (RP/RM = FC1, SP/ SM = FC2, FC = FC1/FC2). The transcripts that were only expressed in one of those treatments compared at a time were provided with FPKM values, and the FPKM values were compared in resistant and susceptible to denote significantly higher expression. To make it clear the transcripts that were detected up-regulated only in resistant or susceptible genotype may be down-regulated or neutrally in other, or may be not significant at P < 0.01, so are not detected in the contrasting genotype in RNA-seq data.

Bioinformatics annotation and methods
The transcripts were annotated using Blast2go program against several databases such as nonredundant protein database, Swiss-Prot, KEGG, etc. for annotation and GO functional classification [39]. The transcripts with P < 0.01 and log 2 FC ! 1 were further retained. All the transcripts that were set to threshold of P < 0.01 were crosschecked with the preliminary annotation for all treatments. Pathway analysis was done by using KAAS (http://www.genome.jp/ tools/kaas/). Arabidopsis thaliana (thale cress) and Oryza sativa japonica (Japanese rice) (RefSeq) were taken as reference. The transcription factors (TFs) were predicted by PlantTFDB webserver [105]. We used the threshold determination as set by the server and the established criteria to identify the candidate, which were considered for predicted TFs and the binding sites. Furthermore, the predictions were used for interpreting the GO annotations.

Expression analysis using quantitative Real-Time PCR (qRT-PCR)
The RNA extracted for transcriptome analysis (RNA-seq) was used for cDNA synthesis. The first strand was synthesized using Affinity Script qPCR cDNA Synthesis Kit (Agilent Technologies, Canada). The qRT-PCR was performed in a reaction volume of 10 μl consisting of 25 ng cDNA, 2 pmole of each primer, Qi-SYBR Green supermix (BioRad, Canada). The reactions were carried out in CF"X384TM Real-Time system (BioRad, ON, Canada). The relative transcript abundance in pathogen treated treatments compared to water (control) treatments were analyzed using 2 −ΔΔC T (C T = cycle threshold) [100]. The transcripts, only those were expressed in pathogen treated treatments; RP>SP were analyzed. The relative transcript abundance was represented as FC values; whereas, in RNA-seq data the values are presented as log 2 FC. The primers used for qRT-PCR analysis were designed using NCBI primer blast software (S4 Table) [106].  Table. High fold change RRI and RRC metabolites detected and putatively identified in rachis tissues of wheat RIL carrying resistant alleles of QTL-Fhb2 following F. graminearum and mock inoculation. Table. The transcripts with higher expression values in R-RIL and S-RIL following F. graminearum and mock inoculation at 48 hpi. (XLSX) S3 Table. The list of all the genes localized within the QTL-Fhb2 region based on the available survey sequence available. The genes within two markers GWM-133 and GWM-644 flanking the QTL-Fhb2 region were taken into account. The list of genes with gene ID, chromosome localization and gene ontology is given. Table. List of primers used for qRT-PCR analysis.