Integrated proteomics, genomics, metabolomics approaches reveal oxalic acid as pathogenicity factor in Tilletia indica inciting Karnal bunt disease of wheat

Tilletia indica incites Karnal bunt (KB) disease in wheat. To date, no KB resistant wheat cultivar could be developed due to non-availability of potential biomarkers related to pathogenicity/virulence for screening of resistant wheat genotypes. The present study was carried out to compare the proteomes of T. indica highly (TiK) and low (TiP) virulent isolates. Twenty one protein spots consistently observed as up-regulated/differential in the TiK proteome were selected for identification by MALDI-TOF/TOF. Identified sequences showed homology with fungal proteins playing essential role in plant infection and pathogen survival, including stress response, adhesion, fungal penetration, invasion, colonization, degradation of host cell wall, signal transduction pathway. These results were integrated with T. indica genome sequence for identification of homologs of candidate pathogenicity/virulence related proteins. Protein identified in TiK isolate as malate dehydrogenase that converts malate to oxaloacetate which is precursor of oxalic acid. Oxalic acid is key pathogenicity factor in phytopathogenic fungi. These results were validated by GC-MS based metabolic profiling of T. indica isolates indicating that oxalic acid was exclusively identified in TiK isolate. Thus, integrated omics approaches leads to identification of pathogenicity/virulence factor(s) that would provide insights into pathogenic mechanisms of fungi and aid in devising effective disease management strategies.

factors. Plant pathogenic fungi employ several pathogenicity/virulence related proteins to infect their host plants. They are classified as components of the signal transduction pathways, enzymes for degrading host defenses, transporters to defend the fungus against host defences, toxins effective against host defences and penetration effectors viz. hydrophobins, glycerol and melanin 5 . The pathogenicity/virulence related proteins may serve as potential biomarkers for screening of resistant wheat genotypes and diagnosis of KB pathogen.
In this context, it is pertinent to differentiate the T. indica isolates based on the degree of virulence/aggressiveness which in turn can be employed for the plant breeding programme. In the present study, the virulence of different T. indica isolates was compared on the basis of their pathogenic variability tested on a set of host differentials of Triticum aestivum, in order to select isolates exhibiting contrasting virulence behaviour for comparative proteomic analysis of mycelial proteins for identification of pathogenicity or virulence related proteins.
In the field of fungal proteomics, 2-DE coupled to tandem mass spectrometry has been extensively used for studying fungal plant pathogens with respect to their proteome maps, pathogenicity proteins and virulence factors [6][7][8] . With this view, the present study was aimed to carry out the comparative proteomic analysis of mycelial proteins from T. indica isolates exhibiting varied virulence behaviour. The results obtained by comprehensive proteomic analysis would be further integrated with the de novo genome sequence generated employing the hybrid approach of Illumina HiSeq. 2000 and PacBio single molecule real time (SMRT) sequencing technology 9 . The identified homologs of candidate pathogenicity/virulence related proteins from T. indica genome were further subjected to sequence and structure based functional annotation. The metabolic/gene regulatory networks involved in conferring virulence behaviour can also be analysed by metabolite profiling through metabolomics. Most of the metabolomic studies conducted till date monitor the metabolomic changes in host plant upon infection by the fungal pathogens 10,11 . Certainly, there is no information available regarding the metabolomic study on this important plant pathogenic fungus. To monitor and identify the difference between the T. indica isolates exhibiting varying virulence levels, metabolomes were carried out that would further complement our results obtained by proteomics and genomics studies.
The prime objective of the present study is to integrate the comparative omics data of T. indica isolates showing varied virulence behaviour and the interpretation of the results would facilitate the identification of pathogenicity or virulence related proteins of this economically important fungus.

Materials and Methods
Fungal isolates, plant material and growth conditions. Karnal bunt infected wheat grain samples were collected from Indian Institute of Wheat and Barley Research (IIWR), Karnal. Cultures of T. indica isolates were grown from a single teliospore using technique described by Warham 12 and further multiplied on modified potato dextrose agar (PDA) in biological oxygen demand (BOD) incubator under alternating dark and light conditions at 20 ± 2 °C for 15-20 days. Fungal isolates were maintained on both liquid and solid modified potato dextrose medium. Ten wheat varieties showing differential disease response to T. indica isolates were generously provided by Dr. M.S. Saharan from IIWR Karnal. Wheat seeds were planted in Crop Research Centre, G.B. Pant University of Agriculture and Technology, Pantnagar for pathogenicity testing studies.
Pathogenic Variability Test. Pathogenic variability of ten T. indica isolates was evaluated on a set of ten host differentials of Triticum aestivum for two consecutive years. All differentials of wheat were sown in a meter row with 25 cm line -to -line spacing in mid-November (the normal wheat sowing time). Need based watering and doses of fertilizer were applied as recommended by Tandon and Sethi 13 and each host was maintained in triplicates About 10 4 /ml sporidial suspension from each isolate was used to inoculate the five tillers of each host differential at boot leaf stage 14 using hypodermic syringe according to the procedure described by Aujla et al. 15 . High-humidity was maintained by using a mist sprayer for at least 4 hours a day to ensure infection. Un-inoculated ear heads for each differential were used as negative control. Bunted ear heads were harvested and manually threshed after maturity. Infected grains were scanned visually and categorized into five different grades of infection based on the conversion of the endosperm region into teliospores sooty mass. Percent coefficient of infection was calculated following Aujla et al. 16 . Based on percent coefficient of infection (CI) and percent incidence for two successive crop seasons, T. indica isolates were categorized as highly aggressive, moderately aggressive and low aggressive. CI data was subjected to two way analysis of variance (ANOVA) using STPR 2.0 software.

Tandem Mass Spectrometry and Database Searching.
For in gel-digestion of proteins, the protein spots were manually excised from CBB stained 2-DE gels and suspended in 10% glacial acetic acid. Excised spots were destained with 50 mM NH 4 HCO 3 in 50% (v/v) methanol at 40 °C for 1 hour. After being completely dried in vacuum centrifuge, gel particles were digested with 5 ng/µl of trypsin at 37 °C for 16 h. The digested peptides were extracted using 0.1% trifluoroacetic acid (TFA) in 50% acetonitrile. Peptides were resuspended in α-cyano-4-hydroxycinnamic acid in 50% acetonitrile containing 0.1% TFA. Mass spectra was obtained using tandem mass spectrometry (ULTRAFLEX III TOF/TOF, Bruker Daltonics). Mass spectrometry spectra were acquired with 1600 laser shots per spectrum and MS/MS spectra with 2500 laser shots per fragmentation spectrum. 10 strongest peaks from MS spectra were chosen as the precursor ions to attain the MS/MS fragmentation spectra. The Flex analysis software 3.0 (Bruker Daltonics) was used for spectra analysis and generation of peak list files. The peak list files were searched in NCBI nonredundant (NR) database (http://www.ncbi.nlm.nih.gov), version 20160114, 79354501 sequences and 28992349963 residues for "Fungi") through the MASCOT search engine (http://www.matrixscience.com). Search parameters were kept to taxonomy, Fungi; proteolytic enzyme, trypsin; max missed cleavages, 1; fixed modifications, carbamidomethyl (C); variable modifications, oxidation (M); peptide mass tolerance, 50 ppm; fragment mass tolerance, 2 Da. Functional annotation of identified proteins was performed by Uniprot database (www.uniprot.org/). Identified proteins were grouped into distinct categories through COG and GO pathway analysis at NCBI. Subcellular localization of identified proteins was carried out by TargetP 19 .
Identification of homologs of pathogenicity/virulence related proteins from T. indica genome sequence data. Genome of T. indica TiK isolate sequenced 9 was used for identification of homologs of candidate pathogenicity/virulence related proteins using BLASTP search. Identified putative pathogenicity/virulence related proteins from T. indica genome were subjected to further sequence and structure based functional annotation.
Sequence and structure based functional annotation of homologs of pathogenicity/virulence related proteins. Identified putative pathogenicity/virulence related proteins homologs from TiK genome were extensively analyzed for the functions by InterProScan 20 and NCBI conserved domain database (CDD) 21 . Other databases such as ScanProsite 22 , SMART 23 , PANTHER 24 and CATH 25 were also used for sequence based functional analysis of proteins. Three dimensional structures of identified putative pathogenicity/virulence related proteins were predicted using RaptorX software 26 . Predicted protein models were validated by various bioinformatics tools such as RAMPAGE 27 , ProQ server 28 , DALI server 29 and ProFunc 30 .
Protein-Protein interaction network. Protein sequences of putative pathogenicity/virulence related genes retrieved from genome sequence of T. indica TiK isolate were used to generate Protein-Protein interaction (PPI) network by STRING (Search Tool for the Retrieval of Interacting Genes) database version 10.5 31 taking Saccharomyces cerevisiae as a model system. The obtained PPI network was further visualized using Cytoscape version 3.5.1 32 .
Gas Chromatography -Mass Spectrometry (GC-MS). GC-MS was performed for metabolic profiling of T. indica isolates showing extreme virulence behaviour, TiK and TiP, at Advanced Instrumentation Research Facility (AIRF), JNU New Delhi. Each sample (20 mg) was extracted using 1 ml of 70% methanol. Extract was sonicated and centrifuged at 5,000 rpm for 10 minutes. Supernatant was collected and filtered through a 0.45 mm filter. Sample (100 µl) was transferred into a GC vial and dried with a nitrogen gas flow at 60 °C for 5 min. For derivatization, each sample (100 µl) was oximated with 50 µl methoxyamine hydrochloride in pyridine. After incubation at 30 °C for 90 min, samples were silylated with 25 µl of N-methyl-N-trimethylsilyl trifluoroacetamide (MSTFA) at 37 °C for 60 min. GC-MS was performed according to Oh, et al. 33 .

Results
Differential disease scoring on a set of host differentials of Triticum aestivum suggests varied virulence behaviour of T. indica isolates. Pathogenicity of ten isolates of T. indica was evaluated on a set of ten host differentials of wheat for two successive years (Year 2013-14 and 2014-15). The average coefficient of infection and percent incidence calculated from the field experiment is presented in Fig. 1a,b. The pathogenic variability among the isolates is indicated by host-pathogen interaction on differential hosts. The T. indica isolates can be categorized into three distinct aggressive types, namely Karnal Bunt-High Aggressive type (KB-HAg), Karnal Bunt-Moderate Aggressive type (KB-MAg) and Karnal Bunt-Low Aggressive type (KB-LAg). Four isolates belonged to KB-HAg type, four isolates to KB-MAg type and remaining two isolates to KB-LAg type. The KB-HAg comprised the isolates from Karnal (TiK), Jagadhari (TiJg), Rohtak (TiR) and Tarori (TiT). Two isolates of the pathogen that belonged to KB-Lag type were represented from Panipat (TiP) and Sonipat (TiS). The KB-MAg type of T. indica was from Jind (TiJ), Hansi (TiH), Pipli (TiPp) and Nilokheri (TiN).
The highest per cent incidence of disease was produced by TiK isolate (21.2%), followed by TiJg (13.8%) and TiT isolate (12.3%). The TiP isolate showed the lowest per cent incidence of disease i.e. 2.1% (Fig. 1b). The pathogenicity variability test also suggested that the TiK and TiP isolate as the most aggressive and the least aggressive SCIEnTIfIC RePoRTS | (2018) 8:7826 | DOI:10.1038/s41598-018-26257-z type as they exhibited the highest (10.4) and the lowest (1.1) coefficient of infection, respectively on most of the host differentials as compared to other isolates (Fig. 1a). Thus, these isolates exhibiting differential virulence behaviour are used in the present investigation for comparative proteomic analysis for identification of pathogenicity/virulence related proteins. ANOVA exhibited the significant differences in disease severity among the isolates, wheat genotypes and their interactions. The means square for the genotypes was relatively lower (56.8) than for the isolates (353.3). The mean square for the genotype X isolate interaction was very high (59.9) ( Table 1).

Comparative proteomic analysis of T. indica isolates exhibiting varied virulence behaviour.
In proteomic research, 2-DE coupled to MS serve as powerful tools for identification of isolate-specific expression of pathogenicity or virulence related proteins in various plant pathogenic fungi such as Botrytis cinerea, Verticillium dahliae, Fusarium oxysporum f. sp. conglutinans [6][7][8]34 .
The growth rates of TiK and TiP isolate cultured on PDA and PDB medium were compared from 1 to 40 days by calculating the fungal biomass (g/100 ml of culture). The mycelium of both TiK and TiP increased exponentially (logarithmic growth phase) up to 21 days followed by decrease in the rate of growth (stationary growth phase) (Fig. 1d). The fungal biomass of TiK isolate at lag phase (14 th day) was approximately 1.8 g/100 ml of culture which increased upto 5.4 g/100 ml at 21 st day of growth (exponential phase) (Fig. 1d). In TiP isolate, fungal biomass 1.2 g/100 ml at 14 th day which rose to 4.7 g/100 ml of culture at 21 st day growth cycle. After 21 st day, decrease in growth in terms of dry weight of mycelia was observed in both the isolates (Fig. 1d). It is due to Figure 1. Comparison of virulence of ten T. indica isolates based on differential disease scoring on a set of ten host differentials of Triticum aestivum in terms of (a) Coefficient of infection; (b) Percentage incidence. Each value represent the average of three biological replicates; (c) the growth of T. indicia highly virulent TiK and low virulent TiP after 21 days of culture on PDA; (d) Growth kinetics of TiK and TiP isolate in terms of total biomass production (g/100 ml) at different time intervals. A proteome map (pH 3-10) of protein samples from the mycelia of two T. indica isolates, TiK and TiP, varying in their aggressiveness or virulence levels was obtained with comparative proteomic analysis of high quality 2-DE gels ( Fig. 2; Supplementary Fig. S1). By comparing the proteome profiles from these two isolates, various qualitative variations between protein spots were revealed. Only those spots that were consistently expressed in three replicate gels of highly virulent TiK isolate were considered for further analysis. Twenty one protein spots that exhibited more than a 1.5 fold change were statistically significant (p < 0.05) and selected for tandem mass spectrometry analysis. Among them, proteins corresponding to spots 3-15 and 17 showed higher expression in TiK isolate and proteins in spots 1, 2, 16, 18, 19, 20, 21 were found to be differentially expressed with respect to TiP isolate (Fig. 2). The proteins of such spots were subjected to MALDI-MS/MS analysis.

Degree of Freedom
All the differentially expressed or up-regulated proteins were identified by MASCOT search as protein with known functions, except the protein from spot 3, 14 and 21, that were identified as a hypothetical protein. The calculated molecular mass of the identified proteins ranged from approximately 13 to 395 kDa. The protein spots identified are indicated by circles and numbers on the gels (Fig. 2). The putative protein identity, protein score and number of matched peptides are given in Table 2.
The functional annotation of protein spots was performed using Uniprot database. On the basis of their putative functions, the identified proteins were classified into metabolic pathways as determined by the COG classification system. About 50% of the identified proteins belong to carbohydrate transport and metabolism (CTM) classification, followed by signal transduction (ST) proteins (18%). Post-translational modification related proteins (PTM) and inorganic ion transport and metabolism related protein (ITM) comprises 4% each of the total classified proteins. 24% of the proteins belong to an unknown function classification (Fig. 3a). The sub-cellular localization of these proteins were primarily classified as any other location (45%) followed by 41% proteins localized in mitochondria and 14% proteins in the secretary pathway (Fig. 3b). For twenty one proteins varying in their abundance in TiK and TiP isolates were classified based on Gene ontology (GO) analysis into three distinct categories (cellular component, biological process and molecular function). The GO analysis for molecular function revealed that most proteins possess catalytic activity (64.7%) followed by oxidoreductase activity (35.3%) (Fig. 3e). Majority of the protein in the biological process GO analysis were involved in carbohydrate metabolic process (17.6%) and cell communication (13.8%) (Fig. 3d). While protein complexes (50%) was most dominant in cellular component category (Fig. 3c).
Proteins involved in carbohydrate transport and metabolism may be indispensable for stronger pathogenicity. Spot1 corresponds to ATPase that is usually produced through carbohydrate metabolism 35 . ATPase is essential for providing energy for growth, development, differentiation and sporulation. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was identified in 2 spots, spot 15 and 16 ( Fig. 2; Table 2). The function of enzyme GAPDH is well known in the glycolytic cycle. However, the role of this enzyme is also implicated in various cellular processes in microbial pathogens and mammals such as DNA repair, regulation of actin filaments in microfilaments, regulation of calcium release from the ER and the nuclear export of tRNA [36][37][38][39][40] . Spot 2 corresponding to Glycoside hydrolases 45 (GH45) was differentially expressed in TiK isolate mycelium. β-1,3-glucanosyltransferases and Fructose-bisphosphate aldolase (FBA) were identified in spot 13 and 7, respectively ( Fig. 2; Table 2).     Malate dehydrogenase (MDH) identified in spot 10, 12 and 17, was highly abundant in the mycelium of TiK isolate ( Fig. 2; Table 2). Enzyme MDH is a key regulatory enzyme of tricarboxylic acid cycle that catalyses the reversible conversion of oxalacetate and malate. The product of the reaction catalysed by MDH, oxaloacetate is a precursor of oxalic acid 41 which acts as a pathogenicity factor in B. cinerea and Sclerotinia sclerotiorum [42][43][44][45][46][47] .
Proteins related to post translational modification and protein turnover may aid in stronger pathogenicity. HSP70 protein was highly abundant in TiK mycelium (Spot 4) ( Fig. 2; Table 2). Heat shock proteins (HSPs), a highly conserved family of molecular chaperones that play a crucial role in correct folding, stabilization and post translational modification of proteins 48 . The role of some HSP 70 family proteins has been demonstrated under various abiotic stress conditions such as osmotic stress and temperature stress 49,50 . Protein involved in signaling pathway may contribute to stronger pathogenicity. Protein spot 6 correponding to Ste7 exhibited higher expression in mycelial protein of TiK isolate ( Fig. 2; Table 2). In S. cerevisiae, Ste7, a MAPKK is involved in mating as well as filamentous growth and responsible for invasive growth in haploid and pseudohyphal development in diploid cells. In our laboratory, three MAP kinase homologues, namely Pmk1, Fus3 and Kpp2 have been cloned and studied using in silico tools, in order to define their role in fungal pathogenesis of T. indica 51 . In the present study, comparative proteomic analysis of mycelial proteins from isolates showing extreme virulence behaviour, TiK and TiP indicates the existence of Ste7 homolog in T. indica and its possible role in TiK stronger pathogenicity.
Sequence and structure-based functional annotation. Functional annotation of putative pathogenecity/virulence related proteins is essentially required for understanding the biological processes at the molecular level. It is essential for predicting the function of protein at a systems level which is required for designing the predictive disease model 52 . There are several new methods, including hidden Markov model and neural network model-based tools available for the functional annotation of proteins which are more intelligent and efficient than classical homology search methods 53 . However, structure-based functional annotation is more reliable and much efficient for identifying protein's biochemical or enzymatic function 54 . In present study, the identified putative pathogenecity/virulence related proteins designated as TiHsp70, TiSte7, TiFBA, TiMDH, TiGT and TiGAPDH were further subjected to comprehensive sequence and structure-based analysis: domain analysis by NCBI's conserved domain database and Interproscan; three -dimensional protein structure prediction by RaptorX software; validation of predicted protein models by Ramachandran plot, ProQ, ProFunc, DALI server.
TiHsp70. The putative TiHsp70 sequence may act as a molecular chaperone. This protein belongs to heat shock protein HSP 70 superfamily and has two domains-the amino terminus ATPase domain and the carboxyl terminus substrate binding domain. It is ATP dependent ( Table 3). The protein structures of PDB I.D. 5e84A and 5tkyA, corresponding to chaperones, as templates were used for the tertiary structure prediction. The p-value and score of 8.83e-15 and 686, respectively, were estimated for the predicted 3-D model as shown in Fig. 4a.
The statistical analysis for our model suggested that 96.0%, 2.6%, 1.4% of the residues in derived curcin model were in the favoured region, allowed region and outlier region, respectively ( Fig. 5a; Table 4). Hence, around 98.6% of the residues are placed into the favoured and allowed categories, indicating the prediction of a high quality protein model ( Table 4). The RMSD of the model with respect to the template (PDB code: 5e84A) is 0.49 A°, indicating similar functionality (Table 4). Moreover, ProFunc predicted HSP70_3 protein and heat shock protein Hsp70 conserved motifs showing a close resemblance with that of Hsp70 protein and Asp126-Lys128 and Ala222-Leu224 were identified as structural motifs (Table 3). We obtained the similar results on DALI server which also identified its significant similarity with Hsp70 chaperone (Z score = 68.6) ( Table 4). All these analyses strongly suggest that TiHsp70 may acts as a molecular chaperone.
TiSte7. Domain analysis results revealed that putative TiSte7 sequence possess catalytic domains of serine/threonine-specific and tyrosine-specific protein kinases (Table 3), capable of catalyzing the transfer of the γ-phosphoryl group from ATP to serine/threonine (ST) or tyrosine residues on the protein substrates. InterProScan results showed that TiSte7 sequence possesses a protein kinase domain with a catalytic core which is common to both serine/threonine and tyrosine protein kinases (Table 3). Based on these observations, it is suggested that TiSte7 may function as dual specificity MAP kinase kinase. RaptorX used protein structures of PDB I.D. 3ornA, 3zlsA, 3vn9A and 3eqcA, corresponding to mitogen-activated protein kinase, as templates for the  tertiary structure prediction. p-value and score of 3.96e-10 and 278, respectively, were estimated for the predicted 3-D model as shown in Fig. 4b. The statistical analysis for our model suggested that 96.0%, 2.6%, 1.4% of the residues in derived curcin model were in the favoured region, allowed region and outlier region, respectively ( Fig. 5b; Table 4). Hence, around 98.6% of the residues were placed into the favored and allowed categories, indicating the prediction of a high quality protein model ( Table 4). The RMSD of the model with respect to template was 0.67 A°, which is quite good, indicating functional similarity (Table 4). ProFunc analysis identified Serine/Threonine Protein kinase and T-X-Y motif. The structural motifs were identified as Thr354-Val356 and a Phe613-Arg615 (Table 3). DALI results showed a significant similarity with dual specificity mitogen-activated protein kinase with a maximum Z score = 43.6 ( Table 4). These results also suggested a probable Serine/Threonine Protein kinase activity of TiSte7.
TiFBA. Domains of putative TiFBA sequence is conserved with FTBP_aldolase_IIA superfamily that includes fructose-1, 6-bisphosphate (FBP) aldolase (Table 3). InterProScan result also suggested that this sequence contains FBA_II domain and has TIM beta/alpha barrel characteristic of aldolases (Table 3). RaptorX used protein structures with PDB I.D. 1c1f:A and 1dos:A as templates for the prediction of 3-D model. Template 1dos:A is Fructose-bisphosphate aldolase from Escherichia coli and 1c1f:A is sugar binding protein of. score for the predicted model were 7.31e-10 and 368, respectively. The 3-D structure as predicted by RaptorX software is shown in Fig. 4c.
Ramachandran plot for the predicted 3-D protein structure showed that 93.5% and 4.6% of the residues are in the favoured and allowed region 1.9% of the residues are in the outlier region. So, 98.1% of the residues were in the favored and allowed categories ( Fig. 4c; Table 4), indicating the prediction of a high quality protein model. The structure of TiFBA showed a close resemblance with the template, Fructose-bisphosphate aldolase from Escherichia coli with an RMSD of 0.16 A° (Table 4). Functional motif search resulted in the identification of FruBisAldo_II_A and Fructose-bisphosphate aldolase class-II conserved motifs in the TiFBA. Profunc server identified Ile226-Val228, Leu311-Thr313 and Gly325-Ser327 as structurally significant motifs (Table 3). DALI search results were consistent with ProFunc finding and showed a significant similarity with the fructose-bisphosphate aldolase II (Z score = 51.4) ( Table 4). These findings clearly indicated the possible fructose-bisphosphate aldolase activity of TiFBA.

TiGAPDH. Conserved domain analysis showed that it belongs to Gp_dh_N superfamily that has
Glyceraldehyde-3-phosphate dehydrogenase, containing N terminal NAD-binding domain and C terminal catalytic domain (Table 3). Interproscan results suggested that TiGAPDH sequence belongs to Glyceraldehyde/ Erythrose phosphate dehydrogenase family possessing two conserved functional domains, a highly conserved catalytic domain and an NAD-binding domain. Protein Structure with PDB I.D. 4o59O, 4iq8A, 4k9dA, 3e5rA and 2vynD were used as templates for 3-D structure prediction. All the templates correspond to Glyceraldehyde-3-phosphate dehydrogenase. The predicted 3-D model had a p-value and score of 4.78e-15 and 353, respectively. This indicated the prediction of high quality 3-D model as shown in Fig. 4d.
The statistical analysis of Ramachandran plot for the predicted 3-D model suggested that 95.3% of the residues are in the favoured region. Residues in the allowed region are 3.2% while 1.5% of the residues are in the outlier region. In all, 98.5% of the residues were placed into the favored and allowed categories ( Fig. 5d; Table 4). This showed that the 3-D model predicted is of high quality with respect to protein folding. The structure of TiGAPDH showed a profound similarity with the Glyceraldehyde-3-phosphate dehydrogenase from S. cerevisiae with RMSD of 0.39 A° (Table 4). Functional analysis using ProFunc database scan revealed the presence GAPDH-I signature in the TiGAPDH with structural motifs Gly114-Ala116 and Gly226-Leu228 (Table 3). DALI results also showed  (Table 4). All these systematic analyses precisely suggested that TiGAPDH may have dehydrogenase-like activity and may be essential for pathogen survival.
TiGT. Domain analysis suggested the similarity with Glycosyltransferase_ GTB_ type superfamily that has domain for glycosyl transferase activity (Table 3). InterProScan result also suggested that the TiGT sequence possesses a motif that is involved in glycosyl transferase activity. Protein Structure with PDB I.D. 2bisA, 5d00A, 3mboA, 3okcA and 4xsoA, corresponding to transferase, were used as templates. The predicted protein structure had a p-value and score of 2.22e-07 and 288, respectively. The predicted 3-D protein structure is given in Fig. 4e.
For Ramachandran plot, the statistical analysis of the predicted protein model showed that 92.9% of the residues are in the favoured region. Residues in the allowed region are 4.5% while 2.6% of the residues are in the outlier region. In all, 97.4% of the residues were placed into the favored and allowed categories ( Fig. 5e; Table 4). These results indicated that a good quality protein model was predicted by RaptorX. The RMSD of the model with reference to the template (PDB code: 2bisA) is 1.49 A° (Table 4) which indicates the functional similarity. Functional anaylsis revealed the presence of Glycosyl transferase like conserved motif in TiGT with Glu306-Asp309 as structurally significant motif (Table 3). These results were consistent with DALI finding that showed similarity with glycosyl transferase (Z score 37.2) ( Table 4). The results of sequence and structure based functional annotation suggested that TiGT possess glycosyl transferase activity.   TiMDH. CDD result suggested that this protein sequence has a malate dehydrogenase activity with NAD binding site and substrate binding site (Table 3). This is consistent with InterProScan results that also suggested the malate dehydrogenase activity. The templates utilized by RaptorX to predict 3-D structure of TiMDH sequence according to their PDB I.D. are 2dfdA and 1sevA, corresponding to Malate dehydrogenase from human and glyoxysomal malate dehydrogenase, respectively. The predicted structure with a p-value and score of 1.52e-08 and 329, respectively indicated that model is of high quality. The 3-D structure predicted by RaptorX is shown in Fig. 4e. The statistical analysis of Ramachandran plot for our model suggested that 96.1%, 3.3%, 0.6% of the residues in derived curcin model were in the favoured region, allowed region and outlier region, respectively. Thus, altogether 99.4% of the residues were placed into the favored and allowed categories; which indicates that the model structure derived from RaptorX was of higher quality in terms of protein folding ( Fig. 5e; Table 4). The RMSD of the model with respect to the template (PDB code: 2dfdA) is 0.42 A°, indicating similar function (Table 4). Furthermore, ProFunc predicted lactate/malate dehydrogenase, alpha/beta C-terminal conserved domain and malate dehydrogenase active site signature motif showing a close resemblance with that of malate dehydrogenase protein. Gly8-Ala10, Lys221-Gly223 and Lys136-Tyr139 were identified as structural motifs (Table 3). We obtained the similar results on DALI server which also identified its significant similarity with malate dehydrogenase protein (Z score = 54.2) ( Table 4). All the analyses strongly suggest that TiMDH may probably be malate dehydrogenase.

Protein-Protein interaction network analysis of putative pathogenicity/virulence related proteins.
We have constructed an extended PPI network using putative pathogenicity/virulence related proteins from TiK genome through STRING database based on a high confidence score of 0.7. This database provides information about both experimental and predicted interactions from different sources based on their neighborhood, cooccurrence, gene fusions, co-expression, experiments and scientific literature. The constructed network implies interactions that only with high level of confidence score were taken from the database and considered as valid links. The obtained PPI network from STRING database was visualized by Cytoscape (Fig. 6). The network was analyzed topologically based on the network parameters such as node degree and betweenness centrality using Network Analyzer (Table 5). Each gene/protein in a given network is represented as a node whereas the interactions present between the nodes are defined as edges. The degree showed the number of edges linked to any node. The nodes having high degree may know as the hub proteins possessing key biological functions. The betweenness centrality represents the importance of the node, which was based on the number of shortest paths that is passing through each node. The network was visualized and analyzed based on these parameters to map   the node degree as the node size and betweenness as the node color in the visualize parameter of the Network Analyzer. Nodes with high degree were displayed as a big circle and these nodes were considered as the hubs. The nodes TDH1, TDH2, TDH3 (isoforms of Glyceraldehyde-3-phosphate dehydrogenase), FBA1 (fructose-1, 6-bisphosphate aldolase), ENO1, ENO2 (isoforms of enolase), PGI (phosphoglucose isomerase), TPI (Triose phosphate isomerase) are found as hub nodes, which can be utilized as a molecular targets for development of the disease prevention and management strategies against Karnal bunt of wheat.

GC-MS based metabolic profiling of T. indica isolates with differing virulence behaviour.
The chromatographic data of the compounds identified by GC-MS in two T. indica isolates are given in Tables 6 and 7. Interestingly, among saturated fatty acids, oxalic acid was exclusively identified from mycelia of highly virulent TiK isolate (Table 6). Oxalic acid is well characterized pathogenicity factor in plant pathogenic fungi, including B. cinerea and S. sclerotium [42][43][44][45][46][47] . It is synthesized from oxaloacetate which in turn is synthesized from malate in the reaction catalyzed by enzyme malate dehydrogenase (MDH). In the present study, enzyme MDH has been identified from the proteome of mycelial protein from TiK isolate. It was further identified from de novo assembled T. indica genome sequence and validated through sequence and structure based functional analysis. The present study suggests the role of oxalic acid as pathogenicity factor in KB pathogen, T. indica. The levels of steric acid and palmitic acid were higher in TiK samples, indicating more active fatty acid synthesis machinery in highly virulent TiK isolate. Among, unsaturated fatty acids, oleic acid level was higher in TiK sample (Table 6).

Discussion
Selection and use of T. indica resistant wheat varieties is a key to successful KB disease management. The present study showed that the host-pathogen reactions were quite consistent, regardless of the environmental conditions. Various T. indica isolates exhibited differential reactions on a set of differential host genotypes. It also confirms the earlier studies regarding the existence of pathogenic variability in T. indica. It is revealed that the high aggressive type (KB-HAg) was present in TiK while the low aggressive type in TiP isolate. Comparative proteomic analysis of mycelial proteins from TiK and TiP isolates revealed that the proteins associated with carbohydrate metabolism, namely Glyceraldehyde-3-phosphate dehydrogenase, Malate dehydrogenase, Fructose-1,6-bisphosphate aldolase and Glucanosyltransferase showed highest percentage (50%) in the mycelium of virulent TiK isolate. Higher expression of proteins associated with carbohydrate metabolism in mycelia of TiK than TiP is consistent with the fact that TiK has higher growth rate than TiP isolate. The fast growing TiK mycelium requires more energy than TiP.
GAPDH was highly abundant in the TiK mycelium. GAPDH serve as virulence factor in several pathogenic fungi, including Candida albicans, Trichomonas vaginalis, Phycomyces Blakesleeanus, Paracoccidioides brasiliensis and Botrytis cinerea [36][37][38][39][40] . Apart from cytosol, a GAPDH is also present on the surface of pathogenic fungi where it binds to host's fibronectin. This facilitates pathogen to localize at intracellular and extracellular matrix environment and aid in invasion and colonization of the host tissue. In the present study, spot 16 identified as GAPDH   was uniquely present in the mycelia of TiK isolate, indicating that the oxidative metabolism in the highly virulent isolate might be much more active than that of the low virulent isolate. It may also be hypothesized the presence of surface GAPDH in virulent T. indica isolate and its role in facilitating fungal invasion and colonization of host wheat tissues. This suggested a putative role of enzyme GAPDH as virulence factor in T. indica. In TiK, FBA that catalyzes the aldol cleavage of fructose bisphosphate, was highly abundant. In Candida albicans, a mutation in FBA gene of perturbs the fungal growth 55 , FBA present on the surface serve as an important virulence factor in Paracoccidioides, by involving in the process of fungal adhesion, invasion and colonization. The FBA protein was also found to be upregulated in Fusarium oxysporum race 4 displaying stronger virulence 56 . FBA may also be present on the surface of TiK isolate that might play crucial role in T. indica adhesion and invasion and act as important virulence factor. The role of β-1,3-glucanosyltransferases (spot 13) have been implicated in the fungal cell wall biosynthesis, morphogenesis and virulence of human pathogenic fungi such as Candida albicans and Aspergillus fumigates 57 . Caracuel et al. 58 examined the role of putative β-1,3-glucanosyltransferase encoding gas1 gene of vascular wilt pathogen, Fusarium oxysporum. They observed that gas1 deletion mutants of F. oxysporum exhibited a significantly lower growth rates, increased resistance to cell wall-degrading enzymes. Li et al. 34 studied the significance of β-1,3-glucanosyltransferase in virulence of F. oxysporum f. sp.conglutinans on cabbage plants. The F. oxysporum f. sp.conglutinans gas1 deletion mutants showed restricted growth and virulence compared to wild type isolates. These results were consistent with the comparative proteomic analysis of two races of F. oxysporum f. sp. conglutinans differing in pathogenicity as glucanosyltransferase showed higher expression in highly pathogenic Race R2. Our comparative proteomic analysis of TiK and TiP isolate also exhibited higher expression of glucanosyltransferase protein in isolate showing stronger pathogenicity. This indicates the essential role of glucanosyltransferase in KB pathogenesis on wheat.
For plant pathogenic fungi, the plant cell wall is the main barrier that is required to breach for colonizing the host plant tissues. They produce an arsenal of cell wall degrading enzymes that are responsible for degrading cellulose. Among them, Glycoside hydrolases of family GH45 (endo-β-1,4-glucanase activity) are of prime importance that cleaves the internal glycosidic bond leaving β-1,4 glucans from the straight chain polymer of glucose, cellulose. The endo-β-1,4-glucanase CelA in Clavibacter michiganensis subsp. michiganensis is a pathogenicity determinant that is necessary for inducing the bacterial wilt in tomato 59 . So far, no information is available regarding the relationship between the endocellulase and T. indica pathogenicity. In this study, Spot 2 corresponding to GH45 endoglucanases was differentially expressed in TiK isolate mycelium may be responsible for its stronger pathogenicity ( Fig. 2; Table 2).
Under biotic stress, both host and pathogen are exposed to various extreme detrimental stress stimuli such as cell wall degrading enzymes, acidic pH and reactive oxygen species. In both host and pathogen, expression of various HSPs are up-regulated with pathogen HSPs being responsible for infection of the host plant and host HSPs in turn provide defense against the pathogen invasion 60-65 . Studies led by Yi et al. 66 showed the significance of HSP70 family proteins LHS1, KAR2 of Magnaporthe oryzae in fungal pathogenicity. The present study suggested that TiK isolate mycelia HSP70 protein abundance may enhance mycelial tolerance to extreme stress condition (viz. acidic pH, cell wall degrading enzymes and reactive oxygen species) that are generated during pathogen invasion of host tissues and ultimately contribute to fungal virulence.
Malate dehydrogenase abundance was high in TiK mycelium (spot 10, 12 and 17). Fernandez-Acero et al. 7 have suggested the role of MDH as a pathogenicity factor in B. cinerea as oxalic acid secretion generates an acidic environment that create an appropriate ecological niche for the fungal pathogenic activities such as secretion of several virulence factors like cell wall degrading enzymes and phytotoxins 67,68 . There is a correlation between oxalic acid biosynthetic capability and S. sclerotiorum virulence as mutants lacking oxalic acid biosynthesis machinery were found to be non-pathogenic while revertant strains that regained their oxalic acid synthesis ability exhibited normal virulence 43 . Likewise, MDH higher expression in virulent TiK isolate mycelia suggested the relation of oxalate synthesis with enhanced fungal virulence. Further, oxalic acid might play a crucial role in invasion of host plant cell by lowering the pH and providing suitable acidic environment for the activity of arsenal of cell wall degrading enzymes. GC-MS based metabolic profiling identified oxalic acid from TiK mycelium, substantiating the role of oxalic acid as potential pathogenicity factor in T. indica.

Conclusion
Despite quarantine significance of T. indica, there is a little knowledge regarding the molecular mechanisms of pathogenesis employed by this important fungus to cause disease. Moreover, all the methods used to manage the disease have proven futile. In order to develop an effective disease management strategy, it is essential to identify pathogenicity or virulence related proteins. However, so far not much information about pathogenic determinants of this economically important fungus is available and certainly not using Hi-throughput proteomics and metabolic tools and platforms. In the present study, we report the first proteome map and comprehensive comparative proteomic analysis of T. indica isolates differing in their virulence/aggressiveness for identification of putative pathogenicity or virulence related proteins expressed in the highly virulent isolate. The identified pathogenecity/virulence related proteins play crucial role in stress response, degradation of host cell wall, adhesion, penetration, invasion, colonization, activation of signal transduction pathway and morphogenesis. Potential pathogenecity/virulence related proteins were further complemented with T. indica genome sequence from hybrid genome assembly, resulting in identification of orthologs of several candidate pathogenecity/virulence related proteins which were annotated through both sequence and structure -based functional analysis. Further, GC-MS based metabolic profiling of T. indica isolates varying in virulence behaviour validates the role of oxalic acid as potential pathogenicity factor in T. indica. The identified pathogenecity/virulence related proteins may serve as potential biomarkers that would have utility for screening KB resistant wheat cultivars, designing target specific fungicides and development of specific, sensitive and rapid field level on -site KB diagnostics.