RNA-Seq-Based Transcriptomics Study to Investigate the Genes Governing Nitrogen Use Efficiency in Indian Wheat Cultivars

Abstract

High NUE (nitrogen use efficiency) has great practical significance for sustainable crop production. Wheat is one of the main cultivated crops worldwide for human food and nutrition. However, wheat grain productivity is dependent upon cultivars with high NUE in addition to the application of nitrogen fertilizers. In order to understand the molecular mechanisms exhibiting a high NUE response, a comparative transcriptomics study was carried out through RNA-seq analysis to investigate the gene expression that regulates NUE, in root and shoot tissue of N-efficient (PBW677) and N-inefficient (703) cultivars under optimum and nitrogen (N) stress. Differentially expressed gene analysis revealed a total of 2,406 differentially expressed genes (DEGs) present in both the contrasting cultivars under N stress. The efficient genotype PBW677 had considerably more abundant DEGs with 1,653 (903 roots +750 shoots) compared to inefficient cultivar PBW703 with 753 (96 roots +657 shoots). Gene ontology enrichment and pathway analysis of these DEGs suggested that the two cultivars differed in terms of adaptive mechanism. Gene enrichment analysis revealed that among the upregulated and downregulated genes the overrepresented and underrepresented gene categories belonged to biological processes like DNA binding, response to abiotic stimulus, photosynthesis, carbon fixation, carbohydrate metabolic process, nitrogen compound metabolic process, nitrate transport, and translation in cultivar PBW677, while the enriched biological processes were nucleosome assembly, chromatin remodeling, DNA packaging, lipid transport, sulfur compound metabolic process, protein modifications, and protein folding and refolding in N inefficient cultivar PBW703. We found several transcription factors (MYB, WRKY, RING finger protein, zinc finger protein, transporters, NRT1, amino acid transporters, sugar), protein kinases, and genes involved in N absorption, transportation, and assimilation to be highly expressed in high NUE cultivar PBW677. In our study, we report 13 potential candidate genes which showed alternate gene expression in the two contrasting cultivars under study. These genes could serve as potential targets for future breeding programs.

Introduction

Nitrogen is one of the essential macronutrients required for plant growth, development, and reproduction. In last few decades, extensive use of N fertilizers resulted in increased biomass and crop yield. It has been estimated that by the year 2050, the application of N fertilizer will increase to 135 Tg N (Good et al., 2004). But excessive use of N causes environmental pollution, deteriorates soil health, and also leads to a higher cost of production. Thus, increased N use efficiency in plants would not only result in higher crop yield under N-stressed conditions, but also benefit farmers via higher net profit under low input and mitigate the environment risks arising due to an excess of fertilizers. Therefore, one of the main aims in agricultural research is to increase the NUE of plants which depends on plant N-uptake efficiency, N-utilization efficiency, and on the remobilization efficiency of nitrogen from dead tissue to growing plant parts (Masclaux-Daubresse et al., 2010).

Grain production in major cereals largely depends upon the application of N fertilizers and cultivars with high NUE (Hitz et al., 2017). Grain production in crop plants could be enhanced by harnessing the genetic variation for improved NUE. Nitrogen use efficiency is strongly affected by genetic as well as environmental factors (Xu et al., 2012). The key genes involved in the pathways involved in NUE can be discovered using the power of omics. Transcriptomics or RNA-seq is one such approach which can be utilized to reveal the key genes responsible for combating N stress under N-deficit conditions using contrasting genotypes or cultivars having different NUE (Kant et al., 2011). RNA-seq via next-generation sequencing platform allows the discovery of genes playing a role in pathways affecting such phenotyping traits. Many studies are available illustrating the potential of transcriptomics to decipher the role of key genes in nitrogen-dependent pathways, for example the genotypic difference in terms of nitrogen use efficiency has been studied in various crops like rice (Vijayalakshmi et al., 2015), soybean (Hao et al., 2011), sorghum (Gelli et al., 2014), and maize (Zamboni et al., 2014).

Wheat is one of the main cultivated crops worldwide (Xu et al., 2013). To meet the needs of the world’s growing population, the grain yield of wheat must be increased at an average annual rate. However, wheat yield is frequently threatened by low-nitrogen stress. Therefore, improving the nitrogen use efficiency for increased biomass as well as increased yield of wheat under N-starvation conditions has great practical significance (Curci et al., 2017). Hence, efforts have been made to understand the molecular and physiological basis of plants grown under N stress, which will help in identification of a large number of differentially expressed genes (DEGs) through RNA sequence-based transcriptomics analysis (Ruuska et al., 2008).

In wheat, attempts have been made to understand the regulatory mechanisms involved in nitrogen metabolism and various studies have been conducted so far to discover the genes via transcriptome sequencing. The RNA-seq studies have been carried out taking various tissue samples, for example leaf tissue, roots, flag leaf, etc., from a single cultivar and differential gene expression was carried out to unravel the key genes and pathways (Meng et al., 2021; Sultana et al., 2021; Tiong et al., 2021).

We have conducted a comparative transcriptomics study using the next-generation sequencing platform Illumina HiSeq 500 using two bread wheat cultivars: PBW677 with high NUE and PBW703 with low NUE (Shamshad, 2019) and studied the whole genome gene expression profile in control as well nitrogen-deficient conditions in root and shoot tissue. This is the first report to be carried out using the two contrasting cultivars and finding the gene expression in root and shoot tissue in normal vs. stress (N-deficient) conditions. As per our knowledge only one study has been conducted in bread wheat taking flag leaf and second leaf tissue after 0 DPA and 10 DPA from three different cultivars having high, medium, and low NUE (Sultana et al., 2020).

Materials and Methods

Plant Materials, Nitrogen Treatments, and Plant Sampling

Plant material included two wheat (Triticum aestivum) cultivars, PBW677 (PFAU/MILAN/5/CHEN/Ae. squarrosa//BCN/3/VEE#7/BOW/4/PASTOR) known to show efficiency at a low input of nitrogen (called N efficient) and PBW703 (called N inefficient) (BWL9250*3//Yr10/6*Avocet/3/BW9250*3//Yr15/6*Avocet) with low NUE (Shamshad, 2019). The field experiments were carried out using a split-plot design with two N treatments (control and no nitrogen) as the main plot and the two cultivars in three replications as sub plots (2.4 m × 5 m). With each subplot, seeds were sown in 12 5 cm-long rows with 20 cm spacing. The nitrogen was supplied at the rate of 120 Kg ha-1 as controlled conditions and no nitrogen was supplied to induce stress conditions. The crop was grown by following the standard agronomic packaging and practices. For RNA-seq, seeds of both the varieties were surface-sterilized and grown as two sets of three replications in a growth chamber maintained at 26/24°C with relative humidity of 75%. One set of genotypes was sown in perchlorate-contaminated soil without adding any nitrogen fertilizer as N-stress (N⁻) while the other set was given nitrogen at a recommended rate of 120 kg/acre as N-control (N⁺). To the N-control set, nitrogen fertilizer was applied in the form of urea 21 days after the germination of seeds. Twenty-4 hours after application of urea, root and shoot samples were collected both from N-stress and N-control sets at the same time. Each biological replicate constituted a pool of three plants and a total of three biological replicates were used. The collected eight tissues, (N-stress root, N-control root, N-stress shoot, and N-control shoot) from both genotypes were immediately frozen in liquid nitrogen and stored at −80°C.

Analysis of Growth, Yield, and N Parameters

Three plots per genotype and N treatment were available for the data collection. Leaf chlorophyll content was measured by randomly selecting five plants from each row using an SPAD meter, plant height was measured using scale in cm, and biological yield and grain yield were recorded as yield per plot in Kg. Grain nitrogen was estimated by using the method proposed by Eastin (2008) from the grounded grain, straw, and leaf tissue samples. The total N content was measured using the distillation and titration method (Pelican Equipment, CIT Nagar, Chennai, Tamil Nadu, India) following the manufacturer’s protocol. NUE, N utilization efficiency, and N harvest index were calculated according to Moll et al. (1982). Seed yield relative to total N accumulation in above-ground tillers was used to calculate NUtE, and NUE was determined as the NUpE and NUtE. NHI was analyzed using the ratio of total N in grains to tillers and grains.

RNA Extraction and Reverse Transcription

Total RNA was extracted from 100 mg of root and shoot samples using Trizol Reagent (Takara kit) as per the manufacturer’s protocol. The concentration and quality of total RNA were monitored on 1.2% agarose gels as well as on a NanoDrop™ 8,000 Spectrophotometer (Thermo Scientific). Approximately 1 µg of RNA was used as a template for first-strand cDNA synthesis (Thermo Scientific Verso cDNA Kit). Oligo (dT)₁₈ primers were used during reverse transcription. The cDNA was stored at −20°C.

RNA-Seq Analysis and High-Quality Read Statistics

To obtain a global overview of the wheat transcriptome and gene activity at a nucleotide resolution, cDNA samples were sequenced by the Illumina Nextseq500 platform. The generated raw reads were submitted to the NCBI sequence read archives (SRA) with accession number PRJNA780342. The raw sequences were assessed for quality using FastQC version 0.11.2 (Andrews et al., 2015). Adaptors clipping and quality trimming of raw reads were performed using Trimmomatic v0.39 software (Bolger et al., 2014). Low quality reads with phred score <30 and read length <50 bp were removed.

Read Alignment and Assessment to Mapping With Reference Wheat Genome

The high-quality transcriptome reads were aligned to the indexed wheat reference genome (http://www.ensembl.org/index.html) using Burrows-Wheeler Aligner (BWA) software (BWA-0.7.17) (Li and Durbin, 2009). The reference Triticum aestivum genome (IWGSC RefSeq v1.0) and the associated annotations (IWGSC RefSeq v1.1) were downloaded from the International Wheat Genome Sequencing Consortium (IWGSC). Indexing of the reference genome of wheat (Triticum aestivum) was done using BWA-0.7.17 (Li and Durbin, 2009). The resulted output sequence alignment map (SAM) files were converted into binary alignment map (BAM) format, sorted, and indexed using samtools v1.10 (Li et al., 2009).

Differential Gene Expression Analysis

The statistical model of the Cufflinks-cuffdiff v2.2.1 package (Trapnell et al., 2012) was used to assemble and quantify differential gene expression in terms of FPKM (fragments per kilobase per million reads). Genes with log2 fold change (FC) values > 2 were considered upregulated whereas FC < −2 were considered downregulated. These genes were further categorized on the basis of statistical significance (p < 0.05).

Functional Annotation and Gene Ontology Term Analysis

Functional Gene Ontology was performed with OmicsBox_windows-64_1_3_11 (http://wego.genomics.org.cn/cgi-bin/wego/index.pl). The FASTA sequences of all the upregulated and downregulated genes were downloaded from the Ensemble Plants database (https://plants.ensembl.org/index.html) and used as input into OmicsBox which contains information of all the genes assigned into three main GO domains, viz., biological process, cellular component, and molecular function. All DEGs were annotated according to wheat IWGSC release 1.1, whereas the sequences lacking annotation in the wheat reference were annotated with Arabidopsis thaliana (Supplementary Tables S1, S2). Volcano plots were made using R 4.0.4 (ggplot2 package) software which relates the observed differences in gene expression associated with Cuffdiff’s statistical model. Venn diagrams were constructed to represent up and downregulated DEGs using the Venny 2.1 tool (https://bioinfogp.cnb.csic.es/tools/venny/) to analyze tissue-specific and tissue-independent genes. R-based 4.0.4 software was used to make heat maps which present hierarchical clustering based on log2 fold changes to visualize the expression patterns of DEGs.

Enrichment Analysis Based on Gene Ontology Terms

To see which class or category of genes were overrepresented and underrepresented among the differentially expressed genes in response to N stress, enrichment analysis was carried out using Fisher’s exact test with FDR <0.05. Enriched bar graphs were made using OmicsBox version 2.0.36.

KEGG Pathway Analysis

In order to reveal the pathways to which DEGs under N-stress belong, KEGG pathway analysis was carried out by using the FASTA sequences of all the differentially expressed genes as input to the OmicsBox 2.0.36 combined pathway analysis plugin.

Real-Time Quantitative Polymerase Chain Reaction Analysis

Real-time PCR was performed using the LightCycler 96 Roche Real-time PCR system and PowerUp™ SYBR™ Green Master Mix (applied biosystem by Thermo Fisher Scientific). Seventeen differentially expressed genes (from shoots and roots of both low N and high N) were selected for validation. Primer3 version 2.4.0 was used to design gene-specific primers and their specificity was verified using the NCBI database through the Blast tool (Supplementary Table S5). The 10 µl RT-qPCR reaction contained 1 µl of template cDNA (20 ng), 1 µl of forward primer, 1 µl of reverse primer, 4 µl of PowerUp™ SYBR™ Green Master Mix, and 3 µl of H₂O. PCR was run at an initial denaturation of 94°C for 3 min followed by 40 cycles of 94°C for 10 s, 60°C for 30 s, 72°C for 30 s, and a final extension at 72°C for 10 min to check the specificity of amplification. The housekeeping gene TaATP (ATP-dependent 26S proteasome regulatory) (Paolacci et al., 2009) was used as the endogenous control and all reactions were performed in triplicate. Relative gene expression was analyzed using the 2−^ΔΔCt method.

Results

Plant Growth, Total Chlorophyll, and Total N Content Analysis

Both the cultivars PBW677 and PBW703 showed differences in growth and yield characteristics under contrasting N conditions (Figure 1A). Plant height showed no differences between the cultivars (based on cultivar) however it showed a significant difference at contrasting N levels (Figure 1B). The chlorophyll content showed a significant difference among the cultivars at different N levels (Figures 1C,D). Further no differences in the spikelets per spike were found under various N levels (Figure 1F). Grain and straw N were lower in PBW703 under N-stressed conditions in comparison to PBW677. The 1,000 grain weight was significantly higher for PBW677 under all the N fertilization regimes (Figure 1J) However, no significant difference was observed in spikelets per spike at different N levels. The nitrogen utilization and harvest index were higher in PBW677 vs. PBW703 supporting the fact that PBW677 is more efficient in acquiring, locating, and using N for grain development. Indeed, total tiller N levels as well as the NHI were significant higher in PBW677 compared to PBW703 (Figures 1M,N).

FIGURE 1

Transcriptomics Quality and Mapping Statistics

A total of 132 Gb 150-bp paired end (PE) reads were generated through Illumina NextSeq⁵⁰⁰ from eight samples (2 genotype × 2 treatment × 2 tissue) including N⁻ and N⁺ treatment. On average, 18 million clean reads were obtained after trimming/clipping from each library with sizes ranging from 10 to 17 GB. The reads were of high quality at 98 and 95% and had a phred score of Q33. Moreover, the average GC% of each library was about 56% (Table 1). After mapping against the wheat reference genome (IWGSC release 1.0), varying contents of perfectly mapped reads in both roots (68.35%) and shoots/leaves (86.58%) were recorded in both genotypes (Table 2). Transcript profiles of the RNA-seq data were analyzed by calculating the read fragments per kilobase per million mapped reads (FPKM).

Differential Gene Expression in Response to Nitrogen Stress

Transcriptomics analysis of N-efficient (PBW677) and N-inefficient (PBW703) genotypes revealed a significant difference of adaptive response in terms of gene expression pattern when there was nitrogen stress. In PBW677, the total number of genes expressed in root tissues was 70,825, out of which 903 genes showed significant differential expression (748 genes upregulated and 155 downregulated) while in shoot tissues, a total of 61,895 genes were expressed, out of which 750 were significant differentially expressed genes (667 upregulated, 83 downregulated) (Figure 2; Table 3). In contrast, the total number of expressed genes in root tissues of PBW703 was 22,046, of which only 96 were differentially expressed genes (92 upregulated and four downregulated) and in shoot tissues, 49,121 genes were expressed, out of which 657 were differentially expressed genes (511 upregulated and 146 downregulated), as shown in Figure 2. Supplementary Tables S1, S2 present the list of DEGs in root and shoot tissues of the cultivars PBW677 and PBW703, respectively. Volcano plots in Figure 3 present the most significant upregulated and downregulated genes in root and shoot tissues of both the cultivars.

FIGURE 2

FIGURE 3

Genes Involved in Primary Nitrogen Metabolism

Genes involved in N absorption and assimilation were found to be differentially expressed in N stress-tolerant cultivar PBW677 compared to N stress-intolerant cultivar PBW703. Genes corresponding to the NRT1/PTR gene family were found to be upregulated in root tissues of PBW677 and downregulated in shoot tissues of PBW677 and PBW703. Glutamate synthase (GLT), amino acid transporter (AVT), 1-aminocyclopropane-1-carboxylate synthase (ACO), IAA-amino acid hydrolase (ILL6), amino acid permease (AAP), and Asparagine synthase (ASNS) were upregulated in both the root and shoot tissues of PBW677 in comparison to PBW703 where these genes showed downregulation. However, some DEGs related to serine--glyoxylate aminotransferase (AGT), glyoxylate aminotransferase (GGAT1), and chloroplast/mitochondrial Glutamine synthetase (GLN) were upregulated in roots of PBW703 compared to PBW677. In addition, Glutamate dehydrogenase (GSH)-encoding genes were upregulated in shoot tissues of both the cultivars (Figure 4A, Supplementary Table S3).

FIGURE 4

Carbon Metabolism Genes

Alternate gene expression was observed for genes involved in carbon metabolism, especially involved in pathways like glycolysis, tricarboxylic acid (TCA) cycle, gluconeogenesis, photorespiration, photosynthesis, and trehalose metabolism as represented in Figure 4B. It has been observed that several DEGs related to glycolysis such as glyceraldehyde-3-phosphate dehydrogenase GAPA1, fructose-bisphosphate aldolase, and fructose-1,6-bisphosphatase were downregulated in both the tissues of PBW677. In contrast, in PBW703, these genes showed upregulation in root tissues. The same pattern was observed for gluconeogenesis metabolism-related genes, such as fructose-bisphosphate aldolase, sedoheptulose-1,7-bisphosphatase, and UDP-glucose 6-dehydrogenase that were found to be downregulated in roots of PBW677 but upregulated in PBW703. Similarly, the majority of DEGs involved in photosynthesis like ribulose bisphosphate carboxylase, chlorophyll a-b binding protein of LHCII, oxygen-evolving enhancer protein 2, photosystem II 5 kDa protein, photosystem I reaction center subunit III, cytochrome b6-f complex iron-sulfur subunit, and ferredoxin-NADP reductase were upregulated in both the tissues of PBW703 (Figure 4B, Supplementary Table S3). Other than this, the resulting DEGs that were involved in energy metabolic pathways other than carbon metabolism such as genes coding NADP-dependent malic enzyme and trehalose-phosphate phosphatase were found to be strongly upregulated under N stress in both tissues of PBW677 but were considerably downregulated in PBW703.

Transcription Factor (TF)-Encoding Genes

Transcription factors play a vital role as molecular switches controlling the expression of specific genes and play crucial roles in plant development, cell cycling, cell signaling, and stress responses. Due to N stress, various unique genes encoding transcription factors were differentially expressed in N stress-tolerant cultivar PBW677 corresponding to different TF families of MYB, WRKY, RING finger protein, Zinc finger protein, TIFY proteins, AP2/ERF, HSP, BTB/POZ domain, NDR1/HIN1domain, bZIP, GATA, bHLH, DIVARICATA factors, Bowman-Birk type trypsin inhibitor isoform, BAG family molecular chaperone, dnaJ protein, ARF, and other families (Supplementary Table S3). In contrast, TF-encoding genes were identified in N stress-sensitive cultivar PBW703 and the number of TF-related DEGs in each family was less than those found in PBW677. We also observed the tissue-specific expression of several TF families. For instance, most of WRKY, HSP and AP2/ERF genes were extensively upregulated in roots only, however genes related to the BAG family molecular chaperone, ARF, DIVARICATA, GATA, and dnaJ protein were upregulated majorly in shoots of PBW703. Additionally, MYB and bHLH genes were expressed in both the tissues of PBW677 (Supplementary Fig. S1, Supplementary Table S3).

Transporter Encoding Genes

Genes encoding for various transporters responsible for uptake and transport of nutrients were found to be differentially expressed under N stress in both the cultivars. In general, it was observed that the genes belonging to different transporter families were shown to be unregulated in both the tissues of PBW677 in comparison to PBW703 where they were downregulated. In particular, calcium-binding protein transporter genes were upregulated in both the tissues of PBW677 in addition to the ABC transporter and SWEET transporter (Supplementary Fig. S2, Supplementary Table S3).

Protein Kinases Encoding Genes

In this study, various unique protein kinases were differentially expressed in both the tissues. The identified PK genes belonged to classes CDK (cyclin-dependrnt kinase), CIPK (CBL-interacting protein kinase), cysteine-rich domains, F-box domains, LRR (leucine rich repeats), MAPK (mitogen activated protein kinase), serine threonine protein kinase, serine arginine protein kinase, and U-box domains (Supplementary Fig. S3, Supplementary Table S3). We found PKs belonging to MAPK, serine threonine protein kinase, CDK, CIPK, and cysteine-rich domains upregulated in roots of PBW677. Besides, most members of F-box domains and LRR were found upregulated only in shoots of PBW677.

Other Stress-Related Genes

In this study, we also observed upregulation of a number of genes related to detoxification and protection from oxidative stress, the majority of which were found in both tissues in N stress-tolerant cultivar PBW677. These mainly belong to classes of glutathione S-transferases (GSTs), cytochrome P450 (CYP450), E3 ubiquitin-protein ligase, peroxidase and phenylalanine ammonia-lyase, and aquaporin PIP1-1. Most of the genes related to glutathione S-transferases (GSTs), phenylalanine ammonia-lyase (PAL), and E3 ubiquitin-protein ligase were upregulated in both root and shoot tissues of PBW677 (Figure 5) (Supplementary Figure S4, Supplementary Table S3).

FIGURE 5

Common Genes Between Root and Shoot Tissues of PBW677 and PBW703

Among the total DEGs studied, common genes showing up and downregulation were detected, also the genes showing alternate expression (that is upregulation in one cultivar and downregulation in another cultivar) were identified. A list of common genes with their log fold change value and involved biological pathways is represented in Supplementary Table S4.

Enriched Gene Ontology Terms in Differentially Expressed Genes Under Nitrogen Stress

In cultivar PBW677, enrichment analysis using Fisher’s exact test depicted the over and underrepresented genes expressed in response to nitrogen stress, as shown in Figure 6. Among the upregulated DEGs in root tissues of cultivar PBW677, 140 genes belonging to different GO categories (biological process, cellular component, and molecular function) were found to be overrepresented. Most of the overrepresented class of genes was found to be involved in sequence-specific DNA binding, response to temperature stimulus, abiotic stimulus, and transcription factors involved in transcription regulation and defense responses. Among the underrepresented class, only 11 genes were found, most of which belonged to the process of translation. The top 10 biological processes that were overrepresented among the upregulated DEGs expressed in root tissues of PBW677 are represented in Table 4. Among the downregulated DEGs in root tissues of PBW677, there were 64 genes which were overrepresented and belonged to the biological processes photosynthesis, carbon fixation, and photorespiration, as represented in Table 5. However, in the shoot tissues of PBW677, 82 upregulated genes were found to be overrepresented and six were underrepresented. In contrast, among the upregulated DEGs in root tissues of PBW703, 163 genes were found to be overrepresented and among them the biological processes photosynthesis, phosphorespiration, carbon fixation, and carbon metabolism, glutamine biosynthetic process, and nitrogen compound metabolic processes were highly enriched and overrepresented, as represented in Figure 7. In the shoot tissues of PBW703, among the upregulated DEGs (Figure 7), the overrepresented biological processes included nucleosome assembly, chromatin assembly, DNA packaging, chromatin remodeling, chromosome organization, response to inorganic substance, abiotic stimulus, lipid transport, sulfur compound metabolic process, glutathione metabolic process, nitrate transport, and response to nitrate with a number of 126 genes. Table 6 represents the top 10 highly enriched biological processes overrepresented in upregulated DEGs in shoot tissues of PBW703. The underrepresented category included macromolecule modification, protein metabolic process, protein modification process, and nitrogen compound metabolic process. Protein folding and refolding and responses to heat and temperature stimulus were among the highly enriched biological processes among downregulated DEGs in shoot tissues of PBW703.

FIGURE 6

FIGURE 7

Validation of DEG Genes Using Real-Time Quantitative Polymerase Chain Reaction Analysis

To validate the expression data obtained by RNA-seq, we performed RT-qPCR analysis of the expression of the 17 selected DEGs (Supplementary Table S5) with both upregulated and downregulated expression in both of the tissues and in both genotypes. Results were found in agreement with the RNA-seq-based gene expression pattern with minor variations in the log2 FC values (Figure 8).

FIGURE 8

Discussion

To find out the underlying genes involved in NUE pathways, a transcriptomics study was conducted in control and N stress conditions in root and shoot tissues of two cultivars PBW677 and PBW703. Differential gene expression revealed the key genes involved in photosynthesis, starch and sugar metabolism, and nitrogen metabolism in both the cultivars. In limited N conditions, N primarily affected plant morphology, accelerating plant flowering time, senesces, decreasing the height, and lowering biomass and harvest index. The overall stresses of the plant were evident from the spikelets per spike while the average seed test grain weight was comparable. This indicated the plants concentrated their energy to produced viable seeds. The nitrogen use efficiency was affected by supplied N levels. The increase of NUE was affected by N regimes in winter and spring wheat (Cormier et al., 2013; Nehe et al., 2018). Under nitrogen stress, the number of upregulated as well as downregulated genes was higher in roots as compared to shoots. This is in accordance with the previous study conducted by Curci et al. (2017) on durum wheat as roots are primary organs responsive to N stress. Also, it was found that the number of DEGs was higher in PBW677 as compared to PBW703. This might be due to the fact that PBW677 responds actively in N-deficient conditions compared to PBW703 as in the case reported by Sultana et al., in 2020 where medium NUE Spitfire had more DEGs in comparison to low NUE cultivar Volcani. The 97 common genes between root and shoot tissues were mostly related to plant hormone signal transduction, MAPK signaling pathway, starch and sucrose metabolism, glutathione metabolism, and chlorophyll and porphyrin metabolism. The 33 DEGs common in shoot tissues were mainly involved in lipid transport, amino acid metabolism, photosynthesis, carotenoid biosynthesis, starch and sucrose metabolism, and nucleotide metabolism. Several previous studies indicated that when plants were subjected to nitrogen stress they exhibited a wide range of responses including molecular expression and biochemical processes involving various genes and pathways. There is a strong relation between nitrogen and carbon metabolism as has been shown in various studies (Duan et al., 2018; Naliwajski and Skłodowska, 2018). Gene enrichment analysis based on Gene Ontology is the method by which key biological processes responsible for stress adaptation can be revealed, for example an N deficiency study in durum wheat reported photosynthesis, carbon metabolism, and nitrogen metabolism as the highly enriched biological processes (Curci et al., 2017). Under N limitation, the metabolic process and oxidation–reduction process in wheat seedlings were enriched significantly (Wang et al., 2019). Under nitrogen stress, catalytic activity, binding, and the metabolic and cellular process were highly enriched GO terms in potato roots, shoots, and stolons (Tiwari et al., 2020). This study revealed that upregulated DEGs in root tissues of PBW677 were involved in pathways like DNA binding, response to abiotic stimulus, and defense response while the downregulated genes were involved in biological processes like photosynthesis, photorespiration, and carbon fixation. While in shoot tissues, the highly enriched biological processes were carbohydrate metabolic process, cell wall biogenesis, and cellular polysaccharide metabolic process. In contrast to this, the upregulated DEGs in root tissues of PBW703 belonged to enriched biological processes like photosynthesis, carbon metabolism, glutamine biosynthetic process, and nitrogen compound metabolic processes while in shoot tissues, the enriched pathways were involved in nucleosome assembly, sulfur compound metabolic process, glutathione metabolic process, nitrate transport, etc. The reason behind the upregulated transcripts belonging to the pathways related to nucleosome assembly and chromatin packaging would be the mechanisms involved towards the stability of the epigenome under N stress in PBW703 (Sudan et al., 2018). N metabolism-related genes were identified in the annotated transcriptome data showing both up as well as downregulation. For example expression of low affinity nitrate transporter NRT1 was found to be upregulated in root tissues of PBW677 and downregulated in shoot tissues of PBW677 and PBW703 under N stress which usually increases the expression of transport systems for nitrate and ammonium (Crawford and Glass, 1998). In bread wheat under N stress, a high level of expression of the NRT1/PTR family was found (Sultana et al., 2020). Also, it was shown that in the low N-treated plants, downregulation of the NRT1/PTR family was observed which is related to low grain yield and grain protein content (Léran et al., 2015). Other genes playing a role in N absorption, assimilation, and remobilization like glutamate synthase (GLT), amino acid transporter (AAP), and asparagine synthase were upregulated in both the root and shoot tissues of PBW677 where as these genes were downregulated in PBW703. The upregulated expression of these genes can be related to the high NUE of PBW677 which is better adapted to N stress. Aminotransferases like serine glyoxylate aminotransferase (AGT), glyoxylate aminotransferase (GGAT1), and chloroplast glutamine synthetase (GLN) were upregulated in roots of PBW703 which is supported by previous reports conducted in sorghum where a transcriptomics study found more N assimilator genes in low NUE cultivars (Singh et al., 1973).

Transcription factors (TFs) are important key regulators that play a significant role in adaptation under environmental stresses (Shahzad et al., 2020). In our study we found a number of transcription factor families up and downregulated under N stress. The number of expressed TFs was higher in PBW677 than PBW703 which might be due to the fact that effective regulatory mechanisms exist under N stress in high NUE cultivar PBW677. Several TF families that were identified in response to N stress were HSF, MYB, WRKY, and ZINC finger. Transcriptomics analysis in durum wheat identified most of the WRKY family in response to N chronic stress in roots (Curci et al., 2017). This family is one of the largest families of plant-specific transcription factors that plays important roles in various abiotic stresses (drought, saline, alkali, temperature, and ultraviolet radiation (Li et al., 2020). Several MYB and bHLH TFs were reported to be involved in regulation of target genes under plant stress (Pireyre and Burow, 2015). In rice, expression of OSMYB48-1 was reported to be upregulated in tolerance toward abiotic stress probably via the regulation of stress-mediated ABA biosynthesis (Xiong et al., 2014).

Zinc finger proteins (ZEPs) enhance tolerance under abiotic stress (Zang et al., 2016). A high level of expression was observed for protein kinases in high NUE PBW677 (either up or downregulated). Protein kinases are known to play an important role in signal transduction and are regarded to be the central regulatory components to major environmental stresses such as drought, heat, cold, and pathogen attack (Wang et al., 2020). PKs play a role in hormone signaling, cell cycle growth, and nutrient signaling as well.

Nitrogen of the plant is invested in making protein and chlorophyll content of photosynthetic apparatus thus carbon metabolism and nitrogen metabolism are interconnected and carbon metabolism is dependent on nitrogen assimilation in plants (Foyer et al., 2001). In our study we found that a number of carbon metabolism genes were downregulated like genes involved in glycolysis, TCA cycle, and gluconeogenesis. This is similar to previous reports in which it was shown that N stress negatively effects plant carbohydrate metabolism (Rufty et al., 1988; Sultana et al., 2020). We found other stress-related genes to be highly expressed under N stress in high NUE cultivar PBW677 like GSTs, cytochrome P 450, E3 protein ligase, etc. Glutathione S transferases prevent cells from oxidative damage by quenching reactive molecules with the addition of glutathione (GSH) (Kumar and Trivedi, 2018).

Conclusion

RNA sequencing of two contrasting cultivars PBW677 and PBW703 for NUE helped in revealing candidate genes which could be utilized in future breeding programs focused on reducing the use of nitrogen fertilizers. The majority of the genes belonged to transcription factor families; protein kinases and stress-related nitrogen metabolism were found to be highly expressive in wheat cultivar PBW677 which might explain its behavior under N stress. There was a difference in highly enriched pathways responsive to nitrogen stress in the contrasting cultivars which might be the cause for their different behavior towards N stress. The 13 common genes showed alternative expression patterns in PBW677 and PBW703 and could be the potential candidates for high NUE-targeted breeding.

References

• 1

  AndrewsS.KruegerF.Seconds-PichonA.BigginsF.WingettS. (2015). FastQC. A Quality Control Tool for High Throughput Sequence Data. Babraham Bioinformatics. Babraham Inst. Available at: https://www.bibsonomy.org/bibtex/f230a919c34360709aa298734d63dca3 (Accessed January 3, 2022).

  • Google Scholar

• 2

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics30, 2114–2120. 10.1093/bioinformatics/btu170

  • CrossRef
  • Google Scholar

• 3

  CormierF.FaureS.DubreuilP.HeumezE.BeauchêneK.LafargeS.et al (2013). A Multi-Environmental Study of Recent Breeding Progress on Nitrogen Use Efficiency in Wheat (Triticum A L.). Theor. Appl. Genet.126, 3035–3048. 10.1007/s00122-013-2191-9

  • CrossRef
  • Google Scholar

• 4

  CrawfordN. M.GlassA. D. M. (1998). Molecular and Physiological Aspects of Nitrate Uptake in Plants. Trends Plant Sci.3, 389–395. 10.1016/S1360-1385(98)01311-9

  • CrossRef
  • Google Scholar

• 5

  CurciP. L.Aiese CiglianoR.ZuluagaD. L.JanniM.SanseverinoW.SonnanteG. (2017). Transcriptomic Response of Durum Wheat to Nitrogen Starvation. Sci. Rep.7, 1–14. 10.1038/s41598-017-01377-0

  • CrossRef
  • Google Scholar

• 6

  DuanW.WangQ.ZhangH.XieB.LiA.HouF.et al (2018). Comparative Study on Carbon-Nitrogen Metabolism and Endogenous Hormone Contents in normal and Overgrown Sweetpotato. South Afr. J. Bot.115, 199–207. 10.1016/J.SAJB.2017.11.016

  • CrossRef
  • Google Scholar

• 7

  EastinE. F. (2008). Use of an AutoAnalyzer for Total Nitrogen Determination in Plants. Commun. Soil Sci. Plant Anal.9, 107–113. 10.1080/00103627809366792

  • CrossRef
  • Google Scholar

• 8

  FoyerC. H.Ferrario-MéryS.NoctorG. (2001). “Interactions between Carbon and Nitrogen Metabolism,” in Plant Nitrogen. Editors LeaP. J.Morot-GaudryJ. F. (Berlin, Heidelberg: Springer), 237–254. 10.1007/978-3-662-04064-5_9

  • CrossRef
  • Google Scholar

• 9

  GelliM.DuoY.KondaA. R.ZhangC.HoldingD.DweikatI. (2014). Identification of Differentially Expressed Genes between Sorghum Genotypes with Contrasting Nitrogen Stress Tolerance by Genome-Wide Transcriptional Profiling. BMC Genomics15, 179. 10.1186/1471-2164-15-179/TABLES/6

  • CrossRef
  • Google Scholar

• 10

  GoodA. G.ShrawatA. K.MuenchD. G. (2004). Can Less Yield More? Is Reducing Nutrient Input into the Environment Compatible with Maintaining Crop Production?Trends Plant Sci.9, 597–605. 10.1016/J.TPLANTS.2004.10.008

  • CrossRef
  • Google Scholar

• 11

  HaoQ. N.ZhouX. A.ShaA. H.WangC.ZhouR.ChenS. L. (2011). Identification of Genes Associated with Nitrogen-Use Efficiency by Genome-wide Transcriptional Analysis of Two Soybean Genotypes. BMC Genomics12, 525. 10.1186/1471-2164-12-525

  • CrossRef
  • Google Scholar

• 12

  HitzK.ClarkA. J.Van SanfordD. A. (2017). Identifying Nitrogen-Use Efficient Soft Red winter Wheat Lines in High and Low Nitrogen Environments. Field Crops Res.200, 1–9. 10.1016/j.fcr.2016.10.001

  • CrossRef
  • Google Scholar

• 13

  KantS.BiY.-M.RothsteinS. J. (2011). Understanding Plant Response to Nitrogen Limitation for the Improvement of Crop Nitrogen Use Efficiency. J. Exp. Bot.62, 1499–1509. 10.1093/JXB/ERQ297

  • CrossRef
  • Google Scholar

• 14

  KumarS.TrivediP. K. (2018). Glutathione S-Transferases: Role in Combating Abiotic Stresses Including Arsenic Detoxification in Plants. Front. Plant Sci.9, 751. 10.3389/FPLS.2018.00751/BIBTEX

  • CrossRef
  • Google Scholar

• 15

  LéranS.EdelK. H.PerventM.HashimotoK.Corratgé-FaillieC.OffenbornJ. N.et al (2015). Nitrate Sensing and Uptake in Arabidopsis Are Enhanced by ABI2, a Phosphatase Inactivated by the Stress Hormone Abscisic Acid. Sci. Signal.8 (375), ra43. 10.1126/SCISIGNAL.AAA4829

  • CrossRef
  • Google Scholar

• 16

  LiH.DurbinR. (2009). Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics25, 1754–1760. 10.1093/BIOINFORMATICS/BTP324

  • CrossRef
  • Google Scholar

• 17

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al (2009). The Sequence Alignment/Map Format and SAMtools. Bioinformatics25, 2078–2079. 10.1093/bioinformatics/btp352

  • CrossRef
  • Google Scholar

• 18

  LiW.PangS.LuZ.JinB. (2020). Function and Mechanism of WRKY Transcription Factors in Abiotic Stress Responses of Plants. Plants (Basel)9, 1–15. 10.3390/plants9111515

  • CrossRef
  • Google Scholar

• 19

  Masclaux-DaubresseC.Daniel-VedeleF.DechorgnatJ.ChardonF.GaufichonL.SuzukiA. (2010). Nitrogen Uptake, Assimilation and Remobilization in Plants: Challenges for Sustainable and Productive Agriculture. Ann. Bot.105, 1141–1157. 10.1093/AOB/MCQ028

  • CrossRef
  • Google Scholar

• 20

  MengX.LiangZ.DaiX.ZhangY.MahboubS.NguD. W.et al (2021). Predicting Transcriptional Responses to Cold Stress across Plant Species. Proc. Natl. Acad. Sci. U. S. A.118 (10), e2026330118. 10.1073/PNAS.2026330118

  • CrossRef
  • Google Scholar

• 21

  MollR. H.KamprathE. J.JacksonW. A. (1982). Analysis and Interpretation of Factors Which Contribute to Efficiency of Nitrogen Utilization1. Agron. J.74, 562–564. 10.2134/AGRONJ1982.00021962007400030037X

  • CrossRef
  • Google Scholar

• 22

  NaliwajskiM. R.SkłodowskaM. (2018). The Relationship between Carbon and Nitrogen Metabolism in Cucumber Leaves Acclimated to Salt Stress. PeerJ6, e6043. 10.7717/peerj.6043

  • CrossRef
  • Google Scholar

• 23

  NeheA. S.MisraS.MurchieE. H.ChinnathambiK.FoulkesM. J. (2018). Genetic Variation in N-Use Efficiency and Associated Traits in Indian Wheat Cultivars. Field Crops Res.225, 152–162. 10.1016/j.fcr.2018.06.002

  • CrossRef
  • Google Scholar

• 24

  PaolacciA. R.TanzarellaO. A.PorcedduE.CiaffiM. (2009). Identification and Validation of Reference Genes for Quantitative RT-PCR Normalization in Wheat. BMC Mol. Biol10, 1–27. 10.1186/1471-2199-10-11

  • CrossRef
  • Google Scholar

• 25

  PireyreM.BurowM. (2015). Regulation of MYB and bHLH Transcription Factors: A Glance at the Protein Level. Mol. Plant8, 378–388. 10.1016/j.molp.2014.11.022

  • CrossRef
  • Google Scholar

• 26

  RuftyT. W.Jr.HuberS. C.VolkR. J. (1988). Alterations in Leaf Carbohydrate Metabolism in Response to Nitrogen Stress. Plant Physiol.88, 725–730. 10.1104/PP.88.3.725

  • CrossRef
  • Google Scholar

• 27

  RuuskaS. A.LewisD. C.KennedyG.FurbankR. T.JenkinsC. L. D.TabeL. M. (2008). Large Scale Transcriptome Analysis of the Effects of Nitrogen Nutrition on Accumulation of Stem Carbohydrate Reserves in Reproductive Stage Wheat. Plant Mol. Biol.66, 15–32. 10.1007/S11103-007-9249-5/FIGURES/7

  • CrossRef
  • Google Scholar

• 28

  ShahzadR.JamilS.AhmadS.NisarA.AminaZ.SaleemS.et al (2020). Harnessing the Potential of Plant Transcription Factors in Developing Climate-Smart Crops: Future Prospects, Challenges, and Opportunities. Preprints. 10.20944/PREPRINTS202010.0532.V1

  • CrossRef
  • Google Scholar

• 29

  ShamshadM. (2019). Genetic and Physiological basis of Nitrogen uptake and utilization in wheat (Triticum aestivum L.). PhD. thesis. Ludhiana, India: Punjab Agricultural University. Available at: https://krishikosh.egranth.ac.in/handle/1/5810137744

  • Google Scholar

• 30

  SinghT.PalegI.AspinallD. (1973). Stress Metabolism I. Nitrogen Metabolism and Growth in the Barley Plant during Water Stress. Aust. Jnl. Bio. Sci.26, 45–56. 10.1071/BI9730045

  • CrossRef
  • Google Scholar

• 31

  SudanJ.RainaM.SinghR. (2018). Plant Epigenetic Mechanisms: Role in Abiotic Stress and Their Generational Heritability. 3 Biotech.8 (3), 172. 10.1007/S13205-018-1202-6

  • CrossRef
  • Google Scholar

• 32

  SultanaN.IslamS.JuhaszA.MaW. (2021). Wheat Leaf Senescence and its Regulatory Gene Network. Crop J.9, 703–717. 10.1016/J.CJ.2021.01.004

  • CrossRef
  • Google Scholar

• 33

  SultanaN.IslamS.JuhaszA.YangR.SheM.AlhabbarZ.et al (2020). Transcriptomic Study for Identification of Major Nitrogen Stress Responsive Genes in Australian Bread Wheat Cultivars. Front. Genet.11, 583785. 10.3389/FGENE.2020.583785

  • CrossRef
  • Google Scholar

• 34

  TiongJ.SharmaN.SampathR.MacKenzieN.WatanabeS.MetotC.et al (2021). Improving Nitrogen Use Efficiency through Overexpression of Alanine Aminotransferase in Rice, Wheat, and Barley. Front. Plant Sci.12, 29. 10.3389/FPLS.2021.628521/BIBTEX

  • CrossRef
  • Google Scholar

• 35

  TiwariJ. K.BucksethT.ZintaR.SaraswatiA.SinghR. K.RawatS.et al (2020). Transcriptome Analysis of Potato Shoots, Roots and Stolons under Nitrogen Stress. Sci. Rep.10, 1152. 10.1038/s41598-020-58167-4

  • CrossRef
  • Google Scholar

• 36

  TrapnellC.RobertsA.GoffL.PerteaG.KimD.KelleyD. R.et al (2012). Differential Gene and Transcript Expression Analysis of RNA-Seq Experiments with TopHat and Cufflinks. Nat. Protoc.7, 562–578. 10.1038/nprot.2012.016

  • CrossRef
  • Google Scholar

• 37

  VijayalakshmiP.VishnukiranT.Ramana KumariB.SrikanthB.Subhakar RaoI.SwamyK. N.et al (2015). Biochemical and Physiological Characterization for Nitrogen Use Efficiency in Aromatic Rice Genotypes. Field Crops Res.179, 132–143. 10.1016/j.fcr.2015.04.012

  • CrossRef
  • Google Scholar

• 38

  WangJ.SongK.SunL.QinQ.SunY.PanJ.et al (2019). Morphological and Transcriptome Analysis of Wheat Seedlings Response to Low Nitrogen Stress. Plants8, 98. 10.3390/PLANTS8040098

  • CrossRef
  • Google Scholar

• 39

  WangP.HsuC.-C.DuY.ZhuP.ZhaoC.FuX.et al (2020). Mapping Proteome-Wide Targets of Protein Kinases in Plant Stress Responses. Proc. Natl. Acad. Sci. USA117, 3270–3280. 10.1073/pnas.1919901117

  • CrossRef
  • Google Scholar

• 40

  XiongH.LiJ.LiuP.DuanJ.ZhaoY.GuoX.et al (2014). Overexpression of OsMYB48-1, a Novel MYB-Related Transcription Factor, Enhances Drought and Salinity Tolerance in rice. PLoS One9, e92913. 10.1371/journal.pone.0092913

  • CrossRef
  • Google Scholar

• 41

  XuG.FanX.MillerA. J. (2012). Plant Nitrogen Assimilation and Use Efficiency. Annu. Rev. Plant Biol.63, 153–182. 10.1146/ANNUREV-ARPLANT-042811-105532

  • CrossRef
  • Google Scholar

• 42

  XuZ.YuZ.ZhaoJ. (2013). Theory and Application for the Promotion of Wheat Production in China: Past, Present and Future. J. Sci. Food Agric.93, 2339–2350. 10.1002/JSFA.6098

  • CrossRef
  • Google Scholar

• 43

  ZamboniA.AstolfiS.ZuchiS.PiiY.GuardiniK.TononiP.et al (2014). Nitrate Induction Triggers Different Transcriptional Changes in a High and a Low Nitrogen Use Efficiency maize Inbred Line. J. Integr. Plant Biol.56, 1080–1094. 10.1111/JIPB.12214

  • CrossRef
  • Google Scholar

• 44

  ZangD.LiH.XuH.ZhangW.ZhangY.ShiX.et al (2016). An Arabidopsis Zinc finger Protein Increases Abiotic Stress Tolerance by Regulating Sodium and Potassium Homeostasis, Reactive Oxygen Species Scavenging and Osmotic Potential. Front. Plant Sci.7, 1–14. 10.3389/fpls.2016.01272

  • CrossRef
  • Google Scholar