Sorghum [Sorghum bicolor (L.) Moench] seeds of three South African cultivars were used, namely Amazi Mhlophe (abbreviated as MHL or M), NS 5511 (referred to as bitter and abbreviated as BTT or B) and NS 5655 (referred to as sweet, abbreviated as SWT or S; Agricol, Pretoria, South Africa). NS 5511 and NS 5655 are both grain sorghum hybrids of the malting class. NS 5511 is a red, tannin type, classified as GH (high levels of condensed tannin, responsible for dark testa), while NS 5655 is classified as GM (no condensed tannins, no dark testa). Both NS 5511 and NS 5655 have a rating of “3” (on a 1–9 scale with 1 being the most resistant) in terms of exhibiting disease resistance against head smut, leaf disease and root rot (Capstone Seeds, Howick, South Africa). In contrast, Amazi-Mhlophe is an open pollinated, white grain sorghum type. The three cultivars are listed in the variety list of the Registrar of Plant Improvement, Department of Agriculture, Forestry and Fisheries, South Africa (www.daff.gov.za/daffweb3/Branches/Agricultural-Production-Health-Food-Safety/Plant-Production/Verietal-Listing).

Seeds were surface-sterilized in a 1.2% sodium hypochlorite solution and rinsed with sterile water before being placed in glass Petri dishes with soaked paper towel and incubated at 28°C in the dark for 48 h to germinate. The seedlings were then planted in horticultural-grade vermiculite under fluorescent lights with a 12 h light (≈85 μmol m⁻² s⁻¹) and 12 h dark cycle. The temperature was kept at 22–27°C. The seedlings were watered (tap water with Multisol N fertilizer, Culterra, Muldersdrift, South Africa) regularly. The study was designed to monitor the plant response to fungal infection over time for 1, 3, 5, 7, and 9 days post infection (d.p.i.) as based on initial optimisation studies. The seedlings were planted in replicas of at least 10 plants per time point, and all plants were grown at the same time under the same environmental conditions. The complete experimental design included three biological repeats.

A pathogenic isolate of C. sublineolum (PPRI 7183) from fodder, grown, and maintained on potato dextrose agar (PDA), was obtained from the National Collection of Fungi, Plant Protection Institute, Agricultural Research Council (ARC), Pretoria, South Africa. The working sub-cultures were maintained on half-strength PDA solid media in Petri dishes. For spore production, the fungus was sub-cultured into 20% aqueous V8 medium (pH 3.9), which was prepared by mixing 100% V8 vegetable juice (Campbell Soup Company Camden, NJ, USA) and distilled water (1:5, v/v), and autoclaved. One hundred milli liter of the autoclaved media in 250 mL Erlenmeyer flasks were then inoculated with fungal mycelia (5 mycelia plugs per flask) from C. sublineolum cultures growing on half-strength PDA plates. The inoculated aqueous V8 media was incubated with constant shaking at 130 rpm, 12 h light cycle, for 7 d. The flask cultures were harvested after 7 d of growth by filtering the medium under vacuum through muslin cloth to remove the mycelial clumps. The spores present in the filtrate were pelleted by centrifugation at 5,000 × g for 15 min, washed by suspension (using autoclaved distilled water) and centrifugation, and diluted to the required concentration. The spore concentration was determined using a haemocytometer and light microscope at 400 × magnification, and adjusted to 10⁶ spores mL⁻¹.

At the four-leaf growth stage (25 days after sowing), the sorghum plant leaves were treated by spraying with the fungal spore suspension, adjusted to 10⁶ spores mL⁻¹, until run-off. The control plants were not sprayed. After inoculum application, the treated plants were incubated at 30°C in an incubator to provide 100% relative humidity, in darkness for 24 h. Following the 24 h incubation period, the plants were then exposed again to the same initial conditions: with cycles of 12 h fluorescence light (≈85 μmol m⁻² s⁻²) and 12 h darkness, and the temperature kept at 22–27°C. Post-treatment harvesting of the plants was done for all cultivars at 1, 3, 5, 7, and 9 d.p.i. by cutting off the leaves and immediate storage at −80°C until metabolite extraction. Similarly, the non-treated (negative controls) plants were harvested at 1, 5, and 9 days.

Metabolites were extracted from treated and non-treated plant leaves using 80% cold aqueous-methanol, in a ratio of 1:15 (w/v), at 4°C. The mixture was homogenized using an Ultra Turrax homogeniser, followed by sonication using a probe sonicator (Bandelin Sonopuls, Berlin, Germany) set at 55% power for 15 s, repeated 3 times. The homogenates were centrifuged at 5,000 × g for 10 min at 4°C. The supernatants were concentrated by evaporating to complete dryness and re-suspending the dried extracts in 300 μL 50% aqueous-methanol. The samples were then filtered through 0.22 μm nylon syringe filters into HPLC glass vials fitted with 500 μL inserts. The filtered extracts were kept at −20°C until analyzed. The methanol used was LC-grade (Romil Pure Chemistry, Cambridge, UK) and ultrapure water. The quality control (QC) samples were pooled samples prepared by pipetting and mixing aliquots of equal volume from all samples.

Ultra-high performance liquid chromatography coupled to high-definition mass spectrometry (UHPLC-MS) was performed on a Waters Acquity UHPLC coupled in tandem to a Waters photodiode array (PDA) detector and SYNAPT G1 Q-TOF mass spectrometer (Waters Corporation Milford, USA). Chromatographic separation of the aqueous-methanol extracts was done using a Waters HSS T3 C18 column (150 mm × 2.1 mm × 1.8 μm) thermostatted at 60°C. Although the T3 column is classified as a C18 reverse phase type, it is able to separate some polar compounds in addition to the non-polar compounds. Elution gradient was carried out with a binary solvent system consisting of 0.1% aqueous formic acid (solvent A) and 0.1% formic acid in acetonitrile (Romil Pure Chemistry, Cambridge, UK; solvent B) at a flow rate of 0.4 mL min⁻¹. The initial conditions were 98% A and 2% B and held for 1 min. A gradient was applied to change the chromatographic conditions to 30% A and 70% B at 14 min; and changed to 5% A and 95% B at 15 min. These conditions were held for 2 min and then changed to the initial conditions at 18 min. The analytical column was allowed to calibrate for 2 min before the next injection. The total chromatographic run time was 20 min and the injection volume was 2 μL. Each sample was analyzed in triplicate to account for any analytical variability. Solvent blanks and the QC samples were also analyzed in parallel with the sample extracts (described below).

High definition mass spectrometry analyses were performed on a Waters SYNAPT G1 Q-TOF MS system in V-optics operated in both positive and negative electrospray ionization (ESI) modes. Leucine encephalin (50 pg mL⁻¹), [M + H]⁺ = 55.2766 and [M–H]⁻ = 554.2615, was used as a reference calibrant, being continuously sampled every 15 s, producing an average intensity of 350 counts scan⁻¹ in centroid mode. Using this reference, the MassLynx™ software automatically correct the centroid mass values in the sample for small deviations from the exact mass measurement, giving typical mass accuracies between 1 and 3 mDa. The capillary and sampling cone voltages were 2.5 kV and 30 V, respectively. The extraction cone was set at 4.0 V. The source temperature used was 120°C and the desolvation temperature 450°C; cone gas flow 50 L h⁻¹ and desolvation gas flow of 550 L h⁻¹. A scan time of 0.2 s was used with a 100–1,000 Da mass range. The nebulisation gas used was nitrogen with a 700 L h⁻¹ flow rate. The data were acquired with different collision energies (MS^E) 0–30 eV to obtain as much structural information as possible of the detected compounds. The software used to control the hyphenated system and perform all data manipulation was MassLynx™ 4.1 (SCN 704, Waters Corporation Milford, USA).

The QC (pooled) samples were used to condition the LC-MS analytical system so as to assess the reliability and reproducibility of the analysis, and for non-linear signal correction (Godzien et al., 2015; Broadhurst et al., 2018). Sample acquisition was randomized and the QC sample (6 injections) was analyzed every 10 injections to monitor and correct changes in the instrument response. Furthermore, 6 QC runs were performed at the beginning and end of the batch to ensure system equilibration. Such sample randomization provides stochastic stratification in sample acquisition so as to minimize measurement bias. In the principal component analysis (PCA) space, the QC samples were clustered closely to each other, confirming thus the stability of the LC-MS system used, the reliability and reproducibility of the analysis. The blank samples (50% aqueous methanol) were randomly run to monitor background noise.

The centroid-raw data obtained from UHPLC-HDMS (both ESI positive and—negative modes) were pre-processed (peak picking, noise filtering, retention time (Rt) alignment, peak integration, and normalization) using MassLynx XS™ 4.1 software (Waters Corporation, Manchester, UK). The MarkerLynx™ application manager of the MassLynx software was used for matrix creation, producing a matrix of (Rt-m/z) variable pairs, with m/z peak intensity for each sample. MarkerLynx software parameters were set to process the 1–15 min Rt range of the chromatograms and m/z domain of mass range 100–1,000 Da. The Rts were allowed to differ by ± 0.2 min and the m/z values by ± 0.05 Da. The mass tolerance used was 0.01 Da, and the intensity threshold was 100 counts. Only data matrices that had noise level <50% (MarkerLynx cut off) were retained for downstream chemometric and statistical analyses. After data pre-processing, the number of metabolite features (Rt, m/z) in the clean data sets were 1536 in ESI positive and 2759 in ESI negative data sets. These data matrices were exported into SIMCA (soft independent modeling of class analogy) software, version 14 (Umetrics, Umeå, Sweden) for statistical analyses. Two unsupervised methods, PCA and hierarchical clustering analysis (HCA), and a supervised modeling, orthogonal projection to latent structures-discriminant analysis (OPLS-DA), were employed. This included the OPLS-DA S plots and the Variable Importance in Projection (VIP) plots. These multivariate methods attempt to highlight trends and groupings within a data set, subsequently facilitating the understanding of the relationships between- and within the samples (Trygg et al., 2007; Tugizimana et al., 2013).

To ensure that the biological question of the study is accurately answered as competently as possible, data scrutiny was meticulously done following the data pre-processing steps. This included assessment of the number of extracted features (<10,000 features, as a rule of thumb), applying the 80% rule (i.e., features found in <20% of the analyzed samples were removed) and monitoring the quality of data and stability of the analysis using QC samples. Data transformation methods such as centering, scaling or transformation were exploratively employed to put all variables on equal footing, minimize variable redundancy and adjust for measurement errors (Tugizimana et al., 2016). A non-linear iterative partial least squares algorithm (in-built within SIMCA; Nelson et al., 1996) was used to manage the missing values, with a correction factor of 3.0 and a default threshold of 50%. A seven-fold cross-validation (CV) method (Bro et al., 2008) was applied as a tuning procedure in computing the chemometric models; and only the components positively contributing to increase the prediction ability of the model (R1 significant components) were considered. Furthermore, thorough model validations were rigorously and consistently applied; and only statistically satisfactory models were examined and used in data mining for knowledge discovery. For PCA, the cumulative R² (explained variation) and Q² (predicted variation) of used models were higher than 0.5. Furthermore, for OPLS-DA models, the analysis of variance testing of cross-validated predictive residuals (CV-ANOVA), p-values were below 0.05. The specific values for these validation parameters (and others) are provided in the results section.

The study design information, LC-MS raw data, analyses and data processing information, and the metadata have been deposited to the EMBL-EBI metabolomics repository—MetaboLights database with the identifier (accession number) MTBLS791.

Total RNA was extracted from harvested leaf tissues, corresponding to the different time-points (1–9 d.p.i.) of each biological repeat, using the Trizol-reagent method (Invitrogen, Carlsbad, CA, USA). The extracted RNA samples were subjected to DNase treatment using DNase I (Thermo Scientific, Waltham, MA, USA). Concentrations were determined using a NanoDrop® ND-1000™Spectrophotometer (NanoDrop Inc., Wilmington, DE, USA). The RNA integrity of all samples were examined by electrophoresis on a 1.5% agarose gel in 1X Tris-Borate-EDTA (TBE) buffer and containing 0.5 μg mL⁻¹ ethidium bromide before use. The gels were visualized under UV light using a Bio-Rad Image Analyser and Quantity One™ Version 4.6.1 Software (Bio-Rad Laboratories, Johannesburg, South Africa). The total RNA samples were aliquoted and stored at −80°C for later use.

Real time PCR (qPCR) was used for sorghum gene expression analysis. Prior to quantification of the expression levels, the DNase-treated RNA were reverse transcribed to cDNA using a RevertAid™ Premium First Strand cDNA synthesis kit (Fermentas, Thermo Scientific, Waltham, MA, USA). The selected genes included: chitinase (PR3), pathogenesis-related protein 10 (PR10), flavonoid 3′-hydroxylase (F3'H), phenylalanine ammonia-lyase (PAL) and polyphenol oxidase (PPO). The gene-specific primer pairs (Supplementary Table S1) were designed using the “Primer Quest” tool (Integrated DNA Technologies, Coralville, IA, USA) from sequences obtained in on-line data bases (GenBank NCBI, www.ncbi.nlm.nih.gov/genbank). qPCR was performed to analyse the expression of each gene using a RotorGene-3000A instrument (Qiagen, Venlo, Netherlands) using the FastStart essential DNA Green Master Kit (Roche, Mannheim, Germany) according to the manufacturer's instructions. Ten micro liter of SYBR (FastStart essential DNA Green Master), 1 μL forward primer (1 μM final concentration), 1 μL reverse primer (1 μM final concentration), and 6 μL of DNase-free water were added to 2 μL of cDNA for amplification in a total volume of 20 μL. The cycling conditions were as follows: initial denaturation for 10 min at 95°C followed by amplification and quantification cycle repeated 40 times each consisting of 5 s denaturing at 95°C, 10 s annealing at primer specific temperatures, 20 s extension at 72°C. Two independent cDNA preparations were used with three technical replicates of each. Quantification of the relative changes in gene expression was performed using the relative standard curve method (Liu and Saint, 2002) with elongation factor 1-alpha (Elf α) and ubiquitin conjugating enzyme 18 (UBC18) as references genes. Data sets were statistically compared with the statistical analysis software GraphPad InStat v3 (GraphPad software, San Diego, CA, USA) using one-way analysis of variation (ANOVA) with Dunnet's post-test comparison of all treated samples vs. non-treated samples (control) at each time point. The confidence level of all analyses was set at 95%, and values with p < 0.05 were considered significant.

The symptomatic observations regarding the development of symptoms and disease severity of the MHL, BTT and SWT cultivars (Supplementary Figures S3, S4 and Supplementary Table S2) point to cultivar-related differential interactions between the sorghum plants and the hemibiotrophic C. sublineolum pathogen.

Hydromethanolic extracts of C. sublineolum-infected and non-infected sorghum plants were analyzed on a reversed phase liquid chromatography (LC) column coupled to a high-resolution quadrupole time-of-flight (QTOF) mass spectrometry (MS) detector system with electrospray ionization (ESI). This LC-ESI-QTOF-MS platform was combined with an untargeted approach to gather information on as many statistically significant metabolites as possible. Considering the inherent chemo-diversity, heterogeneity, and multi-dimensionality of extracted metabolomes, chromatographic separation is an essential step in untargeted metabolomics workflow, providing resolution of sample constituents (Tugizimana et al., 2013). Interfaced in-line with ESI-MS, the resultant LC-MS analytical platform allowed the simultaneous detection of multiple analytes with high sensitivity, providing deeper and more detailed insight into the metabolic composition of a biological sample. Distinct MS chromatograms indicated differential metabolic profiles of the analyzed samples. Figure 1 and Supplementary Figures S5A,B shows typical base peak intensity (BPI) mass chromatograms with differential peak population (presence and intensities), reflecting differences between samples from infected and non-infected plants, as well as cultivar-related differences.

To further elucidate the functional readouts of cellular physiological state(s) related to sorghum responses to C. sublineolum infection, chemometric analyses were applied to the collected LC-MS data. Following data processing (Boccard and Rudaz, 2014; Tugizimana et al., 2016), the created data matrices, with the number of defined features (Rt, m/z) being 1536 in ESI positive and 2759 in ESI negative data sets, were then exported into SIMCA (version 14) software for multivariate data analyses that included PCA, HCA, and OPLS-DA modeling. For the descriptive exploration of the overall structure of the pre-processed multi-dimensional data, unsupervised learning methods—PCA, and HCA—were used. These multivariate methods attempt to highlight descriptively trends and groupings within a data set, subsequently facilitating the understanding of the relationships between- and within the samples (Trygg et al., 2007; Tugizimana et al., 2013). PCA modeling, through the first two principal components (PCs), revealed treatment-related and cultivar-related sample clustering (Figures 2A,B and Supplementary Figure S6). Furthermore, in the PCA space, the QC samples are clustered closely to each other (and more or less in the middle of the 2D-plots), reflecting the stability of the LC-MS system used, and the reliability and reproducibility of the analysis (Godzien et al., 2015; Broadhurst et al., 2018). These sample groupings highlighted by the computed PCA models point to differential metabolic changes in sorghum plants responding to C. sublineolum. The samples from the MHL cultivar formed a clearly different group from the other two cultivars. This observation correlates to the indications from symptomatology (Supplementary Figures S3, S4 and Supplementary Table S2).

The PCA-extracted trends in the data were further examined by applying hierarchical clustering analyses on low-dimensional data generated from the PC analyses. Agglomerative HCA models were computed using Ward's linkage method (incremental sum of squares method) that considers between- and within-cluster distances when forming clusters, and the tree was sorted based on size (Szekely and Rizzo, 2005; Ji and Liu, 2010). The generated hierarchy of clusters was represented graphically on a dendrogram to evaluate whether some natural grouping emerges from the data—i.e., if the “metabolite space” actually contains several distinct subspaces. The computed HC models depicted two major distinct clusters corresponding to the samples from the very susceptible MHL cultivar grouping differentially and separate from the other two cultivars (BTT and SWT). Treatment-related (infected vs. non-infected) and time-related sub-clusters were also formed within each major cluster (Figures 2C,D and Supplementary Figure S6). Thus, both PCA and HCA modeling aided to evaluate descriptively the overall structure of the data, revealing underlying patterns and inner structures and sub-structures within the data: cultivar-related clustering, treatment-dependent groupings (infected vs. non-infected), and time-related variation (Figures 2C,D and Supplementary Figure S6). These observations evidently point to a biological phenomenon in the (extracted) metabolite space—differential metabolite profiles defining temporal cellular events related to the sorghum plants' responses to C. sublineolum infection.

For better biochemical interpretability and detailed assessment of the metabolic changes revealed by PCA and HCA in sorghum responding to the fungal infection, the supervised modeling method, OPLS-DA, was used. Evaluation of this multivariate (binary) classifier helps in extracting the metabolite variables underlying the discrimination between classes or groups (Trygg et al., 2007; Tugizimana et al., 2013). OPLS-DA is an extension to the supervised PLS-DA regression method, featuring an integrated orthogonal signal correction (OSC)-filtering method and, as such, OPLS-DA modeling has added interpretational and discriminatory benefits compared to PLS-DA (Bylesjö et al., 2006). The computed OPLS-DA models (Figure 3 and Supplementary Figure S7) to separate multivariate relationships into predictive (related to C. sublineolum infection) and orthogonal (unrelated to the treatment) variation, were validated with multivariate statistical tools and scrutinized by assessing the robustness, predictive ability, reliability and significance of the models. Some of the multivariate statistical tool used to validate calculated OPLS-DA models included R²–and Q² metrics, the analysis of variance testing of cross-validated predictive residuals (CV-ANOVA, p-value <0.05 as a cut-off), the receiver operator characteristic (ROC) curves, response permutation tests (with n = 50), and predictive testing (Eriksson et al., 2008; Tugizimana et al., 2016). The computed and validated OPLS-DA models (p < 0.05) used in this study were perfect classifiers and statistically reliable, with very good predictive capability: no signs of possible overfitting, as indicated by cross-validation; and none of the permutated models performed better than the original models in separating classes (Figure 3 and Supplementary Figure S7). These binary classifier models allowed to assess explicatively the treatment-related groupings (described by the unsupervised PCA and HCA models above) by extracting features (variables) responsible for differentiating sample groups (e.g., infected vs. non-infected).

Thus, the selection of discriminating features (signatory biomarkers characterized by unique Rt and m/z values) was carried out by evaluating the OPLS-DA loading S-plots (Figure 3B). To avoid overinterpretation of the models and variable selection bias, only features that were statistically significant in contributing to class separation were retained. Therefore, variables that combined both high covariation and correlation (as examined on S-plots) were considered to be statistically relevant as potential discriminant features (Wiklund et al., 2008; Tugizimana et al., 2013). However, since the S-plot is susceptible to data matrix changes due to correlation sensitivity and dependency on data structure, the statistical significance and discriminability of the potential markers derived from the S-plots were further investigated using different tests and tools such as the VIP plots, jackknife confidence intervals (used to estimate standard errors in a non-parametric way as an estimate of bias), variable trends, dot plots and descriptive statistics (Figures 3C,D). The VIP plots display VIP values as a column plot with jackknife uncertainty bars, providing a metric to assess the importance of the variables both to explain X and to correlate to Y, with the jackknife confidence intervals reflecting the variable stability (Galindo-Prieto et al., 2015; Tugizimana et al., 2016). Only S-plot-derived variables with VIP scores exceeding 1.0, with no (or minimal) overlap between groups (as indicated by dot—and trends plots), with positive jackknife confidence intervals and p-value < 0.05 (ANOVA, T-test), and demonstrating stable signals in the QC samples, were selected and retained as statistically significant and chemometrically contributing correctly to class separation. In a logical extension, these selected discriminant variables are regarded as essential chemical repertoires explaining the metabolic changes in sorghum, revealed by PCA—and HCA models (Figures 2C,D and Supplementary Figure S6). Accordingly, such features are fundamental elements for the biochemical interpretation of the chemometrically extracted information. These extracted features (markers) were then annotated (to the Metabolomics Standards Initiative, MI-level 2 annotation), as described in the experimental section and are reported in Table 1.

To identify the most significant metabolic pathways defining the sorghum defense responses, the MetPA (Metabolomics Pathway Analysis)—an integral module of the MetaboAnalyst bioinformatics tool suite (version 3.0; http://www.metaboanalyst.ca/)—was used. MetPA is a pathway analysis and visualization tool that combines several advanced pathway enrichment analysis methods along with the analysis of pathway topological characteristics to facilitate the elucidation of most relevant and altered pathways involved in the conditions under study (Xia et al., 2015; Chong et al., 2018. A representation of all MetPA-computed metabolic pathways displayed according to their significance or pathway impact in shown in Figure 4 and Table 2. The nine most significant pathways were (with some overlap): phenylalanine metabolism, stilbenoid and gingerol biosynthesis, flavonoid biosynthesis, flavone and flavanol biosynthesis, tryptophan metabolism, phenylpropanoid biosynthesis, aromatic amino acid biosynthesis, riboflavin-, and tyrosine metabolism. Furthermore, the topological characteristics of the phenylpropanoid—and flavonoid pathways are shown in Figure 5, illustrating that the two pathways are structurally highly interconnected, with some overlap as also shown in the topological graph-pathways generated from MetPA.

To complement these results, a biochemical and empirical network displaying metabolic relationship patterns between metabolites (indicated by OPLS-DA as signatory biomarkers) were performed. Figure 6 illustrates how the metabolites are connected based on biochemical relationships or structural similarity. The graphic representation and computed network parameters (e.g., clustering coefficient of 0.695; network density of 0.573) revealed a high interconnectivity of the OPLS-DA selected metabolites.

To enrich the metabolomic results with transcriptome insights, a phytoalexin-related gene (F3′H) and some defense-related genes (PAL and PPO, and PR-proteins PR3 and PR10) were selected and analyzed for expression levels in response to C. sublineolum infection. The general observation from these results points to time- and cultivar-related expression profiles; and all genes showed significant expression levels at different time intervals after pathogen inoculation (Figures 7, 8).

The phenotypic observations can be interpreted as BTT exhibiting a stronger defense response than SWT, with hypersensitive response-like lesions and purple color formation around the infected tissue. In contrast, the MHL cultivar exhibited little or no resistance against C. sublineolum. As mentioned, anthracnose development and severity in sorghum vary depending on the interaction between C. sublineolum (variation in virulence within the pathogen population), the genetics-based potential of the host plant to ward off infection, as well as environmental conditions (Tesso et al., 2012). Considering that in this study the sorghum cultivars were infected with a genetically uniform C. sublineolum isolate, in a controlled environment, the symptomatic differences thus reflect cultivar-related responses to the fungal infection. Although the plant defense responses triggered upon fungal infection are broadly similar across cultivars, the kinetics of these biochemical and cellular events and the relative abundance and timing may vary among cultivars (Liu et al., 2010; Tesso et al., 2012). Furthermore, such observation points to the inherent complexity of multi-layered plant innate immunity, of which omics-based studies have barely scratched the surface. Hence, characterizing metabolic phenotypes related to the sorghum defense responses to C. sublineolum infection would provide more insights into cellular pathways linked to underlying biochemical and molecular mechanisms operative in this specific plant-pathogen interaction.

For biochemical interpretation of the post-infection metabolic reprogramming in sorghum plants, as infographically described by the chemometric models, the statistically selected biomarkers/metabolites (Table 1) were further explained in the global metabolic interrelationships. Thus, metabolic pathway mapping and correlation network analyses were performed to elucidate the most relevant pathways and global dynamic metabolic networks involved in sorghum's responses to C. sublineolum infection. Both approaches exploited the relational properties present in the generated metabolomic data. Metabolic pathway analysis (or mapping) uses prior biological knowledge to map and analyse metabolites in an integrative manner, inferring significant pathways related to the study (Barupal et al., 2018; Rosato et al., 2018). On the other hand, metabolic network analysis methodology uses the high degree of correlation (biochemical and/or structural) existing in the generated metabolomic data to construct networks that characterize the complex relationship in measured metabolites (Toubiana et al., 2013; Grapov et al., 2015; Rosato et al., 2018).

Based on the chemometrically extracted metabolites (Table 1), pathway analysis with MetPA revealed that nine significant metabolic pathways out of a total of 24 pathways (impact score >0.10) were uniquely altered during the sorghum: C. sublineolum interactions. These most significant pathways include phenylalanine metabolism, flavonoid biosynthesis, phenylpropanoid biosynthesis, tryptophan metabolism, and riboflavin metabolism, among others (Figure 4 and Table 2). These results suggest that sorghum responses to C. sublineolum involves highly complex cellular reprogramming characterized by altered metabolism spanning a number of metabolic pathways; particularly the 9 significant pathways indicated in Table 2 and Figure 4, with phenylpropanoid and flavonoid biosynthesis pathways showing the highest hits. The constituents of these two metabolic pathways are the most widely occurring secondary metabolites found in the plant kingdom, exhibiting a broad range of biological functions including development, protection against abiotic and biotic stresses, modulation of essential physiological, and biochemical processes such as signal transduction, and transcriptional regulation (Cheynier et al., 2013; Petrussa et al., 2013).

Structurally, metabolites of phenylpropanoid and flavonoid pathways are phenolic compounds comprising an aromatic ring, with one or more hydroxyl groups, and include simple phenolic molecules to highly polymerised and conjugated compounds (Quideau et al., 2011; Petrussa et al., 2013). The two pathways are thus structurally highly interconnected, with some overlap as also shown in the topological graph-pathways generated from MetPA (Figure 5). Both pathways start with the conversion of phenylalanine to p-coumaroyl-CoA by phenylalanine ammonia-lyase (PAL), cinnamate-4-hydroxylase (C4H) and 4-coumaroyl:CoA-ligase (4CL) (Petrussa et al., 2013). The differentiation of the pathways rises with the formation of basic ring systems, with flavonoids showing a common three ring chemical structure (C₆-C₃-C₆). The wide diversity of both phenylpropanoids and flavonoids is then brought about by efficient decoration, modification, and amplification of core structures by a set of enzymes that are spatially and temporally coordinated and highly regulated (Quideau et al., 2011; Petrussa et al., 2013). The widespread presence of the phenolic compounds (both from phenylpropanoid and flavonoid pathways) at cellular, tissue and organ level is a clear indication of the multiple biological—and biochemical functions in plants Quideau et al., 2011; Petrussa et al., 2013) and sorghum is found to contain an array of these phenolic compounds (Althwab et al., 2015; Kang et al., 2016). Studies have shown that these phenolics play crucial roles in plant-fungal interactions as protectants: either as pre-formed antifungal compounds (phytoanticipins) or induced antifungal molecules (phytoalexins; Lattanzio et al., 2006). As presented in Figure 5, the topological characteristics of phenylpropanoid and flavonoid pathways display some of the phenolic compounds (in red) that exhibited dynamic changes, also providing some insights into the relational properties.

Thus, differential metabolite changes (Table 1 and Supplementary Figures S8, 9) observed in this study, evidently indicate that following the perception of C. sublineolum invasion, sorghum launched a complex arsenal of chemical defenses. These involved changes in phytohormone levels, indole compounds and modulation and activation of the pre-existing antifungal (phenolic) metabolites, as well as de novo biosynthesis and translocation of (flavonoid) phytoalexins. The metabolic reprogramming was cultivar dependent and was typically exemplified by the identification of quantitative changes in jasmonic acid (JA)- and salicylic acid (SA) conjugates, abscisic acid (ABA) and in the constitutive metabolites such as naringin, quercetin, and its conjugates (e.g., quercetin 3-O-rhamnoside), kaempferol di-glucoside and coumarate and sinapoyl conjugates, among others (Figure 5, Table 1, Supplementary Figures S8, 9). The phytohormones coordinate multiple physiological and biochemical processed in plants, such as growth, development and responses to abiotic and biotic stresses. The intricate web of crosstalk between various plant hormones, either through synergistic or antagonistic interactions, fine-tunes the regulation of plant immune responses, and is linked to biotrophic- or necrotrophic pathogen lifestyles (Shigenaga and Argueso, 2016). The occurrence of and cultivar-related changes in JA, SA, and ABA in response to infection by hemibiotrophic C. sublineolum, point to a multicomponent sorghum defense response (Table 1, Supplementary Figure S8).

In previous studies of different phytopathosystems, some of the phenolic metabolites, found to differentially accumulated in this study, have been reported to be involved in defense mechanisms: naringin has been documented to show antifungal activity, acting as a defense barrier. Moreover, quercetin -, naringenin -, and kaempferol conjugates have been shown to exhibit biological activity against fungal pathogens and a significant inhibition of spore germination (Lattanzio et al., 2006; Cheynier et al., 2013). If the levels and types of pre-existing antifungal phenolics are not enough to effectively limit the infection process, plant cells would respond by altering the pool of these phenolics, by increasing the levels or structurally modifying these phenolics to biochemically activate the compounds. This alteration in phenolics metabolism provides adequate substrates to polyphenol oxidase-catalyzed reactions that produce an environment unfavorable to fungal pathogenicity (Dicko et al., 2005; Lattanzio et al., 2006; Constabel and Barbehenn, 2008). The responses of sorghum to C. sublineolum infection by altering the phenolic matrix (as shown in this study—Figure 5, Tables 1, 2, Supplementary Figure S9) thus demonstrates the onset of active defense mechanisms by sorghum to stop the fungal invasion. This correlates to previous studies that observed deposition of phenolics in sorghum leaves and stalks of both resistant and susceptible varieties post C. sublineolum infection (Dicko et al., 2005; Anjum et al., 2013). Furthermore, the quantitative assessment of these phenolics indicates cultivar-dependent responses: overall, the NS5511 (bitter, BTT) cultivar had higher levels of phenolics compared to other cultivars, whereas the Mhlope (sweet white, MHL) cultivar showed low levels of these defense-related metabolites (Supplementary Figure S9). This correlates to the symptomatology observations: the MHL cultivar appeared to be severely affected by the C. sublineolum infection (Supplementary Figure S4).

In addition to the other observed changes in the flavonoid pathways, the de novo biosynthesis of 3-deoxyanthocyanidin phytoalexins, apigeninidin, and luteolinidin, together with the related flavonoids such as apigenin, apigenin conjugates, luteolin and luteolin di-glucoside (Figures 5, 6; Table 1, Supplementary Figure S9) were observed. As infographically captured by the unsupervised color-coded-PCA scores plot (Figure 7A), for instance, the presence and abundance (expressed as integrated peak area in the data matrix X) of the apigeninidin molecular feature started appearing from 3 d.p.i. with a clear increase over time. This evidently showed that in non-infected sorghum plants (at 0–9 d) there was no detectable presence of apigeninidin, and only infected plants were seen to have synthesized this phytoalexin, with increasing levels over time. The accumulation of these antifungal phytoalexins—apigeninidin and luteolinidin—was further confirmed by gene expression analysis of flavonoid 3′-hydroxylase (F3′H), which showed a time-related increasing expression profile (Figure 7B). The F3′H gene encodes one of the key enzymes responsible for the biosynthesis of these 3-deoxyanthocyanidin phytoalexins, unique antifungal compounds synthesized by sorghum (and related plants such as sugar cane) after fungal infection. The F3'H enzyme is responsible for the multi-step biochemical formation of 3-hydroxyantocyanidins from naringenin (Boddu et al., 2004; Liu et al., 2010; Poloni and Schirawski, 2014). Flavonoid hydroxylases are microsomal cytochrome P450 enzymes responsible for hydroxylation patterns of flavonoids. The expression of F3'H gene is induced by fungal infection and responsible for the biosynthesis of 3-hydroxylated flavonoids, as shown in Supplementary Figures S1, S2 (Boddu et al., 2004; Petrussa et al., 2013). Thus, the expression of the F3'H gene in sorghum responding to C. sublineolum supports and confirms the metabolic results, that is the presence and accumulation of these 3-deoxyanthocyanidin phytoalexins. Previous studies have also indicated that sorghum responses to fungal infection are characterized by the accumulation of an array of phenolic compounds, with major components being the phytoalexins—apigeninidin and luteolinidin—and their conjugates, arabinosyl-5-O-apigeninidin, 7-methylapigeninidin, and 5-methoxyluterolinidin (Liu et al., 2010; Poloni and Schirawski, 2014). These experiment-based observations correlate with our results regarding the biosynthesis of phytoalexins apigeninidin, luteolinidin, and other related flavonoids (Figure 7, Table 1, Supplementary Figure S9). Functionally, these unique class of flavonoid phytoalexins are initially synthesized in the cytoplasm of epidermal sorghum cells following fungal infection, and accumulate in inclusion bodies. These are translocated toward the site of fungal invasion, where they are then released in active form and kill both the fungus and cells that synthesized them (Poloni and Schirawski, 2014; Meyer et al., 2016). The kinetics and intensity of this response appear to be cultivar-dependent, as reflected by the results in Figure 7C (and Supplementary Figure S9), showing comparative cultivar-related differences in the fold change of apigeninidin and luteolinidin, respectively.

The observations from this study and supported by the literature, point clearly to a complexly coordinated and highly regulated dynamic (and cultivar-dependent) metabolic reprogramming in sorghum responding to C. sublineolum infection. The functionally altered sorghum metabolism involved a range of different metabolic pathways (Table 2), which exhibit a complex interconnection as illustrated by the phenylpropanoid and flavonoid pathways (Figure 5). Pathway enrichment and overrepresentation analyses with the MetPA tool thus facilitated biochemical interpretation by integrating biological domain knowledge (i.e., biochemical pathways) with experimental results, revealing underlying relevant metabolic pathways (Figure 4 and Table 2). However, these pathway-based methodologies rely on predefined pathways, and fail for instance to capture linkage information of metabolites belonging to multiple pathways. This may not accurately represent the complexity of biological systems, subsequently providing limited insights into underlying mechanisms in the crowded cellular milieu, and spatial and temporal regulation of organismal reprogramming (Kruger and Ratcliffe, 2012; Toubiana et al., 2013; Barupal et al., 2018).

Hence, to gain more insights into possible global biochemical and molecular frameworks that choreograph the response of sorghum to C. sublineolum infection, a network analysis approach was adopted. This methodology uses the high degree of correlation (biochemical and empirical, in this study) existing in the experimentally generated metabolomic results to construct networks that characterize the complex relationship in measured metabolites (Kruger and Ratcliffe, 2012; Toubiana et al., 2013). Unlike pathway analysis, correlation-based approaches build metabolite networks according to relational patterns observed in the experimental data, and help identify altered graph neighborhoods, which do not depend on any predefined biochemical pathways. Such mathematically constructed cartography allows the characterization of the molecular and cellular states induced by pathway interconnections under given experimental conditions. In the computed network, each metabolite is represented by a node, and in contrast to pathway analysis, the links between nodes correspond to the level of mathematical correlation between each pair of metabolites (Kruger and Ratcliffe, 2012; Toubiana et al., 2013; Rosato et al., 2018). A biochemical/chemical similarity network analysis was accordingly applied to calculate and display relationship patterns between precursor and product metabolite reactant pairs, and molecules sharing a high degree of structural similarity, with Tanimoto coefficient ≥0.7. The resultant graphic representation and computed network parameters (e.g., clustering coefficient of 0.695; network density of 0.573) revealed a high interconnectivity of the OPLS-DA selected metabolites (Figure 6). These were statistically shown to explain class separation (e.g., infected vs. non-infected), thus describing the metabolic changes in sorghum revealed by the PCA and HCA models (Figures 2 and Supplementary Figure S6). Phenylpropanoids and flavonoids appear to form a central hub in the topology of the network, also maintaining high interactions with other metabolite categories, such as indolic or amino acid related metabolites (Figure 6).

Furthermore, the network analysis showed that separate from the major cluster, there were other small clusters that were formed: (i) luteolinidin, indole, indole-3-acetalamide, and hesperidin, (ii) kaempferol-diglucoside and p-coumaroylquinate, and (iii) jasmonate and jasmonic acid methyl ester (Figure 6). This differential clustering may possibly suggest different regulation of these metabolites in the concerted metabolic reprogramming of sorghum defense responses to the fungal infection. The computed high interconnectivity of nodes in the network demonstrates highly correlated biochemical and structural metabolic relationships that coordinate the altered sorghum metabolism in sorghum: C. sublineolum interactions. The significant metabolic pathways underlying sorghum responses to the fungal infection as revealed by pathway analysis with MetPA (Figures 4, 5 and Table 2), are highly interconnected, as demonstrated by the network topology (Figure 6). This points to regulatory hubs in the biochemical network (different clusters in the network), because the correlation matrix of metabolite pairs is a fingerprint of the enzymatic and regulatory reaction networks (Kruger and Ratcliffe, 2012; Toubiana et al., 2013). Furthermore, the results from both pathway and network analyses emphasize the centrality of the phenylpropanoid and flavonoid metabolic pathways in the sorghum responses. Although more sophisticated approaches (e.g., Gaussian graphical models and Bayesian networks; Kim et al., 2011; Kayano et al., 2013) may be needed to decouple direct from indirect variable associations, thus helping the identification of conditionally independent pairwise metabolic relationships, the methodology used in this study comprehensively captured essential features of sorghum defensive metabolism against C. sublineolum infection.

The observed expression of these genes evidently indicates a multiphase defense state in the plants over time. Furthermore, these gene expression results clearly corroborate the information from metabolomics analyses and chemometric models that sorghum responses to C. sublineolum infection is time- and cultivar-related. The significant metabolic pathways underlying the response to the fungal infection were found to be the same in all the three cultivars investigated in this study (Figure 4 and Tables 1, 2). However, it is apparent that the kinetics, magnitude, and timing of the responses vary with cultivars: different levels of significant metabolites (expressed as fold change—Table 1) and differential levels and kinetics of gene expression (Figure 8).

Both PAL and PPO are essential enzymes in phenylpropanoid and flavonoid metabolism, leading to the biosynthesis of structural barrier components such as lignin, formation of antimicrobial phenolics such as phytoanticipins and de novo biosynthesis of phytoalexins in plant defense events (Mengiste, 2012; Anjum et al., 2013). In the case of the BTT (B) cultivar, the PAL gene expression exhibited a bi-phasic pattern, with an initial response at days 3–5, decreasing at day 7, followed by a strong increase at day 9. Similar bi-phasic responses have been reported in plants exhibiting a high level of resistance to pathogen attack (Ding et al., 2011) and in this case might be related to the biotrophic vs. necrotrophic stages of the infection. In the case of the SWT (S) cultivar no similar early PAL response was observed. For both cultivars, the highest expression of the PAL gene was observed from 7 d.p.i. onwards, corresponding to the necrotrophic stage of infection that is accompanied with cell destruction and eventual death. In the case of MHL (M), PAL expression levels were very low compared to other two cultivars (Figure 8), suggesting that this cultivar is unable to launch an effective defense response which correlates with the phenotypic observations as described earlier. A clear expression of PPO was also observed, with highest levels at 1 d.p.i. in the SWT cultivar, 5 d.p.i. in BTT, and 9 d.p.i. in MHL (Figure 8). Here, no obvious correlation with the type of phytopathogenic interaction could be deduced. However, the substantial expression of both PAL and PPO genes in sorghum plants responding to C. sublineolum clearly correlates and supports the metabolic alterations elucidated by metabolomic analyses, particularly changes and accumulation of phenolic compounds, and including de novo biosynthesis of the phytoalexins: apigeninidin and luteolinidin (Figures 4, 5, Tables 1, 2). The general relative low levels of PAL and PPO (up to 5 d.p.i) expression in the MHL cultivar (Figure 8) corresponds to the low level of observed defense-related phenolics (Table 1, Figure 7). Some previous studies have also reported significant expression levels of both PPO and PAL in sorghum plants responding to fungal infection, with substantial accumulation and channeling of phenolic compounds to combat the fungal proliferation and colonization (Basavaraju et al., 2009; Anjum et al., 2013).

The two pathogenesis-related genes (PR3 and PR10) were substantially expressed in all cultivars after the pathogen inoculation, with time- and cultivar-dependent expression profiles: early increased expression levels of PR3 in all cultivars followed by a decrease over the time-course of infection (expect in MHL-cultivar); a Gaussian-type expression profile for PR10 in both BTT (B) and SWT (S) cultivars, with BTT showing the highest expression levels of PR10 (at 3–5 d.p.i.) compared to other cultivars, which then decreased (Figure 8). Generally, plants are known to express chitinases (PR3) soon after infection for endolytical hydrolysis of microbial cell walls (Heil, 2002). Studies have shown that in plant: hemibiotrophic pathogen interactions, early defense mechanisms are characterized by a transient induction of chitinase (PR3) genes with the onset of the biotrophic interaction, which are suppressed with the progression of the infection into the necrotrophic phase (Münch et al., 2008; Vargas et al., 2012). This evidently correlates to the expression profiles of PR3 in this study—an early significant expression of the gene followed by a remarkable decrease at 5 d.p.i. (Figure 8), which points to an asymptomatic interaction phase. However, the MHL cultivar showed a different profile with a second increase of the PR3 expression, which could imply a multiphase response or simply other underlying cellular processes. On the other hand, PR10 is a member of a group of intracellular defense-related proteins with ribonuclease-like activity. The PR-10 group is a multigene family having cis regulatory elements responsive to various signals like ABA, SA and JA (Sudisha et al., 2012). In general, the PR-10 gene family shows non-specific induction patterns to pathogen and pathogen-derived molecules. These defense-related proteins are widely spread and conserved within the plant kingdom, and are induced following pathogen attack in a wide variety of plant species (Mcgee et al., 2001). Furthermore, the PR10s also exhibit similar amino acid sequence to food and pollen allergens, which point to diverse biological functions (Mcgee et al., 2001; Edreva, 2005). Activation and accumulation of PR10 gene transcripts have previously also been observed in sorghum infected by C. sublineolum, exhibiting increasing expression levels over time (Anjum et al., 2013). Such gradual increase of PR10 expression levels corresponds to the upregulation patterns observed in this study (Figure 8).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.