Overview of the transcriptomic dataset differences between the three wheat genotypes
To investigate the global transcriptomic profile of the three wheat genotypes, leaves from 11-day-old seedlings were collected from genotypes with the same phenology (two-leaf stage; Figure S1). A comparison of transcriptomic data was performed with annotated gene models found in the Chinese Spring reference genome sequence [20,38], and the mapped sequence reads showed 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 1, the PCA plot indicated that samples from each genotype were clustered with one another, while the genotypes were totally separated from each other, with component 1 (45%) showing a separation of Chinese Spring from the tetraploid wheat genotypes, and component 2 (36%) separating the wild (Zavitan) and cultivated (Svevo and Chinese Spring) wheat genotypes. Overall, the transcriptomic analysis of the three genotypes showed a unique pattern for each one.
Clustering of expression patterns and pathway enrichment and gene ontology in the transcriptomic data
In order to detect differentially expressed genes in the RNA-seq data, the number of sequences associated with each gene (counts) for each sample was used to statistically compare between- and within-condition variability, by a negative binomial generalized linear model, using the R package DESeq2. The DESeq2 output, compared using the likelihood ratio test (LRT), 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), presented in Figure S2. The analysis of differentially expressed genes (q-value < 0.05; |log2 (fold change)| >= 1) resulted in 8,735 unique transcripts. We estimated the cluster number of the results using clusGap [44], 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 in the three wheat genotypes (Figure 2). 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 higher level in Zavitan than in 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 S3.
To elucidate the metabolic processes, an over-representation pathway enrichment analysis was performed on each cluster using MetGenMAP [45], using rice orthologues (LOC gene ID; Table S4). Table 1 describes the significantly enriched pathways of each cluster. 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, phenylpropanoid biosynthesis, cell structure biosynthesis (cellulose) phospholipases, and ascorbate biosynthesis were associated with 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 [46], with default parameters. The International Wheat Genome Sequencing Consortium (IWGSC) database gene IDs were used as background in the SEA. In Table S5, the GO terms of the biological processes of all eight clusters are presented in pairs. In Cluster 5 for example, the comparison between Zavitan-Svevo and Zavitan-Chinese Spring revealed biological processes related to isoprenoid biosynthesis, lipid metabolism, oxidation-reduction process, photosynthesis, and tetrapyrrole metabolism. The results of Cluster 5 are partially redundant with the pathway enrichment results (Table 1), regarding the lipid/13-LOX biosynthesis and the isoprenoid biosynthesis. The functional enrichment terms in this cluster are the most consistent with the differences in defense adaptations between the three genotypes. Other clusters included functions such as the organic acid metabolic process, the cellular amino acid metabolic process, cell redox homeostasis and metal ion transport (Cluster 7), the oxidation-reduction process (only differentially expressed genes of Svevo vs. Chinese Spring; Cluster 6), and the phenylalanine catabolic process (only differentially expressed genes of Svevo vs. Zavitan; Cluster 3). The remaining clusters included functional enrichment terms, such as protein phosphorylation, lipid transport, and recognition of pollen.
Characterization of metabolic and physiological differences between the three wheat genotypes
Several pathways related to central metabolism, such as carbohydrate, amino, and fatty acids and the TCA cycle, showed variation in gene expression between the three genotypes (Table 1). Therefore, we performed a metabolic analysis of 11-day-old wheat seedlings using gas chromatography-mass spectrometry (GC-MS). This analysis revealed the relative levels of 72 metabolites, including amino acids, organic acids, sugars, and sugar alcohols (Table S6). Table 2 shows the relative levels of 24 significantly different metabolites using a one-way-ANOVA. The results indicated that three aromatic amino acids, organic acids, and sugars were significantly lower in Svevo. Other metabolites showed no significant differences and are presented in Table S6.
Both the pathway enrichment analysis of Cluster 2 (Table 1) and the GO terms of Cluster 5 (Table S5) revealed that genes related to polyamine biosynthesis and the oxidation-reduction process are differently expressed between the three genotypes. Polyamines serve as substrates for oxidation reactions that produce hydrogen peroxide (H2O2) both intra- and extracellularly [47]. 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 [48]. 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 [49,50]. In addition, both oxidation-reduction processes and photosynthetic processes were significantly represented in the GO term analysis (Table S5, Cluster 5). As an indication of the differences in photosynthesis, we measured the total chlorophyll content by quantifying the levels of chlorophyll a and b [51]. As described in Figure 3B, Zavitan leaves 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 [52] and water content [53]. Therefore, we also measured the water content in the leaves, which is a critical indication of plant response to different environmental stresses [54–57]. As shown in Figure 3C, Zavitan leaves have significantly higher water content than Svevo and Chinese Spring leaves. The leaves’ fresh, turgid, and dry weights are presented in Table S7. Overall, these results support our transcriptomic analysis and show some similarities to the cultivated wheat genotypes relative to the wild emmer.
Differences in gene expression related to chemical and physical defense adaptations
The pathway enrichment suggested that several pathways related to plant defense mechanisms differ between the genotypes. These include the 13-LOX and 13-HPL pathway, which is related to JA biosynthesis [58], isopentenyl diphosphate biosynthesis (terpenoids) [59], both in Cluster 5, and flavonoid biosynthesis [60] in Cluster 3 (Table 1). Therefore, we further investigated the gene expression of the benzoxazinoids as chemical toxic compounds and trichomes as a physical defense mechanism. To generate a gene list, we searched the BREADWHEATCYC2.0 database (www.plantcyc.org) and the literature (the benzoxazinoid biosynthetic genes named Bx1 through Bx14) [16,61,62]. In the case of the unknown Bx genes in wheat (Bx6, Bx7, and Bx10-14), we aligned the maize protein sequences [63–65] to the wheat Ensemble Plant database (see Table S8 for the full gene list). The expression values of Bx genes in each genotype are presented in the heatmap, and the samples are sorted by hierarchical clustering (Figure 4A). The heatmap indicated that the samples of the cultivated wheat Svevo were closer to Chinese Spring than to Zavitan. 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 glucosidases and O-methyltransferases, were higher in Svevo and Chinese Spring. This suggested some similarities between Svevo and Chinese Spring regarding the benzoxazinoid biosynthesis genes.
We also compared the gene expression of trichome-formation-related genes as representative of a physical barrier [66,67]. To generate a gene list, we searched the literature and found some evidence for trichome formation genes in rice, including Glabrous Rice 1, encoding a homeodomain protein [68], and the pubescence gene GL6 [69]. We also identified maize protein homologs, including i) HD-ZIP IV transcription factor OCL4, which is necessary for trichome patterning [70]; 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; (see Table S9 for the full gene list). The expression values of the trichome formation genes in each genotype are presented in the heatmap, and the samples are sorted by hierarchical clustering (Figure 4B). The heatmap did not show a clear pattern, as two Svevo samples were clustered with Zavitan and one was clustered with Chinese Spring.
Quantification of changes in benzoxazinoid levels during plant development
We focused our analysis of defense metabolites on benzoxazinoids (BXDs), which are major deterrent compounds and have been demonstrated to play a role in chemical defense in wheat leaves [71,72]. The BXDs are abundant in wheat seedlings and show the highest activity in plants at the juvenile stage [14–16]. In this experiment, we measured, by HPLC-UV, the levels of BXDs in young wheat plants: 11, 15, and 18 days after germination. Overall, three BXD compounds were detected and identified, including 4-dihydroxy-7-methoxy-1,4-benzoxazine-3-one (DIMBOA) ,4-dihydroxy-7,8-dimethoxy-1,4-benzoxazin-3-one glucoside (DIM2BOA-Glc), and 2-hydroxy-4,7-dimethoxy-1,4-benzoxazin-3-one glucoside (HDMBOA-Glc) as presented in Figure 5A-C. These compounds were further annotated by comparing the pattern of their fragment using UPLC-qTOF-MS and previous studies (Table S10; [73,74]). All compounds were detected in at least one sampling time point, in the domesticated wheat genotypes, while they were below detection levels in the wild emmer Zavitan. The genotypes Svevo and Chinese Spring showed the highest levels of DIMBOA on day 11, which gradually declined until day 18. In Chinese Spring, 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. Overall, the two-way-ANOVA of both DIMBOA and DIM2BOA-Glc demonstrates significant differences in time (day after germination) and in the two wheat genotypes. The HDMBOA-Glc levels were much lower than the other two BXD compounds and only detected in Chinese Spring 11 and 18 days after germination and only at 18 days after germination in Svevo leaves.
Calculating the trichome density
To explore the physical adaptation of wheat plants against aphid invasion, we evaluated the trichome density on the leaf edge and surface as a physical barrier (Figure 6). As presented in Figure 6A, the two-way-ANOVA suggests a significant difference in trichome number on the edges between the three wheat genotypes (F genotype (2,224) = 498.36, P value < 0.0001), the day after germination (F time (2,224) = 5.30, P value = 0.0056), and a cross-effect (F genotype*time (2,224) = 3.05, P value = 0.0178). The highest number of trichomes on the leaf edges was in Zavitan, while Svevo and Chinese Spring had a similar number of trichomes, except 18 days after germination, when Chinese Spring showed the lowest trichome number. The number of trichomes slightly increased over time, mainly on Zavitan and Svevo. In Figure 6B, the trichome density on the leaf surface is presented. The two-way-ANOVA suggests a significant difference in trichome number on the surface between the three wheat genotypes (F genotype (2,372) = 268.32, P value < 0.0001), the day after germination (F time (2,372) = 15.99, P value < 0.0001), and a cross-effect (F genotype*time (2,372) = 8.70, P value < 0.0001). The highest number of trichomes on the leaf surface was in Zavitan, while Chinese Spring possessed the lowest trichome number. Images of the trichomes demonstrated that the trichomes are also different in their lengths and angles. As shown in Figure 6C, the trichome lengths observed on the leaf surfaces of Zavitan and Chinese Spring were longer than those on Svevo. Additionally, the trichomes on the edges faced one direction only on Zavitan and Chinese Spring, while on Svevo, they faced both directions. Altogether, this suggested that Zavitan has more trichomes as a physical barrier than the other two genotypes.
Evaluation of aphid reproduction on different wheat genotypes
The aphid performance on wheat seedlings infested with aphids for 96 h was evaluated at two time points: 15 and 18 days after germination. The two-way ANOVA suggested a significant difference in the aphid progeny between the three wheat genotypes (F genotype (2,45) = 20.60, P value < 0.0001), the day after germination (F time (1,45) = 161.93, P value < 0.0001), and a cross-effect (F genotype*time (2,45) = 4.24, P value = 0.021), indicating the effect of both genotype and age (day after germination) on aphid reproduction (Figure 7). In the two measurements (15 and 18 days after germination), the Chinese Spring genotype was more aphid-resistant than the other two genotypes, while Zavitan and Svevo were not significantly different. Additionally, 18-day-old wheat seedlings were more susceptible to aphids than the 15-day-old plants.