Overview of the transcriptomic and metabolomic dataset differences between the three wheat genotypes
To investigate the global transcriptomic profile of the three wheat genotypes, we grew the plants under a controlled environment for 11 d in which the three genotypes had the same phenology (Figure S1). Our recent investigation of the young wheat Svevo and Zavitan showed a massive chemical variation [29]. Therefore, we focused our transcriptomic assays on young seedlings and sampled them to the top of the second leaf. A comparison of transcriptomic data (Illumina RNAseq) was performed with annotated gene models found in the Chinese Spring reference genome sequence [19,35]. The mapping sequence reads, as compared to this reference, show a high similarity between the genotypes: Svevo, 95.85%; Zavitan, 94.68%; and Chinese Spring, 96.72% from the total mapped reads (Table S1). Due to differences in the number of subgenomes, we eliminated all the transcripts that were annotated to the D subgenome or an unidentified subgenome (U). This analysis revealed 42,474 transcripts (Table S2). The total transcript levels were used to conduct a principal component analysis (PCA) plot. As presented in Figure 1A, the PCA indicated that samples from each genotype were clustered with one another, while the genotypes were totally separated from each other, with component 1 (52.5%) showing a separation of Zavitan from the other cultivated wheats. To examine the metabolic diversity of the three wheat genotypes, we performed an untargeted metabolic analysis using liquid chromatography-time of flight-mass spectrometry (LC-QTof-MS). The PCA plot was performed using 878 mass features of the negative ion mode and showed massive metabolic differences between the three genotypes while component 1 (36.2%) showed a separation of Zavitan from the other cultivated wheats (Figure 1B). The mass feature dataset is presented in Table S3. Overall, both transcriptomic and metabolic analyses of the three genotypes showed a unique pattern for each one.
Clustering of expression patterns and gene ontology in the transcriptomic data
The DESeq2 output of the LRT process was subjected to an rlog transformation, and the resulting heatmap clearly divided the overall transcriptional profiles of the two cultivated genotypes (Svevo and Chinese Spring) from the wild genotype (Zavitan) (Figure S2). The analysis of differentially expressed genes (q-value < 0.05; |log2 (fold change)| >= 1) resulted in 8,735 unique transcripts. Following the LRT process, DESeq2 was performed on the results with the genotypes as the differentiating factor. We estimated the cluster number of the results using clusGap [41], which suggested dividing the data into eight clusters. The k-means analysis was performed on scaled and centered rlog values, and each cluster is represented by the Z-score (standard score) of gene expression from the set of genes showing similar response patterns of the three wheat genotypes. The eight clusters were divided into five expression groups, derived from trends observed in the standard scores: i) clusters 1 and 2: genes with a high level in Zavitan relative to Svevo and Chinese Spring. Cluster 2 presents genes with a higher Z-score than cluster 1. ii) Clusters 3 and 4: genes with a lower level in Zavitan than in Svevo and Chinese Spring. Cluster 4 presents genes with a lower Z-score than cluster 3. iii) Clusters 5 and 6: genes with moderately lower levels in Svevo than in Zavitan and Chinese Spring. Cluster 5 presents genes with a higher Z-score than cluster 6. iv) Cluster 7: genes that have lower levels in Zavitan than in Svevo and Chinese Spring. v) Cluster 8: genes that have lower levels in Chinese Spring than in Svevo and Zavitan (Figure 2). The distribution of genes into the eight clusters is presented in Table S4.
Pathway enrichment of the transcriptomic data
To elucidate the metabolic processes that are involved in each gene cluster shown in Figure 2, an over-representation pathway enrichment analysis was performed using MetGenMAP [42] to compare them to rice orthologues (LOC gene ID; Table S5). In Table 1, the significantly enriched pathways of each cluster are described. The pathways that were significantly enriched in clusters 1 and 2 (high in Zavitan) were mainly associated with polyamine biosynthesis and sugar degradation. The genes related to threonine and homoserine biosynthesis, phenylpropanoids biosynthesis, cell structure biosynthesis (cellulose) phospholipases, and ascorbate biosynthesis were associated with a low expression in Zavitan (clusters 3 and 4). The pathways that were significantly enriched in clusters 5 and 6 (slightly low in Svevo) were mainly associated with the isoprenoid phosphate pathway, glycine and glycerol degradation, brassinosteroid, jasmonic acid (13-lipoxygenase; 13-LOX) and gibberellin phytohormone biosynthesis, N-acetylgalactosamine biosynthesis, TCA-, Calvin-, and γ-glutamyl cycles and sugar degradation. The pathways that were significantly enriched in cluster 7 (slightly low in Zavitan) were mainly associated with glutathione, β-alanine, methionine, homocysteine, and cysteine biosynthesis, as well as gluconeogenesis, tyrosine degradation, phospholipid desaturation, and fatty acid β-oxidation. Additionally, pathways that were significantly enriched in cluster 8 (slightly low in Chinese Spring) were mainly associated with phytohormone cytokinin biosynthesis. Together, these observations indicate a unique gene expression for each wheat genotype that involves diverse pathways from primary and secondary metabolites, phytohormones, oxidation state, and cell wall.
To determine which gene ontology categories were represented, we conducted a Singular Enrichment Analysis (SEA) using agriGO v2 [43], with default parameters and IWGSC IDs as the background. In Figure S3, the GO terms of cluster 5 are presented. The significant processes are related to isoprenoid biosynthesis, the oxidation-reduction process, lipid metabolism, and photosynthesis. These results are partially redundant with the pathway enrichment results (Table 1) referring to the lipid/13-LOX biosynthesis, as well as isoprenoid biosynthesis.
Metabolic and physiological characterization of the three wheat genotypes
Several pathways related to central metabolism, such as carbohydrate, amino- and fatty acids and the TCA cycle, showed lower expression in Svevo than in Zavitan and Chinese Spring (Table 1). Therefore, we analyzed the metabolic profile of 11-day-old seedlings using gas chromatography-mass spectrometry (GC-MS) and quantified a total of 73 metabolites, including amino acids, organic acids, sugars, and sugar alcohols. The normalized data of the 24 significantly different metabolites (one-way-ANOVA) are presented in Table 2. The results showed that three aromatic amino acids, organic acids, and sugars were significantly lower in Svevo. Other metabolites that showed no significant differences are presented in Table S6.
The pathway enrichment analysis revealed genes related to polyamine biosynthesis, which serve as substrates for oxidation reactions that produce hydrogen peroxide (H2O2) both intra- and extracellularly [44]. One of the earliest signaling roles in many environmental stresses involves reactive oxygen species (ROS). As the most stable ROS, hydrogen peroxide plays a crucial role in physiological processes in plants such as growth, development, and reproductive growth, and it is also involved in defense against pathogens and diseases [45]. We measured hydrogen peroxide in the leaves of 11-day-old plants using DAB staining (3,3′-Diaminobenzidine). The leaves of Zavitan generated more dark brown precipitate than the other two genotypes. The sodium phosphate control treatment showed no precipitate (Figure 3A). This suggested that the oxidative status of Zavitan is higher than Svevo and Chinese Spring.
The effect of various abiotic and biotic stresses on the photosynthetic apparatus is inevitably associated with the formation of harmful ROS [46,47]. As an indication for the differences in the photosynthetic apparatus, we measured the total chlorophyll content by quantifying the levels of chlorophyll a and b [48]. As described in Figure 3B, Zavitan had significantly lower levels of total chlorophyll than Svevo and Chinese Spring leaves.
Genes related to sucrose and starch degradation were expressed at higher levels in Zavitan than in the other two genotypes (cluster 1). Previous studies linked the effect of carbohydrate metabolism to the condition of the photosynthetic apparatus [49] and water content [50]. Therefore, we also measured the water content in the leaves, which is a critical indication of plant response to different environmental stresses [51–54]. The analysis indicated that Zavitan leaves have significantly higher water content than Svevo and Chinese Spring (Figure 3C). Overall, these results support our transcriptomic analysis and show some similarity to the cultivated wheat genotypes relative to the wild emmer.
Differences in gene expression related to chemical and physical defense strategies
The pathway enrichment suggested that several pathways related to plant defense mechanisms differ between the genotypes. Therefore, we decided to further investigate the gene expression of the chemical and physical defense mechanisms. The major class of specialized metabolites in several of Gramineae grass family belongs to the benzoxazinoids, which are abundant in Zea mays (maize) and wheat species [55]. We searched the BREADWHEATCYC2.0 database (https://www.plantcyc.org/databases/breadwheatcyc/2.0) and the literature, and we generated a list of known and putative Bx genes (the benzoxazinoid biosynthetic genes named Bx1 through Bx14) [15,56,57]. In the case of unknown Bx genes in wheat (Bx6, Bx7, and Bx10-14), we aligned the protein sequences known in maize [58–61] to the wheat database Ensemble Plant. For the full list, see Table S7. The heatmap of the Bx genes indicated that the samples of cultivated wheat Svevo were closer to Chinese Spring than to Zavitan (Figure 4A). It also showed that some genes, such as cytochrome P450, Bx3 and Bx5, were expressed at higher levels in Zavitan, while other genes, such as the downstream glucosidase and O-methyltransferase were higher in Svevo and Chinese Spring. This suggested some similarity between Svevo and Chinese Spring regarding the of the benzoxazinoid biosynthesis genes.
Structures such as trichomes, present on the leaf surface, may affect aphid probing behavior, hinder movement, and serve as a source of toxic metabolites [62,63]. Therefore, we compared the gene expression of trichome formation-related genes as representative of a physical barrier. We searched the literature and found some evidence for trichome formation genes in rice including Glabrous Rice 1, encoding a homeodomain protein [64], and the pubescence gene GL6 [65]. We also identified maize protein homologs including i) HD-ZIP IV transcription factor OCL4, which is necessary for trichome patterning [66]; ii) UMC1273, which is a protein trichome birefringence-like 39; iii) UMC1601 protein trichome birefringence-like 28; iv) AY110056 protein trichome birefringence-like 26; and v) prm5 (powdery mildew resistant protein 5). For the full list, see Table S8. The heatmap did not show a clear pattern, as two Svevo samples were clustered with Zavitan and one was clustered with Chinese Spring (Figure 4B).
Quantification of changes in benzoxazinoid levels during plant development
We focused our analysis of defense metabolites on BXDs, which are major deterrent compounds and have been demonstrated to play a role in chemical defense in wheat leaves [67,68]. The most abundant BXDs in maize are 4-dihydroxy-7-methoxy-1,4-benzoxazine (DIMBOA), and 2-hydroxy-4,7-dimethoxy1,4-benzoxazine (HDMBOA) [67–69]. The end product of benzoxazinoid biosynthesis is a glucoside that reduces toxicity [68]. BXDs are abundant in wheat seedlings and show the highest activity in plants at the juvenile stage [13–15]. In the experiment, plants were grown for 11, 15, and 18 days after germination; their second leaf was harvested, and their levels of DIMBOA and DIM2BOA-Glc were measured using HPLC (Figure 5A-B). Both DIMBOA and DIM2BOA-Glc were detected only in the cultivated wheat genotypes Svevo and Chinese Spring, while they were below detection levels in the wild emmer Zavitan. The genotypes Svevo and Chinese Spring showed the highest levels of DIMBOA at day 11, which gradually declined until day 18. In Chinese Spring, the DIMBOA levels declined more rapidly than in Svevo. DIM2BOA-Glc levels showed different accumulation patterns than DIMBOA, as the highest levels were detected on day 15 when Chinse Spring was much higher than Svevo. We also identified five unknown BXDs according to their UV spectra (Figure S4) and quantified their peak area. The total unknown BXD peak areas are presented in Figure 5C. The results revealed that the unknown BXDs were more abundant in Chinese Spring, mainly in the youngest seedlings. In Svevo, only mild levels of unknown BXDs were accumulated, while in Zavitan, unknown 8.3 was detected only after 11 days. These data revealed a significant variation of BXDs, which vary between wheat genotypes and plant age.
Calculating the trichome density
To explore the physical defense strategy of wheat plants against aphid invasion, we evaluated the trichome density on the leaf surface and edge as a physical barrier (Figure 6). Among the three wheat genotypes, the highest number of trichomes on leaf edges and surface density was in Zavitan (F2,224 = 498.36, P value < 0.0001 and F2,372 = 268.32, P value < 0.0001 respectively). The number of trichomes on leaf edges was similar in Svevo and Chinese Spring, and Svevo possessed a higher trichome surface density than Chinese Spring (Figure 6A). The effect of time after germination on trichome density and on the amount on the edge was also significant (F2,224 = 5.30, P value <0.0056 and F2,372 = 16.00, P value <0.0001, respectively). In Zavitan and Svevo, the number of leaf edge trichomes slightly increased over time, while their surface density was reduced (Figure 6B). The length of trichomes observed on the leaf surfaces of Zavitan and Chinese Spring were longer than on Svevo (Figure 6C). Altogether, this suggested that Zavitan has more physical barriers, such as trichomes, than the other two genotypes.
Evaluate aphid reproduction on different wheat genotypes
To evaluate aphid performance on wheat seedlings, the Rhopalosiphum padi aphid was selected [32,70]. Ten adult aphids were applied to the second leaf of wheat 11- and 14-day-old seedlings for 96 h, and then the total number of aphids was counted on days 15 and 18, respectively (Figure 7). In general, the cultivated Chinese Spring genotype showed different rates of aphid reproduction (F2,59 = 23.40, P value <0.0001). Although the wild emmer wheat Zavitan has been shown to be more susceptible to aphids than Svevo, these differences were not significant. The same trend of aphid performance was repeated on infested seedlings at both days 15 and 18, while the 18-day-old seedlings were more susceptible to aphids than the 15-day-old ones (F1,59 = 168.03, P value <0.0001). Overall, the results indicated that the cultivated bread wheat was more resistant to R. padi at both infestation time points (Figure 7). Furthermore, we measured the BXD levels upon aphid infestation. The results indicated that no significant differences in BXDs were observed between the control and the R. padi-infested plants for 96 h (Table S9). This suggested that the BXD levels are constitutive and not induced by aphids [60].