Abstract High nitrogen use efficiency (NUE) in bread wheat is pivotal to sustain high productivity. Knowledge about the physiological and transcriptomic changes that regulate NUE, in particular how plants cope with nitrogen (N) stress during flowering and the grain filling period, is crucial in achieving high NUE. Nitrogen response is differentially manifested in different tissues and shows significant genetic variability. A comparative transcriptome study was carried out using RNA-seq analysis to investigate the effect of nitrogen levels on gene expression at 0 days post anthesis (0 DPA) and 10 DPA in second leaf and grain tissues of three Australian wheat (Triticum aestivum) varieties that were known to have varying NUEs. A total of 12,344 differentially expressed genes (DEGs) were identified under nitrogen stress where down-regulated DEGs were predominantly associated with carbohydrate metabolic process, photosynthesis, light-harvesting, and defense response, whereas the up-regulated DEGs were associated with nucleotide metabolism, proteolysis, and transmembrane transport under nitrogen stress. Protein–protein interaction and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis further revealed that highly interacted down-regulated DEGs were involved in light-harvesting and photosynthesis, and up-regulated DEGs were mostly involved in steroid biosynthesis under N stress. The common down-regulated genes across the cultivars included photosystem II 10 kDa polypeptide family proteins, plant protein 1589 of uncharacterized protein function, etc., whereas common up-regulated genes included glutamate carboxypeptidase 2, placenta-specific8 (PLAC8) family protein, and a sulfate transporter. On the other hand, high NUE cultivar Mace responded to nitrogen stress by down-regulation of a stress-related gene annotated as beta-1,3-endoglucanase and pathogenesis-related protein (PR-4, PR-1) and up-regulation of MYB/SANT domain-containing RADIALIS (RAD)-like transcription factors. The medium NUE cultivar Spitfire and low NUE cultivar Volcani demonstrated strong down-regulation of Photosystem II 10 kDa polypeptide family protein and predominant up-regulation of 11S globulin seed storage protein 2 and protein transport protein Sec61 subunit gamma. In grain tissue, most of the DEGs were related to nitrogen metabolism and proteolysis. The DEGs with high abundance in high NUE cultivar can be good candidates to develop nitrogen stress-tolerant variety with improved NUE. Keywords: transcriptomics, nitrogen use efficiency, Australian wheat cultivars, nitrogen stress, RNA-seq Introduction Over the past several decades, application of nitrogen fertilizer has been a practiced way to gain optimal crop yield. N fertilizer usage is predicted to reach 105 Tg N by 2030 and 135 Tg N by 2050 ([37]Good et al., 2004). However, overuse of fertilizers can cause significant environmental issues such as erosion, soil quality depletion, and contamination of water supplies at local, regional, and global scales ([38]Ahmad et al., 2008; [39]Guo et al., 2010). Thus, it is important to develop new varieties with high nitrogen use efficiency (NUE). A better understanding of gene expression and regulation under nitrogen stressed conditions will help achieve this goal. Response to nitrogen scarcity in plants is controlled by changes in gene expression involved in different molecular mechanisms that are mainly related to plant developmental processes and yield ([40]Zhang et al., 2006; [41]Kant et al., 2011). In particular, wheat grain production largely depends on the provision of N fertilizer and cultivars with high N uptake and utilization efficiency ([42]Nyikako et al., 2014; [43]Garnett et al., 2015; [44]Cormier et al., 2016; [45]Hitz et al., 2017). The biological pathways related to NUE are known to be strongly influenced by genetic variation as well as environmental factors such as N availability ([46]Moll et al., 1982; [47]Xu et al., 2012). Studies showed that N limitation can negatively affect wheat growth, morphology, and agronomic traits ([48]Chandna and Ahmad, 2015; [49]Curci et al., 2017; [50]Wen et al., 2018; [51]Wang J. et al., 2019). Identifying key genes to improve stress tolerance in low N conditions is a feasible way to raise NUE. It is important to select cultivars that have contrasting NUEs for a comparative understanding of gene expression and regulation in response to N stressed conditions ([52]Hirel et al., 2007; [53]Kant et al., 2010). There are a number of approaches that have been undertaken by researchers to unravel how plants adapt to stressed conditions ([54]Shrawat et al., 2008). In recent years, next-generation sequencing techniques have provided opportunities to study the gene expression and their regulations at the transcriptome level, and they have significantly enhanced the success rate of gene discovery ([55]Diao et al., 2019). A number of studies also reported on transcriptome profiling by using Illumina’s RNA-sequencing ([56]Dai et al., 2015). Most of the studies demonstrated how a single genotype performed using contrasting environmental and growth conditions. In Arabidopsis, N response-related genes were identified using microarray analysis of gene expression changes in response to short-term and long-term treatments for nitrate with different concentrations ([57]Wang et al., 2001; [58]Price et al., 2004). Likewise, transcriptome study on different tissues with short-term N stress in rice also revealed a significant number of N responsive genes ([59]Lian et al., 2006). Transcriptome study on long-term N stress was also reported in rice ([60]Ym et al., 2009). However, a comprehensive transcriptome investigation by combining contrasting tissue, developmental stage, genotype, and N treatment is still lacking. Nitrogen stress has a significant impact on the overall plant physiological process ([61]Zhao et al., 2005) related to plant height, dry matter, grain yield (GY), and grain protein content (GPC) ([62]Baligar et al., 2007; [63]Wang W. et al., 2003). Nitrogen strongly influences photosynthesis through a large deposition of leaf N to ribulose bisphosphate carboxylase/oxygenase (Rubisco) and its involvement in stomatal opening ([64]Evans, 1989). Approximately 75% of leaf N is allocated to chloroplasts, with about 27% of this utilized in Rubisco to ensure high photosynthetic activity ([65]Evans, 1989; [66]Makino, 2003). Nitrogen also influences photosynthesis via its impact on CO[2] assimilation and sugar partitioning ([67]Drew et al., 1989; [68]Foyer et al., 2011; [69]Ishikawa-Sakurai et al., 2014). The decreased photosynthesis ultimately resulted in decreased biomass production and yield ([70]Poorter and Evans, 1998; [71]Long et al., 2006; [72]Jin et al., 2015). The regulation of plant photosynthetic activity is reported to be associated with brassinosteroids (BRs), a class of steroid hormones ([73]Sakamoto et al., 2006; [74]Komatsu et al., 2010). BRs are known to regulate stress responses and play important roles in regulating plant growth and development ([75]Wang et al., 2012; [76]Zhao and Li, 2012; [77]Hayes, 2019). Several studies in Arabidopsis and rice showed involvement of BRs in controlling flowering, leaf senescence, chloroplast development, plant height, tiller numbers, and biomass, which are important agronomic traits affecting GY ([78]Chono et al., 2003; [79]Mussig et al., 2003; [80]Sakamoto et al., 2006; [81]Wu et al., 2008; [82]Jeong et al., 2010). In wheat, BRs were also reported to be involved in promoting root growth and water stress tolerance ([83]Hayes, 2019; [84]Hou et al., 2019). However, the correlation of N stress on steroid biosynthesis has not been well studied. Thus, response to N stress is a rather complex process, and a better understanding of genes involved in different pathways is needed to develop stress-tolerant wheat varieties. This study investigated three Australian bread wheat varieties, Mace, Spitfire, and Volcani, which are known to have high, medium, and low NUEs, respectively ([85]Alhabbar et al., 2018b). Since gene expression in plants is controlled in a temporal and tissue-specific manner ([86]Koltunow et al., 1990; [87]Maizel and Weigel, 2004) and the N demand is subject to plant developmental stages, the current study used different tissues at different growth stages to unravel the broad picture of transcriptome profile with the objectives of identifying novel genes that are differentially expressed under long-term N stress compared to high N treatment, and by characterizing the underlying physiological and molecular mechanisms of tolerance to N stress. Materials and Methods Plant Material, Growth Conditions, and Sample Collection Three Australian wheat cultivars, Mace, Spitfire, and Volcani, were used in this study. Plants were grown in a glasshouse with a complete randomized block design (RCBD) including three replicates and using pots (190 mm height × 200 mm top diameter × 180 mm bottom diameter) without holes to avoid leaching. Plants were grown under controlled temperature and sunlight conditions [20/11°C (day/night)] for an 8 h light and 16 h dark photoperiod. The pots were watered manually based on soil water capacity. All plants were supplied with a basal nitrogen dose of 25 kg N ha^–1 after 1 week of sowing. Nitrogen-free Hoagland solution^[88]1 was applied to all plants once every 2 weeks to meet the nutrient demand of plants except N. Two N rates—low (LN)/0 kg N ha^–1 and high (HN)/100 kg N ha^–1—were applied at mid-tillering (Zadoks scale 22–25) and booting (Z43–Z45) stages for plants considered as low and high N treated, respectively. The timing for N applications was adjusted according to Zadoks (Z) decimal growth stage for wheat. Flexi-N (containing 50% urea, 25% nitrate, and 25% ammonium) was used as a source of N because of its high N content (42.2% N). Flexi-N was used since it contains nitrate that is directly available to plants while the urea and the ammonium become available more slowly, enabling a controlled release of N over an extended period ([89]CSBP, 2012). Times for N application, recording of flowering time, measurement of chlorophyll content, and tissue collection were adjusted according to each cultivar’s growth stage. For RNA extraction, the whole flag and second leaf samples were collected at the start of the flowering [0 days post anthesis (DPA)], 10, 20, and 30 DPA, while the developing grains were collected at 10, 20, and 30 DPA from the middle section of the main head, then snap-frozen in liquid nitrogen, and then stored at −80°C for later RNA extraction. Anthesis dates were estimated by the appearance of anthers on approximately 50% of all heads. Plant height was measured from soil surface to the top of the plant, and peduncle length was measured from the peduncle bottom to the joint with the stem. Chlorophyll content was measured using a handheld chlorophyll meter (IC-CCM-200—Chlorophyll Concentration Meter CCM-200 plus). One value per plant was taken from the flag leaf and second leaf on the main stem at four different growth stages: flowering (0 DPA), 10, 20, and 30 DPA. Each value was the average of three measurements recorded from the middle of the leaves. The main stem of each plant was individually labeled to ensure the same leaves were always measured. All plants in a pot (main stem plus tillers) were hand-harvested to measure yield components and the head number per plant counted. The heads were cut off and the seed number per head was counted. Grain samples were oven-dried in a forced-air circulating dryer at 60°C for 72 h. GPC was measured by near-infrared reflectance (NIR) spectroscopy using a FOSS NIR Systems model 5000 spinning cup. NIR data collection used DPIRD wheat calibrations developed over many years with the WinISI software (FOSS NIR Systems Inc., Laurel, MD, United States). RNA Isolation, Library Construction, and Sequencing Leaf and grain samples from three biological replicates were ground in liquid nitrogen, and the total RNA was extracted using a pre-chilled Trizol reagent (Invitrogen, Carlsbad, CA) following the manufacturer’s directions, with some modifications. Proteins were removed with a protein extraction buffer (1 M Tris–HCl, 5 M NaCl, 10% SDS, 0.125 M EDTA, and 1 M DTT). After the protein removal, the acid phenol/chloroform/isopropanol (49:49:2), Trizol, and chloroform were added sequentially for the extraction of total RNA. Isopropanol was used for the precipitation of total RNA, which was subsequently treated with the Qiagen DNase kit to remove potential genomic DNA contamination. Concentration and purity were checked by Nanodrop, with 260/280 absorbance ratios of approximately 2.0, and the degradation and potential contamination was tested by agarose gel electrophoresis. RNA integrity was confirmed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). The mRNA was enriched using oligo (dT) beads and then fragmented randomly in a fragmentation buffer, followed by cDNA synthesis using random hexamers and reverse transcriptase. After first-strand synthesis, a custom second-strand synthesis buffer (Illumina) was added together with dNTPs, RNase H, and Escherichia coli polymerase I to generate the second strand by nick-translation. The final cDNA library was ready after a round of purification, terminal repair, A-tailing, ligation of sequencing adapters, size selection, and PCR enrichment. Library concentration was first estimated using a Qubit 2.0 fluorometer (Life Technologies) and then diluted to 1 ng μl^–1 before checking the insert size on an Agilent 2100 Bioanalyzer. The concentration was then quantified at greater accuracy by quantitative PCR (Q-PCR) (library activity >2 nM). Each library with an individual barcode was sequenced by Illumina HiSeqTM PE125/PE150 (Illumina Inc., United States). Transcriptome Analysis A total of 90 different samples, including 30 each from three cultivars, Mace, Spitfire, and Volcani, were used for RNA-seq analysis. The samples were subjected to low and high nitrogen treatments to study a broad range of cell responses under nitrogen stress. For both treatment conditions, the replicates showed a high correlation coefficient (r > 0.8) between samples. For the RNA-seq downstream analysis, three samples (VAScLNR3, VEScHNR1, and SEScLNR2) were excluded due to sample quality. A total of 2070.85 million raw reads were filtered. A total of 1963.99 million clean reads were aligned against IWGSC RefSeq v1.0 gene models that produce 1750.09 million total mapped reads (TMRs), of which 128.89 million were mapped to multiple sites (MMR) and 1621.21 million were uniquely mapped. Among the TMRs, 810.66 million were mapped with a positive strand and 810.55 million were mapped with a negative strand ([90]Supplementary Table 1). The average leaf Q20, Q30, and GC (Base G + Base C) contents were 96.93, 92.31, and 57.21%, respectively. Similarly, the average grain Q20, Q30, and GC (Base G + Base C) contents were 96.78, 92.16, and 57.79%, respectively. For both tissues, 95% of the total reads were filtered as cleaned reads, which confirms the fine quality of the sequencing results. Approximately, an average of 89% of clean reads were mapped for N-treated leaf samples, whereas 86% were mapped for grain tissue ([91]Supplementary Table 2). For each sample, the percent of reads mapped to exon regions was above 90%, intron reads less than 5%, and intergenic reads less than 3%. The distribution of mapped reads of each sample in chromosome 3B was the highest, while the lowest reads were mapped in chromosome 6A. The gene expression level was measured by calculating the reads mapped to exons. Read count was proportional to the actual expression level as well as to the gene length and the sequencing depth. In order to make comparable gene expression levels estimated from different genes and experiments, fragment per kilobase of transcript per million mapped reads (FPKM) was used for normalization. Considering the influence of various gene lengths and sequencing intensity, FPKM is commonly used to make comparison of gene expression levels among different samples. Analysis of Differentially Expressed Genes For the FPKM, a value of 1.0 was set as the threshold for determining whether a gene was expressed or not. HiSeq v0.6.1 (a Python package for high-throughput sequencing data analysis) was used to analyze gene expression levels in this experiment using the union mode. The correlation between samples was justified by the square of the Pearson correlation coefficient. The DESeq (version 1.10.1, R Bioconductor package) was used to conduct the differential expression analysis. The normalized data were fitted to a negative binomial generalized linear model. The threshold of the p-value after normalization (padj, q-value) was set as ≤0.05 for filtering accurate differentially expressed genes (DEGs). The clustering of DEGs was analyzed based on the FPKM value with the use of ggplot2 (version 2.1.0) and pheatmap (version 1.0.8) ([92]Anders and Huber, 2010, [93]2012; [94]Robinson et al., 2010; [95]Trapnell et al., 2012). The DEGs were identified using the functional annotations of the IWGSC RefSeq v1.0 gene annotation. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis of DEGs Gene ontology (GO) analysis was performed using ShinyGO v0.61^[96]2. GO with a false discovery rate (FDR) corrected at p ≤ 0.05 was regarded as significant enrichment ([97]Young et al., 2010). KOBAS (version 2.0^[98]3, a web server for annotation and identification of enriched pathways and diseases, was applied for Kyoto Encyclopedia of Genes and Genomes (KEGG^[99]4) pathway enrichment analysis. Pathways with an FDR corrected at p ≤ 0.05 were considered as significant enrichment ([100]Mao et al., 2005; [101]Kanehisa et al., 2007). Protein–Protein Interaction Analysis To predict the interaction of DEGs at the protein level under N stress and further confirmation of association of DEGs with biological pathways at the protein level, deduced amino acid sequences of DEGs were used to make a protein–protein network using the STRING (version 11.0) tool^[102]5, a database for known and predicted protein interactions and functional associations predicted in common pathways. Due to the lack of detailed annotation of the wheat protein data available in STRING, we used two well-annotated species, rice, and Arabidopsis, as the reference to get protein–protein interaction information of the homologous wheat proteins. The global network graph of these interactions was constructed using the experimentally evident interacted proteins, and disconnected nodes (proteins) were removed to show the advanced view of highly connected proteins. MCL clustering using the inflation parameter 1.70 was used to show the association of clusters in KEGG pathways. Hierarchical Cluster Analysis Hierarchical cluster analysis was performed using the Morpheus package^[103]6 Complete linkage analysis was performed using the Spearman rank correlation values. Statistical Analysis All data generated from the glasshouse experiments were analyzed by SPSS (version 24). Univariate analysis of variance (UNIANOVA) was used to determine the significance of different factors on each agronomic trait and protein parameter. The significant statistical difference was judged at p ≤ 0.05. Results Agronomic Performance of Wheat Cultivars Under Low and High Nitrogen Conditions Under N stress (0 Kg N/ha), most of the agronomic traits were affected negatively in all three cultivars. In general, tiller number, GY, and chlorophyll content were mostly affected by N stress, whereas flowering days and GPC were less affected. A strong variation in grain weight per plant has been observed, which is considered as a yield component for small-scale glasshouse experiments. Grain weight per plant was dropped by 78% in Mace, 81% in Spitfire, and 80% in Volcani ([104]Figure 1A) due to N stress. Similarly, under N stress, the tiller number ([105]Figure 1B) was decreased by 72.4, 84.2, and 81.2%, and the chlorophyll content of both flag leaf and leaf 2 ([106]Figures 1C,D) were decreased by approximately 85, 80, and 68% for Mace, Spitfire, and Volcani, respectively. In addition, a significant reduction in plant height ([107]Figure 1E), main head length ([108]Figure 1F), and spikelet number per head ([109]Figure 1G) has also been observed under N stress. Flowering days ([110]Figure 1H) were decreased by 3.2, 4.9, and 7.9%, and the GPC ([111]Figure 1I) decreased by 4.1, 9.3, and 29.5% in Mace, Spitfire, and Volcani, respectively. A significant negative impact of low N on leaf area and peduncle length has also been noticed ([112]Figures 1J–L). Overall, the influence of N stress on growth and agronomic traits was variable across the cultivars, where Spitfire and Volcani were more affected compared to Mace. FIGURE 1. [113]FIGURE 1 [114]Open in a new tab Growth and agronomic traits of Triticum aestivum cultivars (Mace, Spitfire and Volcani) under low (0N) and high (100 N) treatments. (A) grain weight/plant (B) tiller numbers/plant (C) flag leaf chlorophyll content (D) leaf 2 chlorophyll content (E) plant height (F) main head length (G) spikelet number/head (H) flowering days (I) grain protein percentage (J) flag leaf area (K) leaf 2 area (L) peduncle length. The values are presented as means ± standard deviation (SD) of three independent biological repeats. Error bars were calculated from three biological replicates and one-way ANOVA was used to test for significance of nitrogen treatment effects on different parameters at P≤0.05 level. ^∗, ^∗∗, ^∗∗∗ Significant at the 0.05, 0.01 and 0.001 probability level, respectively. Overview of RNA-Seq Transcriptome Profile in Response to Nitrogen Stress A total of 12,108 DEGs in leaf tissue and 276 DEGs in grain tissue were identified under the N stressed condition. Mace, Spitfire, and Volcani had 699, 10,535, and 1671 DEGs in second leaf and another 25, 252, and 16 DEGs in grain tissue, respectively. In the second leaf tissue, under N stress, the total up- and down-regulated DEGs at two different time points were variable across the cultivars. In Mace, at 0 DPA, the down-regulated and up-regulated DEGs were 434 and 102, respectively. At 10 DPA, the up-regulated and down-regulated DEGs were 109 and 74, respectively. Similarly, in Volcani, the down-regulated DEGs at 0 and 10 DPA were counted as 753 and 430, respectively, whereas the up-regulated DEGs were 354 at 0 DPA and 261 at 10 DPA. Cultivar Spitfire showed 536 up-regulated and 39 down-regulated DEGs at 0 DPA, whereas it showed 6624 up-regulated and 3830 down-regulated DEGs at 10 DPA. On the other hand, in grain tissue at 10 DPA, the down-regulated DEGs were 5, 237, and 8, while the up-regulated DEGs were 0, 15, and 8 identified in Mace, Spitfire, and Volcani, respectively. Variation in up- and down-regulated genes across the cultivars can be related to the difference in their response to N stress ([115]Figure 2). FIGURE 2. [116]FIGURE 2 [117]Open in a new tab Number of up- and down-regulated expressed at 0 and 10 DPA DEGs in wheat cultivars (A) Mace, (B) Spitfire, and (C) Volcani. Differentially expressed genes (DEGs), 0 days post anthesis (0 DPA), 10 days post anthesis (10 DPA). Common DEGs Between Leaf and Grain A total of 50 common DEGs were identified between the second leaf and grain tissue under the N stressed condition, of which 30 were down-regulated and 7 were up-regulated in both tissues. Thirteen DEGs showed inconsistent up- and down-regulation ([118]Supplementary Table 3). Several stress-related genes have been identified among those common DEGs with >log2 fold change, including plasma membrane ATPase, Serine protease HtrA-like, transcription factor AS2/LOB, etc. Several transmembrane transport-related proteins including sulfate transporter, glycosyltransferase, and WAT1-related protein were also common in second leaf and grain tissues. On the other hand, NUE-related glutamine synthetase and glutamine dumper were significantly up-regulated in second leaf tissues but down-regulated in the grain tissue of Volcani, indicating their tissue-specific expression. Another gene related to amino acid metabolism, isoaspartyl peptidase/L-asparaginase, was up-regulated in the second leaf and grain tissue of Spitfire and Mace, indicating non-specific tissue expression. In general, the common down-regulating DEGs were largely involved in carbohydrate metabolic process (chitinase, trehalose-6-phosphate synthase) and oxidation–reduction process (aldehyde dehydrogenase, peroxidase, methyl sterol monooxygenase 1-2, gibberellin 20 oxidase 2, catalase). The up-regulating common DEGs are involved in N compound metabolic process (glutamine synthetase), sulfate transmembrane transport (sulfate transporter), and amino acid metabolism (aminotransferase like protein, isoaspartyl peptidase/L-asparaginase). Common DEGs Between 0 and 10 DPA Under N stress, some DEGs showed consistent up- or down-regulation at both 0 and 10 DPA ([119]Figure 3) despite the fact that they were variable between the cultivars. For example, in Mace, a total of 28 DEGs were commonly expressed at both 0 and 10 DPA, and of them, 23 were down-regulated and 5 were upregulated under N stress. Some of those DEGs showed high fold change (> + 2 or < −2) including plant protein DUF1589 of uncharacterized protein function, gibberellin receptor GID1A, catalase, a two-component response regulator, and cytochrome P450 was down-regulated, whereas RADIALIS-like TF, glycosyltransferase, and receptor-like protein kinase were up-regulated. In contrast, in Spitfire, among the 310 commonly expressed DEGs at 0 and 10 DPA, 16 showed down-regulation and 294 showed up-regulation under N stress. Among the DEGs in Spitfire, the Dof zinc finger protein, two-component response regulator, glycine-rich protein-A3, and calcium-dependent protein kinase 15 were down-regulated (log2 fold change < −4.0), and the cinnamoyl CoA reductase, receptor-like kinase, protein kinase-like, translation initiation factor IF-2, aspartate-tRNA ligase, and a beta-glucosidase were up-regulated under N stress. In Volcani, a total of 127 DEGs were found to be expressed both at 0 and 10 DPA under N stress. Of them, 86 were down-regulated and 38 were up-regulated commonly at both DPA, whereas three DEGs were down-regulated at 0 DPA but up-regulated at 10 DPA. The top down-regulated DEGs included a chlorophyll a-b binding protein, methyltransferase, endo-1,3 beta-glucanase, and plant protein DUF1589 of uncharacterized protein function, whereas the top up-regulated DEGs included cinnamoyl CoA reductase, MYB TF, glycosyltransferase, and beta-glucosidase ([120]Supplementary Table 4). FIGURE 3. [121]FIGURE 3 [122]Open in a new tab Consistently expressed DEGs at two time points in the leaf tissue of wheat cultivars (A) Spitfire, (B) Mace, and (C) Volcani. Differentially expressed genes (DEGs), 0 days post anthesis (0 DPA), 10 days post anthesis (10 DPA). Common DEGs Among Cultivars Venn diagram analysis was used to identify the number of common DEGs among the cultivars ([123]Hulsen et al., 2008). In the second leaf tissue, down-regulated 4 DEGs at 0 and 10 DEGs at 10 DPA whereas only 3 up-regulated DEGs at 10 DPA were found in common. The common down-regulated DEGs were identified as glycine-rich protein A3, calcium-dependent protein kinase 15, etc. The common up-regulated DEGs were identified as sulfate transporter and L-allo-threonine aldolase, which is related to amino acid metabolism. However, in grain tissue, only two down-regulated DEGs identified across the three cultivars were annotated as LOB-domain containing proteins ([124]Table 1). TABLE 1. List of DEGs common among wheat cultivars: Mace, Spitfire, and Volcani. Tissue Up-/down-regulation Stage Gene_id Annotation Log2 fold change __________________________________________________________________ Volcani Spitfire Mace Leaf Down 0 DPA TraesCS3A02G439500LC Glycine-rich protein A3 −5.816 −4.5188 −5.5826 TraesCS3D02G150400LC Glycine-rich protein A3 −5.4414 −3.2827 −4.9171 TraesCS4A02G245300 Protein DETOXIFICATION −6.9194 −2.5502 −4.1621 TraesCS3D02G150300LC Calcium-dependent protein kinase 15 −4.3632 −4.3883 −5.9832 Down 10 DPA TraesCS2D02G555300 ARM repeat superfamily protein −2.7152 −3.2171 −3.3102 TraesCS2D02G259200 Two-component response regulator −1.9826 −2.7094 −2.2783 TraesCS1B02G388700 Methyltransferase −6.0313 −7.4915 −4.9876 TraesCS6B02G051800 Glycerol-3-phosphate acyltransferase −2.6867 −2.1339 −2.0549 TraesCS7D02G516800 Chaperone protein dnaJ −1.8551 −3.0128 −2.0382 TraesCS3D02G144900 Protein DJ-1 −3.295 −5.0529 −3.237 TraesCS5A02G472500 Amino acid transporter, putative −4.416 −2.385 −2.4211 TraesCS2B02G277300 Two-component response regulator −2.7784 −3.5708 −2.8394 TraesCS3D02G316900LC Nucleoside triphosphatase I −4.4166 HN −4.1874 TraesCS7D02G388400 Tryptophan synthase beta chain −2.9456 −1.4694 −1.9753 Up 10 DPA TraesCS7D02G084100 Sulfate transporter 2.0942 3.4823 2.2596 TraesCS7B02G128800 Epoxide hydrolase 2 4.7394 1.3931 2.5759 TraesCS2D02G379000 L-allo-threonine aldolase 3.7648 2.4275 2.1093 Grain Down 10 DPA TraesCS2A02G194500 LOB domain-containing protein, putative −1.3364 −2.2599 −1.6826 TraesCS2D02G193400 LOB domain-containing protein −1.5747 −2.1909 −1.8922 [125]Open in a new tab Differentially expressed genes (DEGs), 0 days post anthesis (0 DPA), 10 days post anthesis (10 DPA). While considering the common DEGs between two cultivars, in most cases, the highest number of DEGs was common between Spitfire and Volcani among all combinations ([126]Figure 4). The major common down-regulated DEGs between Spitfire and Volcani in second leaf were identified as methyltransferase, chlorophyll a-b binding protein, methyltransferase, and aquaporin, whereas the common up-regulated DEGs were identified as aminotransferase, early light-induced protein, F-box domain-containing protein, and glycosyltransferase. In grain tissues, among the four common down-regulated DEGs between Spitfire and Volcani, cysteine proteinase inhibitor and Ureide permease-like protein are related to N metabolism. The top down-regulated DEGs common between Mace and Volcani were identified as photosystem II 10 kDa polypeptide family proteins, chlorophyll a-b binding protein, and plant protein DUF1589 of uncharacterized protein function, whereas plasma membrane ATPase and glycosyltransferase were found as the common top up-regulated DEGs. Similarly, the common top up-regulated DEGs between Mace and Spitfire were identified as a vacuolar-processing enzyme, a boron transporter, a nuclease S1, and cytochrome P450 family protein, whereas down-regulated DEGs were identified as haloacid dehalogenase-like hydrolase (HAD) superfamily protein and thaumatin-like protein. FIGURE 4. [127]FIGURE 4 [128]Open in a new tab Venn diagrams of DEGs shared among Mace, Spitfire, and Volcani in leaf and grain tissues at two development stages. (A) The number of down-regulated genes in leaf. (B) The number of down-regulated genes in grain. (C) The number of up-regulated genes in leaf. (D) The number of up-regulated genes in grain. Differentially expressed genes (DEGs), 0 days post anthesis (0 DPA), 10 days post anthesis (10 DPA). Gene Ontology Reflects the Function of DEGs in Response to Nitrogen Stress The top 10 biological process GO terms characteristic to the DEGs are presented in [129]Figure 5. The frequency of the GO term is shown as percentage of the genes compared to the total gene number related to the GO term. FIGURE 5. [130]FIGURE 5 [131]Open in a new tab Top biological process GO terms in leaf and grain tissue of wheat cultivars (A) Mace, (B) Spitfire, and (C) Volcani. The frequency of the GO term is shown as percentage of the genes related to the GO term. GO, gene ontology. In the case of the second leaf tissue, transmembrane transport GO term appeared as the top group within up-regulated DEGs in all three cultivars. Notably, another three top GO terms were common in Spitfire and Volcani, which were ion transport, lipid metabolic process, and metal ion transport, indicating that these two cultivars have some common physiological response mechanisms to N stress. In contrast, Mace did not have any other top 10 GO common with either cultivar. DNA metabolic process and organelle organization are the next top GO terms for cultivar Mace. On the other hand, genes showing decreased expression under N stress in Mace second leaf tissue were mostly involved in defense response and carbohydrate metabolic process. In Spitfire, decreasing gene expression was largely related to photosynthesis and light harvesting, organonitrogen compound biosynthesis process, and small molecule biosynthetic process. Similarly, in Volcani, genes with decreased expression patterns were also related to photosynthesis, carbohydrate metabolic process, and response to external stimulus. In the grain tissue, the transmembrane transport process GO term was the top enriched group among the up-regulated DEGs in Mace and Spitfire. Mace also showed enrichment in carbohydrate metabolic process and ion transport. However, proteolysis and negative regulation of catalytic activity were common in Spitfire and Volcani among the top 10 enriched GO terms. Nitrogen compound transport appeared as the common GO term in all three cultivars among the down-regulated DEGs. Mace did not show enrichment for negative regulation of endopeptidase activity and proteolysis like Spitfire and Volcani. However, in Mace and Spitfire, glycolysis process was the top enriched down-regulated GO term. KEGG Analysis Spanned Function of DEGs in Response to Nitrogen Stress Using the well-annotated rice genome as a reference, KEGG pathway enrichment analysis identified significantly enriched metabolic pathways and signal transduction pathways associated with DEGs. The top 20 most significantly enriched pathways were selected to produce the KEGG scatter plot ([132]Supplementary Figure 1). Results for the KEGG pathway terms that were significant at adjusted p-value q-≤0.5 are shown in [133]Tables 2, [134]3 for second leaf and grain, respectively. Under N stress, a total of 41 KEGG pathway terms were significantly associated with 12,108 DEGs in the second leaf, and 14 KEGG pathway terms were associated with 276 DEGs in the grain tissue. Among the 41 significant KEGG terms for the second leaf tissue, 3, 15, and 6 KEGG terms were specific to Mace, Spitfire, and Volcani, respectively, whereas one KEGG term was common between Mace and Spitfire, six KEGG terms between Mace and Volcani, and five KEGG terms between Spitfire and Volcani ([135]Table 2). There were five KEGG terms common in all three cultivars under N stress, namely, phenylpropanoid biosynthesis, biosynthesis of secondary metabolites, flavonoid biosynthesis, metabolic pathways, and starch and sucrose metabolism. Among the 14 significant KEGG terms associated with DEGs in grain, eight, four, and one KEGG terms were specific to Mace, Spitfire, and Volcani, respectively ([136]Table 3). The DEGs in the grain of all three cultivars were commonly associated with the KEGG pathway term glycolysis/gluconeogenesis. Among the significant KEGG terms for DEGs in the second leaf, zeatin biosynthesis, arginine and proline metabolism, and sulfur metabolism were specific to Mace, with terms like plant–pathogen interaction, photosynthesis, pentose phosphate pathway, porphyrin, and chlorophyll metabolism specific to Spitfire and beta-alanine metabolism, tryptophan metabolism, ubiquinone, and other terpenoid-quinone biosynthesis pathways found only in DEGs of Volcani. In the grain tissue, the KEGG pathways specific to Mace were glycerolipid metabolism, sphingolipid metabolism, porphyrin and chlorophyll metabolism, and galactose metabolism, whereas the pathways specific to Spitfire were alanine, aspartate and glutamate metabolism, glycine, serine, and threonine metabolism, and ribosome biogenesis in eukaryotes. The pathways specific to Volcani were cysteine and methionine metabolism. In addition, some KEGG pathways were common between two cultivars only, e.g., Mace and Spitfire had a MAPK signaling pathway common in the second leaf and glycolysis/gluconeogenesis common in the grain tissue. The KEGG terms common in Mace and Volcani included glutathione metabolism, galactose metabolism, and ABC transporters, whereas biosynthesis of amino acids, photosynthesis–antenna proteins, and circadian rhythm–plant pathways were common between Spitfire and Volcani. TABLE 2. Significantly (adjusted p-value ≤ 0.5) enriched KEGG pathways in Mace, Spitfire, and Volcani under nitrogen stress in second leaf tissue. Cultivar KEGG pathway KEGG id Mace Zeatin biosynthesis osa00908 Arginine and proline metabolism osa00330 Sulfur metabolism osa00920 Spitfire Plant–pathogen interaction osa04626 Carbon metabolism osa01200 Glyoxylate and dicarboxylate metabolism osa00630 Photosynthesis osa00195 Glycine, serine, and threonine metabolism osa00260 Carbon fixation in photosynthetic organisms osa00710 Glycolysis/gluconeogenesis osa00010 Fructose and mannose metabolism osa00051 Alanine, aspartate, and glutamate metabolism osa00250 Taurine and hypotaurine metabolism osa00430 Pentose phosphate pathway osa00030 One carbon pool by folate osa00670 Porphyrin and chlorophyll metabolism osa00860 Histidine metabolism osa00340 Ascorbate and aldarate metabolism osa00053 Volcani Beta-alanine metabolism osa00410 Terpenoid backbone biosynthesis osa00900 Ubiquinone and other terpenoid-quinone biosynthesis osa00130 Tryptophan metabolism osa00380 Butanoate metabolism osa00400 Phenylalanine, tyrosine, and tryptophan biosynthesis osa00651 Mace–Spitfire MAPK signaling pathway—plant osa4016 Mace–Volcani Stilbenoid, diarylheptanoid, and gingerol biosynthesis osa00945 Plant hormone signal transduction osa04075 Phenylalanine metabolism osa00360 Galactose metabolism osa00052 ABC transporters osa02010 Glutathione metabolism osa00480 Spitfire–Volcani Biosynthesis of amino acids osa01230 Photosynthesis-antenna proteins osa00196 Circadian rhythm—plant osa04712 Cyanoamino acid metabolism osa00460 Cysteine and methionine metabolism osa00270 Mace–Spitfire–Volcani Phenylpropanoid biosynthesis osa00940 Biosynthesis of secondary metabolites osa01110 Flavonoid biosynthesis osa00941 Metabolic pathways osa01100 Starch and sucrose metabolism osa00500 [137]Open in a new tab TABLE 3. Significantly (adjusted p-value ≤ 0.5) enriched KEGG pathway in Mace, Spitfire, and Volcani under nitrogen stress in grain tissue. Cultivar KEGG pathway KEGG id Mace Galactose metabolism osa00052 Oxidative phosphorylation osa00190 Glycerolipid metabolism osa00561 Sphingolipid metabolism osa00600 Glycosphingolipid biosynthesis—globo series osa00603 Porphyrin and chlorophyll metabolism osa00860 Biosynthesis of secondary metabolites osa01110 RNA degradation osa03018 Spitfire Alanine, aspartate, and glutamate metabolism osa00250 Alanine, aspartate, and glutamate metabolism osa00250 Glycine, serine, and threonine metabolism osa00260 Ribosome biogenesis in eukaryotes osa03008 Volcani Cysteine and methionine metabolism osa00270 Mace–Spitfire Glycolysis/gluconeogenesis osa00010 [138]Open in a new tab Protein–Protein Interaction Network Analysis of DEGs MCL clustering using the inflation parameter 1.70 was used to show the association of clusters in KEGG pathways ([139]Figures 6, [140]7). Networks showed that a large number of proteins were involved in multiple interactions and grouped into seven major clusters. Among the seven clusters, two large clusters were enriched in photosynthesis and steroid biosynthesis. All the interacting DEGs identified as photosynthesis-related, and photosynthesis antenna proteins were down-regulated while some DEGs related to steroid biosynthesis were up- or down-regulated. Among the other clusters, the majority of down-regulated DEGs were involved in carbohydrate metabolism, amino sugar, and nucleotide metabolism, whereas up-regulated DEGs were mostly related to amino acid metabolism and signaling. In the biosynthesis of secondary metabolites, both the up- and down-regulated DEGs were involved. Overall, the number of down-regulated DEGs was higher in the network and was mainly involved in photosynthesis and photosynthesis-antenna proteins. FIGURE 6. [141]FIGURE 6 [142]Open in a new tab Protein–protein interaction network analysis of DEGs under N stress using Oryza sativa as reference. The different highlighted color indicates the different clusters of DEGs involved in different KEGG pathways. Differentially expressed genes (DEGs), Kyoto Encyclopedia of Genes and Genomes (KEGG). FIGURE 7. [143]FIGURE 7 [144]Open in a new tab Protein–protein interaction network analysis of DEGs under N stress using Arabidopsis thaliana as reference. The different highlighted color indicates the different clusters of DEGs involved in different KEGG pathways. Differentially expressed genes (DEGs), Kyoto Encyclopedia of Genes and Genomes (KEGG). Identification of Nitrogen Metabolism-Related Genes in Response to Nitrogen Stress N metabolism is a vital biological process in plants that determines crop productivity and yield ([145]Stitt et al., 2002; [146]Cai et al., 2009). The DEGs involved in N uptake, transport, and assimilation were listed separately and are presented in [147]Table 4. Most of the N metabolism-related DEGs showed up-regulation under N stress. Among the most significant DEGs (fold change >2.0), 65% were up-regulated and 35% were down-regulated ([148]Table 5). Spitfire showed abundancy for N metabolism-related DEGs compared to Mace and Volcani. The top up-regulated N metabolism-related DEGs included amino acid permease, glutamate dehydrogenase, low-affinity nitrate transporter protein NRT1/PTR family 1.1, tyrosine aminotransferase, and high-affinity nitrate transporter, whereas the top down-regulated DEGs included amino acid transporter family protein, nitrate transporter 1.1, nitrate transporter 1.2, nitrate reductase, and tryptophan aminotransferase. Spitfire showed the most induction ratio for protein NRT1/PTR FAMILY 1.1 (log2 fold change 6.4) and tyrosine aminotransferase (log2 fold change 5.74). Mace showed up-regulation of cationic amino acid transporter and down-regulation of amino acid transporter family proteins, amino acid permease, and protein NRT1/PTR FAMILY 1.1. Volcani showed up-regulation of amino acid permease, nitrate transporter protein NRT1/PTR FAMILY 5.5, and ammonium transporter and down-regulation of isoaspartyl peptidase/L-asparaginase, nitrate transporter 1.1 and 1.2, and tryptophan aminotransferase. TABLE 4. Differentially expressed genes (DEGs) involved in nitrogen uptake, transport, and assimilation. Name NUE class DEG count (current study) NUE related effect Host References