Comparative transcriptome and metabolic analysis of wild and domesticated wheat genotypes reveals differences in chemical and physical defense responses against aphids

Background: Young wheat plants are continuously exposed to insect herbivorous attack. To maintain their fitness, plants have evolved different defense mechanisms, including the biosynthesis of deterrent compounds named benzoxazinoids, and/or trichome formation that provides physical barriers. It is unclear whether both of these mechanisms are equally critical in providing an efficient defense for wheat seedlings against the bird cherry-oat ( Rhopalosiphum padi ) aphid—an economically costly pest in cereal production. Results: In this study, we compared the transcriptomic, metabolomic, chemical, and physical defenses of three selected wheat genotypes to aphid performance. We chose diverse wheat genotypes, two tetraploid wheat genotypes, domesticated durum ‘Svevo’ and wild emmer ‘Zavitan,’ and one hexaploid bread wheat, ‘Chinese Spring.’ The full transcriptomic analysis revealed a major difference between the three genotypes, while the clustering of significantly different genes suggested a higher similarity between the two domesticated than the wild wheat. A pathway enrichment analysis indicated that genes associated with primary metabolism, as well as the pathways associated with defense such as phytohormones and specialized metabolites, were altered between the three genotypes. Measurement of benzoxazinoid levels at the three time-points (11, 15 and 18 days-after-germination) revealed high abundance levels in the two domesticated genotypes, while the levels were very low in the wild emmer wheat. The Chinese Spring showed a more diverse benzoxazinoid (known and putative) composition than the other two genotypes. In contrast to the benzoxazinoid levels, the trichome density was dramatically higher in the wild emmer than in the domesticated wheat. Evaluation of aphid reproduction indicated that the domesticated bread wheat is more resistant than the tetraploid genotypes. Conclusions: We compared the amount of benzoxazinoids, the trichome number, and aphid reproduction at three time-points, as well as performing a tryptophan decarboxylase the resistance against cereal cyst altering the downstream secondary metabolite and germination selected during wheat domestication are co-located on


Abstract
Background: Young wheat plants are continuously exposed to insect herbivorous attack.
To maintain their fitness, plants have evolved different defense mechanisms, including the biosynthesis of deterrent compounds named benzoxazinoids, and/or trichome formation that provides physical barriers. It is unclear whether both of these mechanisms are equally critical in providing an efficient defense for wheat seedlings against the bird cherry-oat ( Rhopalosiphum padi ) aphid-an economically costly pest in cereal production. Results: In this study, we compared the transcriptomic, metabolomic, chemical, and physical defenses of three selected wheat genotypes to aphid performance.
We chose diverse wheat genotypes, two tetraploid wheat genotypes, domesticated durum 'Svevo' and wild emmer 'Zavitan,' and one hexaploid bread wheat, 'Chinese Spring.' The full transcriptomic analysis revealed a major difference between the three genotypes, while the clustering of significantly different genes suggested a higher similarity between the two domesticated than the wild wheat. A pathway enrichment analysis indicated that genes associated with primary metabolism, as well as the pathways associated with defense such as phytohormones and specialized metabolites, were altered between the three genotypes. Measurement of benzoxazinoid levels at the three time-points (11,15 and 18 days-after-germination) revealed high abundance levels in the two domesticated genotypes, while the levels were very low in the wild emmer wheat. The Chinese Spring showed a more diverse benzoxazinoid (known and putative) composition than the other two genotypes. In contrast to the benzoxazinoid levels, the trichome density was dramatically higher in the wild emmer than in the domesticated wheat. Evaluation of aphid reproduction indicated that the domesticated bread wheat is more resistant than the tetraploid genotypes. Conclusions: We compared the amount of benzoxazinoids, the trichome number, and aphid reproduction at three time-points, as well as performing a 3 transcriptome analysis. Overall, the results suggested that although wheat seedlings possess both chemical and physical defenses, the chemical defense plays a more significant defensive role than the physical barriers.

Background
Crop plants are continually exposed to different environmental stress conditions, such as herbivore infestation [1]. Insect herbivory is a crucial factor in yield loss and quality degradation in agricultural crop production. Average losses can reach 20-30% of yield, and in some cases, they can cause a total yield loss [2]. One of the most economically significant herbivores is found in the aphid family (Hemiptera: Aphididae) [3], a piercingsucking pest that feeds on the phloem sieve. Aphid infestation causes direct damage by consumption of water and nutrients, as well as indirect damage by plant virus transmission [4][5][6]. In response to insect infestation, plants produce constitutive and inducible defenses to reduce damage and enhance their own fitness [7]. Although many plant defenses are produced constitutively during a specific developmental stage, regardless of insect attack, others are inducible in response to insect damage. Examples of herbivore-induced defense strategies are the accumulation of chemical defenses such as benzoxazinoids, glucosinolates, and alkaloids, which are classes of specialized metabolites that function as deterrents. Another strategy is mechanical defense, including physical barriers such as the increased density of thorns, spikes, or trichomes [8][9][10][11][12].
Most of the defense strategies are abundant in young seedlings and decrease during development toward the juvenile stage [13][14][15]. The herbivore-induced strategies are mediated by the modification of signaling (i.e., jasmonic and salicylic acid) [16], which allows the plants to conserve metabolic resources and energy to be directed toward growth and reproduction in the absence of insect herbivory.
Wheat is a staple crop that provides 20% of the world population's caloric and protein 4 intake [17]. It was first domesticated more than 10,000 years ago, making it one of the earliest domesticated crops [18]. The process of domestication, centuries of cultivation, and modern wheat breeding have led to the reduction or narrowing of genetic variation compared to their wild ancestors [19]. This reduction is due to the small initial crop population, coupled with intense selection for agronomic traits, considered as a ''domestication bottleneck'' [20,21]. Adaptation of domesticated wheat and wheat varieties to local conditions intensified the reduction, giving rise to landrace cultivars [22]. Moreover, cultivation of a germplasm with a narrow genetic base entails a risk due to genetic vulnerability to biotic and abiotic stresses, possibly resulting in severe crop losses. Using wild relatives of wheat as a proxy, the genetic diversity of agriculturally important traits can be contrasted before and after domestication [20]. It has been assumed that several agriculturally important traits, such as biotic and abiotic stress resistance, significantly decreased during wheat domestication [23,24]. This lowered resistance, reported for responses to herbivore attack [25,26], bacterial blight [27], and fungal disease [28], was revealed in domesticated members of the plant families Fabaceae and Brassicaceae [24]. In contrast, increased resistance of domesticated wheat to aphids was also reported [29]. It was suggested that the tuning of plant domestication defense mechanisms is dependent on pest feeding habits [30].
In this work, we investigated the differences in the transcriptome and metabolome of wheat seedlings, and we studied the effect on the chemical and physical defenses against Rhopalosiphum padi aphids. The variation between Triticum species for cereal aphid resistance has demonstrated the potential of using tetraploid wheat to reveal plant defense mechanisms [31][32][33]. Therefore, we focused our research on three representative wheat genotypes: i) Svevo, a tetraploid durum wheat cultivar (Triticum turgidum ssp. durum), ii) a wild emmer, Zavitan, a tetraploid ancestor of modern domesticated tetraploids (Triticum turgidum ssp. dicoccoides); both tetraploids have been intensively investigated as potential sources for resistance genes and markers [29,34]; and iii) the spring wheat genotype Chinese Spring, which is widely used as a reference genome and for cytogenetic studies [19,[35][36][37][38][39][40]. We focused on young seedlings (11-18 days old) which produce high levels of chemical and physical defenses. The chemical deterrent metabolites (benzoxazinoids) and the physical defenses (trichome) of tetraploid and hexaploid wheat species were analyzed and compared to R. padi aphid reproduction. The benzoxazinoids were highly abundant in the domesticated wheat while trichomes were highly present in the wild emmer. Comparing the effect of aphid reproduction suggested that the chemical defense plays a more significant defensive role than the physical barriers. Altogether, this study evaluates the dynamics of defense mechanisms in response to aphid attack.

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 6 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 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; |log 2 (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 7 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 (13lipoxygenase; 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 8 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 (H 2 O 2 ) 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][52][53][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][59][60][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 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][14][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 (F 2,224 = 498.36, P value < 0.0001 and F 2,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 (F 2,224 = 5.30, P value <0.0056 and F 2,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 14day-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 (F 2,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 (F 1,59 = 168.03, P value <0.0001). Overall, 13 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].

Discussion
In this research, we selected three wheat genotypes, including wild emmer, durum, and bread wheat, to address the fundamental question of the effect of domestication on plant resistance against aphids. To understand the overall gene levels and the differential gene expression between genotypes [71], we compared the transcriptome of 11-day-old seedlings of the three genotypes. While the PCA plots relying on the transcriptomic or untargeted metabolite data suggested a unique pattern for each genotype (Figure 1), the heatmap of the differentially expressed genes (8,735 unique transcripts) indicated a higher similarity between the domesticated genotypes than the wild emmer wheat ( Figure   S2). This pattern was similar when we compared the gene expression of specific biosynthetic genes of the benzoxazinoid pathway ( Figure 4). The similarity between Svevo and Chinese Spring is supported by a recent study that compared the exome of approximately 500 wheat genotypes, and suggested that the wild emmer wheat is the progenitor of the A and B subgenomes of all the modern tetraploid and hexaploid genotypes. The T. durum lineage was found to be the most likely ancestor of the bread wheat cultivated germplasm [19]. The similarity between Svevo and Chinese Spring might be due to a "domestication bottleneck" [20][21][22].
A previous study that explored the variation between 19 hexaploid bread wheat pangenomes reported that genes involved in the response to environmental stress and defense against biotic stress were variable between the genomes [72]. Similarly, in our 14 results, several enriched functional genes involved in the response to environmental stress and defense responses are variable between the three genotypes, including polyamine biosynthesis and 13-lipoxygenase, which are the first enzymatic steps of the phytohormone, jasmonic acid (Table 1). Additionally, the over-representation pathway enrichment analysis demonstrated that genes associated with amino acid metabolism and the biosynthesis of carbohydrates, cell structures, fatty acids and lipids, phytohormones, and specialized metabolites were significantly different between the three genotypes.
These observations indicated the major differences in the basal gene expressions between the genotypes. This was supported by several experiments where we measured the water content, and primary metabolites of the 11-day-old leaves ( Table 2). The main results suggested that Zavitan has higher levels of amino acids, organic acids, and sugars than the other two genotypes (Figure 8). Interestingly, the total chlorophyll, and chlorophyll a and b levels were positively correlated with the amount of DIMBOA measured at 11-days ( Figure 8). This might be supported by previous studies that revealed additional role of BXDs as iron chelators in the roots [73,74].
A recent study suggested that domesticated wheat has maintained its defense traits against specialized herbivores that have coexisted with the crop throughout its domestication, but that it is less efficient against generalist herbivores [30]. We compared the reproduction of the R. padi aphid on the three wheat genotypes (Figure 7), which is among the most economically important aphid associated with a host range of well over 100 species [75,76]. Our results show that the wild emmer, which possessed a very low amount of chemical defenses and a higher number of trichomes, is more susceptible to aphids than the cultivated genotypes ( Figures 5-7). In contrast to our results, a previous study by Migui and Lamb (2003) tested 41 accessions of wild and cultivated wheat in the field for resistance to three species of aphids, including R. padi. The results showed the highest aphid resistance was in the diploid species and the lowest was in the hexaploid species [32]. This discordance might arise from accession-specific differences or may be due to the performing of the experiment on mature wheat plants in field conditions, compared to young seedlings under controlled conditions ( Figure S1). This was supported by another study that suggests that the differences in aphid preference depend on the plant developmental stage in the field compared with seedlings in the laboratory [77].
In wheat, genes encoding defense mechanisms are found in hexaploid bread wheat (genome BBAADD), tetraploid wheat (genome BBAA), and in the three diploid progenitors of hexaploid wheat (genomes AA, BB, and DD) [15,78], while several of the enzyme steps are not yet identified. A comparison between the transcript levels from the three subgenomes in hexaploid wheat indicated that the homeologs on the B subgenome are the main contributors to the benzoxazinoid biosynthesis pathway, especially in shoots [15,79].
It was suggested that gene expression in diploids and tetraploids are far more varied than in the hexaploid population, indicating that the former species likely encode a greater allelic variation that, in turn, could facilitate breeding for pest resistance [80]. In this research, Bx gene expression was detected in all three genotypes, while the two domesticated wheat genotypes showed a closer pattern of gene expression and BXD levels than the wild emmer ( Figure 4A).
The BXDs are a diverse class of specialized metabolites, which are mainly known for their deterrent functions. They are known to have a crucial effect on plant resistance to insects such as aphids [67, 81,82], chewing herbivores [83][84][85], fungal infection [86], and other insects, diseases and weeds [69]. However, a recent study reported that the molecular functions of the BXDs are diverse beyond toxicity. It was suggested that the BXDs function in shaping the root microbiome by selectively attracting the plant-beneficial rhizobacterium Pseudomonas putida [87,88], and are also used as iron chelators [89].
Also, it has been estimated that naturally occurring DIMBOA may govern a physical defense mechanism against aphid feeding by the accumulation of callose [90].
Furthermore, it was demonstrated that in maize, plant-aphid infestation caused secretion into the apoplast of DIMBOA-Glc, but not HDMBOA-Glc [67], which may indicate a unique function of DIMBOA-Glc. Our results indicated that the domesticated bread wheat Chinese Spring, which showed the lowest amount of aphid progeny, possesses a varied array of BXDs including known compounds (DIMBOA and DIM2BOA-Glc), and five other unknown BXDs ( Figure 5). This suggested that some of the unknown BXDs play a defensive role, which requires further investigation.
BXD abundance in wheat seedlings decreased during development; these metabolites show the highest activity in plants at the young stage [13][14][15], and tend to decline as the plant ages [91]. Our results demonstrate that aphid progeny numbers increased over time, while the BXD levels decreased (Figures 5, 7). The leaf surface is the first defense barrier against insects; it includes trichomes, thorns, silica, and wax [92,93]. Trichomes are used as a physical defense by disturbing herbivore movement, development, and oviposition, while glandular trichomes are used for exudate storage and secretion [94]. The surface density of trichomes and their dispersal on the leaf edges were not directly proportional to the number of aphid progeny. This suggests higher effectiveness of the chemical defenses rather than the defensive physical barrier. Additionally, the metabolic profile and water content indicated that 11-day-old Zavitan seedlings have higher water content, as well as higher levels of some amino acids, organic acids, and sugars, than the cultivated wheats.
This lack of chemical defense and better water and nutrients status might be beneficial for aphid reproduction on Zavitan relative to the domesticated wheats.

Conclusions
Our results suggest that benzoxazinoids provide a better defense mechanism than trichomes against R. padi aphids. Comparisons of significant gene expression, phenotypic characterization, and the chemical defense and physical responses indicated a higher similarity between the domesticated wheats than between either of them and the wild emmer. This suggested that under insect pressure, wheat plants might have undergone evolutionary convergence, which results in similarities of defense mechanisms by the biosynthesis of defense metabolites and, to a lesser extent, trichome formation.

Plant growth condition
Wheat seedlings were grown on moistened HR2 (soil mix). Wheat seeds were planted 1.5-

Aphid bioassays
A bird cherry-oat aphid (Rhopalosiphum padi) was collected in the field, identified and colony was maintained on two-week-old wheat plants, Triticum aestivum variety named Rotem (Rotem seeds were obtained from Agridera Seeds & Agriculture LTD, Israel). the colony was grown under a 16-hr-light/8-h-dark photoperiod at a constant room temperature of 25 0 C. Apterous adult R. padi aphids were used for the experiments. The aphids were confined to the upper part of the second leaf of a seedling for 96 h using clip cages (4.5 cm in diameter), and empty cages without aphids were used as a control. The number of aphids was counted after 96 h of infestation, and tissue samples were harvested at the same time.

RNA extraction, transcriptome sequencing, and analysis
Tissue from 11-day-old plants was combined into one experimental replicate, and three replicates were collected for each genotype. Total RNA was extracted using an SV Total RNA Isolation Kit with on-column DNaseI treatment (QIAGEN). The purified RNA was quantified, and 2.5 mg of each was used for next-generation sequencing using an Illumina HiSeq 4000 instrument with a 150 bp paired-end read length conducted by the GeneWIZ Company (www.genewiz.com). Quality control was performed using FASTQC, and adapters, and low-quality sequences were trimmed and removed using Trimmomatic v0.36. Mapping was performed using a STAR aligner v2.5.2b against the Triticum aestivum reference genome v1.1 [35]. Reads aligning to exons were retrieved using Subread v1.5.2. For differential gene analysis, DESeq2 v1.22.2 [100] was used to perform a likelihood ratio test to evaluate multiple genotypes at once (adjusted p-value < 0.05). The IWGSC refseq v1.0 high confidence (HC) functional annotation was used without sequences from the D genome, as this genome is not present in the wild and cultivated tetraploid wheat, and from the unknown genome. Next, DESeq2 was performed on the statistically significant genes from the LRT, and a regularized log was applied to the results. MetGenMAP was used to perform a pathway enrichment analysis [42], using rice ortholog IDs, which were converted from Phytozome wheat transcript IDs. In order to convert between IWGSC and Phytozome wheat transcript IDs, a reciprocal BLASTp comparison was performed, and only transcript IDs with mutual hits were retrieved.

Benzoxazinoid metabolite analyses using high-performance liquid chromatography (HPLC)
For measuring the benzoxazinoid levels, we used the HPLC analysis as previously described [70,101,102]. In brief, the wheat tissue within the clip cage was collected, weighed and extracted in 80% methanol, 19.9% DDW and 0.1% formic acid, with a 1:10 (w:v) ratio. The mix was vortexed briefly; then, the tubes were shaken for 40 min at 4 °C, and centrifuged for 5 min at 14,000 g. A filtration was performed by centrifuging the samples on a 22-μm filter plate (EMD Millipore Corp., Billerica, MA, USA) at 3,000 g for 5 min. The extracted mixtures were covered with a WebSeal Mat and kept at 10 °C. Samples were injected into a DIONEX UltiMate 3000 HPLC system using a C18 reverse-phase Hypersil GOLD column (3 µm pore size, 150x4.60 mm; Thermo Fisher Scientific, Germany).
The column oven temperature was 40 0 C and UV-VIS absorbance spectra at 266 nm. For metabolite separation, a gradient of 0.1% formic acid in DDW (solvent A) and 0.1% formic acid in acetonitrile (solvent B) with a flow rate of 1 ml x min -1 was used. The following linear gradient was used: 0 to 10 min, gradient from 5-55% B; 10 to 13 min, gradient from 55-100% B; 13-14 min, gradient from 100-5%; and 14-18 min, 5% solvent B. For benzoxazinoid identification, we compared the chromatograms with the authentic standards of DIMBOA and DIBOA (Toronto Research Chemicals, Toronto, Canada) and DIMBOA-Glc, DIM2BOA-Glc; HDMBOA-Glc, and HDM2BOA-Glc from a plant crude extract (Table S9). UV unique spectra were used to confirm the known BXD in addition to five unknown BXD metabolites ( Figure S4). Chromeleon software (Thermo Fisher Scientific Inc.) version 7.2 was used for system control and data acquisition. During the sample running, the mobile phase consisted of 95% water: 5% acetonitrile: 0.1% formic acid, and 0.1% formic acid in acetonitrile [103]. MassLynx software (Waters) version 4.1 was used for system control and data acquisition (Table S3).

Metabolic analysis using gas chromatography-mass spectrometry (GC-MS)
Metabolites were extracted using approximately 200 mg of ground frozen plant tissue mixed with solvents containing a ratio of methanol/water/chloroform of 55:23:22, v/v following the protocol as previously described [104]. After phase separation, 100 μl of the top hydrophilic layer was collected and dried in a vacuum. Samples were derivatized by adding 40 μl of 20 mg/ml metoxyamine hydrochloride (Sigma-Aldrich) dissolved in pyridine following incubation for 2 h in an orbital shaker at 1,000 rpm at 37 °C. Next, N-methyl-N-(trimethylsilyl) tri-fluoroacetamide (MSTFA) including standard mix (Alkanes) in a volume of 77 μl was added to each sample followed incubation for 30 min in an orbital shaker at 37 °C. For the last step, 1 μl of the sample was injected into the Agilent 5977B GC-MS instrument, and data acquisition was conducted by Mass Hunter software, and the NIST mass spectral library normalized to D-sorbitol 13 C 6 (Sigma-Aldrich) as an internal standard and normalized to the sample fresh weight.

Trichome density
The upper parts of second leaves (same as used for applying clip cages) were collected, and chlorophyll was bleached by 70% ethanol at 85 0 C for 8 min, then rinsed with water.
The tissue was placed on glass microscope slides facing to the adaxial side (leaf surface up). Trichome density images were photographed with a digital camera (Axiocam 305 color) connected to a Zeiss Axioplan 2 Upright Light Microscope (Zeiss, Oberkochen, Germany). For each leaf, nine photos were taken including three from each side (left, right), and three from the middle (top, medium, and bottom positions). We counted the number of trichomes on the edges per mm and the density in mm 2 using ImageJ software (https://imagej.nih.gov/ij/).

Water content in the leaves
The second leaf (10 cm of tissue from the leaf tip) of 11-day-old wheat seedlings was collected. The total fresh weight was measured, and then the samples were placed in 2-3 ml of 5 mM CaCl 2 solution for 8 h followed by drying in a 60 0 C oven for 3-4 days followed by measuring the dry weight. The calculation of the relative water content of the leaf was previously described [105].

Detection of hydrogen peroxide
To examine the basal level of hydrogen peroxide in wheat leaves, a 3,3'-diaminobenzidine (DAB) staining was used for in situ detection [106]. In accordance with the aforementioned protocol, the second leaves of 11-day-old wheat seedlings were gently vacuum-infiltrated with DAB solution. As a control treatment, 10 mM of sodium phosphate was applied to replicate leaves. Following vacuuming, samples were incubated in DAB solution for 4 h, which was then replaced with bleaching solution (ethanol: acetic acid: glycerol 3:1:1) to remove the chlorophyll and to visualize the precipitate formed by hydrogen peroxide 22 (renders precipitates in dark brown).

Chlorophyll content
Fresh leaf tissues (50 mg) were incubated in 5 ml of ice-cold 80% acetone for 48 h, centrifuged at 5,000 × g for 5 min, and absorbance was recorded at 663 nm (chlorophyll a) and 645 nm (chlorophyll b). The amount of chlorophyll was calculated following the procedure of Arnon (1949) and expressed in mg g -1 FW [48].

Statistical analysis
For the principal component analysis (PCA) plot, the data were normalized as previously described [107], and graphs were plotted using MetaboAnalyst 3.0 software [108]. The one-and two-way ANOVAs (analysis of variance) used JMP software (SAS; www.jmp.com) and Microsoft Excel for figure representation. In order to test the three genotype groups, we performed a likelihood ratio test (LRT) using DESeq2. LRT is a statistical test, similar to ANOVA, which allows the comparison of all levels of a factor at once. We estimated the cluster number (K) of the results using clusGap [41] and performed K-means clustering with the k-means base function. The resulting gene clusters were evaluated for over-and under-representation with agriGO v2 [43].

Availability of data and material
The datasets used and/or analyzed during the current study are available in this manuscript in the supplementary files or from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.

This research was supported by the Binational Agricultural Research and Development
Fund (BARD grant number IS-5092-18R). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.  Table 1. Enrichment analysis of metabolic pathways grouped by k-means clustering. Gene expression patterns were sorted into eight clusters, as determined by a k-means analysis of transcripts detected in the three wheat genotypes.  Overview of gene expression clusters calculated by k-means clustering. A oneway ANOVA was performed, which revealed a total of 8,735 unique transcripts with significant expression profile changes for at least one genotype. The total 41 number of transcripts in each cluster is indicated, and the data for individual genes are shown in light gray. In black, the expression responses for each cluster are shown in black. Heatmap of genes from two selected pathways: A) benzoxazinoids and B) trichome formation. 44 45 Figure 5 The BXD levels in three selected wheat genotypes. Leaves were harvested 11, 15 and 18 days after germination, and the samples were analyzed using HUPLC (mean ± SE, n = 5-8). ND, not detected. Asterisks above the bars indicate significant differences (Statistics Student's t-test (P value < 0.05). A) 2,4dihydroxy-7-methoxy-1,4-benzoxazin (DIMBOA), B) 2,4-dihydroxy-7,8-dimethoxy-1,4-benzoxazin-3-one glucoside (DIM2BOA-Glc), C) Total unknown BXD. Trichomes on wheat seedlings. Bars represent the average number of edge trichomes (A) and density ± SE of surface trichomes (B) (mean +/-SE, n = 6-8).
Different letters above the bars indicate significant differences, P value < 0.05, ANOVA followed by Tukey's HSD test. Aphid reproduction on wheat seedlings. The total number of R. padi aphids was counted and calculated per apterous adult aphids after 96 h of infestation. Bars represent the average number of aphid reproduction per adult (mean ± SE, n = 5-21). Different letters above the bars indicate significant differences, P value <0.05, one-way analysis of variance (ANOVA) followed by Tukey's HSD test.   Figure