1. Introduction
The foundation of modern agriculture is based on the use of high-yielding cultivars that are responsive to heavy doses of nitrogenous-fertilizers because nitrogen (N) is the single most critical macronutrient used for commercial crop production. The application of nitrogenous fertilizers has doubled with the development of high-yielding cultivars to increase food production for the growing global population [1]. Despite the ecological unsustainability of heavy N fertilization in agriculture, this practice continues unabated. Most applied nitrogen is lost to the environment raising both environmental and health concerns due to air and water pollution [2]. As a result, nitrate contamination of ground and surface water has become both a regulatory and social issue. Over-fertilization with N leads to increased disease incidence and yield loss due to lodging. An increased efficiency of N fertilizers, which is the single largest variable cost for crop production, can enhance the profitability of farming by reducing fertilizer usage [3]. A reduction in use of N fertilizers can be achieved by developing sustainable management practices and crop cultivars with improved nitrogen use efficiency (NUE) to minimize the impact on the global environment [2]. In contrast to the earlier green revolution strategy that focuses on developing genotypes responding to higher doses of N, current research is now directed at genotypes that perform well under N limited environments without compromising crop productivity. In rice growing areas of the USA, an alternate wetting and drying technique is slowly gaining popularity to reduce the cost of cultivation by reducing water usage. However, this strategy will require the development of rice cultivars with improved water and nitrogen use efficiency to reduce N loss during the process. Improvements in N utilization can be enhanced through an improved nitrogen transport, assimilation, and remobilization [4]. Therefore, deciphering the molecular basis of the natural variation for sensing and responding to N availability has great potential to design crop varieties which can ensure a healthy balance between the crop yield and environmental consequences [5,6].
More than half of the world’s population consume rice (Oryza sativa L.) as a staple food. Regardless of the physiological and genetic complexity of N stress tolerance, rice is an attractive model organism due to its small genome size and the availability of many genomic and germplasm resources [7,8]. As a result of domestication and artificial selection, there is a significant difference in N uptake and assimilation between the two subspecies of rice [9]. The indica varieties have long been known to be superior in NUE compared to japonica varieties [10,11]. Based on this premise, Hu et al. (2015) [9] cloned NTR1.1B, which explained the divergence in nitrogen use between indica and japonica subspecies and the superiority of the indica allele. Therefore, the mining and utilization of superior variants from indica germplasm should help design varieties with improved NUE [12]. Although there were studies focusing on mapping and gene expression in rice, none of these studies involved the USA adapted rice genotypes. Since it is widely known that quantitative trait loci (QTL), gene expression and stability varies with the genetic backgrounds due to the difference in allelic combinations, there is a need to discover the genetic determinants that differentiate US cultivars from the indica cultivars for use in rice breeding programs in the United States.
Genetic engineering has been advocated to manipulate nitrogen assimilatory pathway genes [13,14,15,16], but no cultivar has reached farmers’ fields to date despite the claims of an increased biomass and grain yield in transgenic plants. The most plausible approach is to exploit natural genetic variations in existing germplasm [4]. Specifically, the superior ability of indica rice for N uptake and mobilization [9,10] offers a unique opportunity to design nitrogen use-efficient cultivars by providing insights into the molecular genetic basis of NUE. However, the real gap between research and commercial product development lies with the difficulty of assembling all beneficial alleles in an advanced breeding line.
RNA-seq is a highly sensitive and powerful transcriptomic tool for identification of differentially expressed genes (DEGs) [17]. It has been applied to identify genes responsive to variable N regime. In most studies, only one genotype is used irrespective of its level of NUE [18,19,20]. Yang et al. [19] identified 1650 DEGs in both root and leaf sheaths under N deficient condition and 86 were transcription factors (TF) from 28 families such as AP2/EREBP, WRKY, GATA, Dof, and MADS. A few studies used comparative RNA-seq approaches using two rice varieties [21,22]. A transcriptome profiling study involving IR64 and N22 under optimal and N starved conditions revealed that most DEGs were associated with starch and chloroplast metabolisms and signaling, and many of these genes were localized within the NUE-related QTL intervals [21]. A comparison of the shoot apical meristem transcriptome between an indica variety “YD6” and a japonica variety “Nipponbare” under a variable N regime revealed that majority of the DEGs are involved in general stress responses, stimulus responses, and hormonal signaling process [22].
Although many genes involved in N uptake and assimilation are known, information on transcriptional regulation and signal transduction processes, associated with improvements in NUE, is limited. Many QTLs from NUE mapping studies in rice and other crops overlapped with genes that were involved in nitrogen uptake and assimilation [23,24,25]. For a complex trait like NUE, which is controlled by many genes with several interacting environmental factors, we applied a genome-wide expression profiling technique, using a japonica genotype “Bengal” and an indica genotype “Pokkali”, with contrasting responses to nitrogen availability to discover candidate genes with the potential to improve the NUE in rice. Our analysis revealed many regulatory genes and genes involved in different pathways which can be used for developing strategies to improve the NUE in rice.
2. Results 2.1. Chlorate Experiment
Chlorate and nitrate are analogs and the uptake of chlorate during the seedling stage results in growth inhibition due to its toxicity to plants. Hence, the nitrogen uptake capacity was inferred from the shoot length measurements in a chlorate assay. Among the eight genotypes used for this assay, the shoot inhibition rate was significantly higher in indicaa genotypes compared with others (Figure 1 and Figure S1). Pokkali and Bengal showed the highest (48%) and the least (2%) shoot reduction, respectively. Based on the result of this initial experiment, Pokkali (indica) and Bengal (japonica) were selected for N stress and RNA-seq experiments.
2.2. Nitrogen Stress Experiment
Limiting the amount of N in nutrient solution was able to differentiate the agro-morphological responses between Pokkali (PK) and Bengal (BG). Generally, a reduction in shoot length, leaf greenness, and biomass was observed in both genotypes (Figure 2). However, Pokkali showed a smaller reduction in shoot length, chlorophyll content, and biomass than Bengal.
Changes in root architecture under nitrogen starvation was evident in Pokkali and Bengal (Table 1 and Figure 3). The response to N stress (low N) lead to an increased lateral root length, total number of roots, total root surface area, and root volume except for root diameter/thickness in both genotypes. The root attributes were more pronounced in Pokkali compared with Bengal under normal and low N. The root length and root counts were positively correlated with the total root surface area and root volume (Table S1 and Figure S2).
2.3. Transcriptome Analysis under Nitrogen Stress
RNA sequencing of root tissues in different treatments resulted in 421 and 394 million raw reads in Bengal and Pokkali, respectively (Table S2). Treatments were low N (LN), full N (FN) and 1 h after the transfer from low N to full N. All treatments had three biological replicates. The mean assembly percentage to the reference genome was 70 and 83 for Bengal and Pokkali, respectively. Only high-quality reads were used for all downstream analyses. In total, 208,719 and 208,338 transcripts were identified in Bengal and Pokkali, respectively. On average, there were 2652 and 2371 novel transcripts in Bengal and Pokkali, respectively.
The differentially expressed genes were identified for each pairwise comparisons of Pokkali and Bengal, and N treatments (Figure S3). The number of common, treatment-specific, and genotype-specific DEGs are illustrated in Figure 4 and Figure 5. The total number of DEGs in response to low N was 2937 which was comprised of 1753 (611 upregulated and 1142 downregulated), 1150 (136 upregulated and 1014 downregulated), and 756 (98 upregulated and 658 downregulated) for PKLN (Pokkali-low nitrogen) vs. BGLN (Bengal-low nitrogen), PKLN vs. PKFN (Pokkali-full nitrogen) and BGLN vs. BGFN (Bengal-full nitrogen), respectively (Figure 4A and Figure 5). On the other hand, there were 2183 early N responses (1 h after transfer from low N to full N) DEGs which were comprised of 1549 (531 up and 1018 down), 492 (393 up and 99 down), and 424 (264 up and 160 down) for PK1H (Pokkali 1 h after transfer from low N to full N) vs. BG1H (Bengal 1 h after transfer from low N to full N), PK1H vs. PKLN, and BG1H vs. BGLN, respectively (Figure 4B and Figure 5). Overall, the number of DEGs in Pokkali was higher than in Bengal during low N stress and early recovery. The DEGs in both genotypes grown under low N condition (PKLN vs. PKFN, BGLN vs. BGFN and PKLN vs. BGLN) were mostly downregulated, whereas DEGs were mostly upregulated in early N recovery (PK1H vs. PKLN, BG1H vs. BGLN and PK1H vs. BG1H (Figure 5)).
Furthermore, the Pokkali-specific N response DEGs were 112 (18 up and 94 down), and 38 (23 up and 15 down) for low N [(PKLN vs. PKFN ∩ PKLN vs. BGLN) – (BGLN vs. BGFN)], and early N recovery [(PK1H vs. PKLN ∩ PK1H vs. BG1H) – (BG1H vs. BGLN)], respectively (Figure 5 and Figure 6). Due to the overlapping of some genes, there were 148 Pokkali-specific DEGs (Table S3). Few genes were common between both N conditions. For example, LOC_Os06g05450 and LOC_Os12g28015 were both regulated during N stress and early recovery. In addition, a literature review of Pokkali-specific DEGs for gene ontology (GO) functions indicated that 9% of genes were involved in plant growth and development whereas the rest were associated with regulation and metabolism (15%), stress responses (17%) and unknown functions (59% (Figure S4)).
2.4. Gene Ontology (GO) and Pathway Analysis of Differentially Expressed Genes under Nitrogen Stress
GO enrichment was performed to ascertain the functions of DEGs and the related biological processes, whereas the pathway analysis was done to determine the major and specific plant pathways involved in each N treatment. Combined Pokkali and Bengal DEGs from each treatment were used in a singular enrichment analysis (SEA) and Plant Reactome pathway analysis. In the SEA analysis, the N-responsive DEGs from both Pokkali and Bengal were classified within the three domains: “biological process”, “molecular function” and “cellular component” (Figure 7). There were 11 and 6 significant GO terms identified for low N and early N recovery, respectively. The “biological process” had the greatest number of enriched GO terms among the three classified domains, and this was true for both low and recovery N treatments. The “molecular function” domain revealed “oxygen binding”, “binding”, and “catalytic activity” GO terms in low N while the early N response showed “binding” GO terms. These differences in GO terms showed diversity in response to different N conditions. In addition, the GO enrichment also revealed that only low N conditions had the “cellular component” terms: “the extracellular region”, “external encapsulating structure”, and “cell wall”. These confer a root morphological response for foraging N under low N condition as indicated by root induction.
Plant Reactome pathway analyses of DEGs showed major plant pathways involved in “cellular process”, “growth and developmental process”, “metabolism and regulation”, “abiotic and biotic stimuli and stresses” (Table 2). The majority of the DEGs in low N and early N recovery (84% and 81%, respectively) were involved in the “metabolism and regulation” plant pathway (Figure S5). The common basic “metabolism and regulation” pathways between low N and immediate recovery conditions were hormone signaling, transport and metabolism, amino acid metabolism, carbohydrate metabolism, cofactor biosynthesis, cytokinin biosynthesis, inorganic nutrient metabolism, and secondary metabolism. Similarly, the “cellular process”, “growth and developmental process” and “response to abiotic and stimuli and stresses” had common basic pathways as well. In addition, unique pathways for each N response were observed: vegetative structure developments, amine and polyamine biosynthesis, detoxification, fatty acid and lipid metabolism, generation of precursor metabolites and energy, photorespiration, response to phosphate deficiency, recognition of fungal and bacterial pathogens and immunity response for a low N response, and DNA replication and activation of the pre-replicative complex for early N recovery. This indicates that there are common basic pathways for both Pokkali and Bengal in response to low N and early N recovery and also unique pathways to cope with each N condition.
2.5. Expression Pattern of Nitrogen Utilization and Long Distance Signaling
The nitrogen availability in the nutrient solution affected certain N utilization genes in the roots of Pokkali and Bengal which were essential for N uptake, translocation, and assimilation. Some of these genes were found to be differentially expressed (Table 3). The N uptake genes were nitrate reductases (LOC_Os08g36480, LOC_Os08g36500, and LOC_Os02g53130). Nitrate reductase 1 and 2 were not differentially expressed in Pokkali, regardless of N conditions. However, these reductases in Bengal were upregulated at early recovery. The NADH/NADPH-dependent NO3-reductase 2 was differentially expressed in Pokkali and Bengal at low N, and recovery and both cultivars showed similar regulation pattern. Among the N transport/translocation-related genes, the ammonium transporter (LOC_Os02g40730) and high-affinity nitrate transporter 2.1 (LOC_Os02g02170) showed different expression patterns in Pokkali and Bengal under low N conditions, whereas the other nitrate transporter genes (LOC_Os10g40600, LOC_Os02g02190, and LOC_Os01g50820) showed similar patterns for both cultivars. Lastly, a few N assimilation genes, such as glutamine synthetase 1;2 (LOC_Os03g12290) and NADH-GOGAT (LOC_Os01g48960), were regulated in Pokkali and Bengal with similar expression patterns.
Certain proteins such as C-terminally encoded peptides (CEPs) and CLV3/ESR-related (CLE) peptides were known to be involved in long distance systemic signaling in plants during N deprivation [26]. Long distance signaling between the shoot and roots is important in the root growth response during N stresses. In our study, a few CEPs (LOC_Os05g11620, LOC_Os05g11580, and LOC_Os06g43080) and CLEs (LOC_Os01g23705, LOC_Os01g32560 and LOC_Os03g48570) showed differential expressions in Pokkali and Bengal roots at various N treatments (Table S4).
2.6. Transcription Factors and Signaling for Nitrogen Utilization and Root Growth and Development Genes
Transcription factors are known to control plant growth and modulate responses to biotic and abiotic stimuli through the regulation of gene expression. There were 152 and 125 TFs encompassing 30 TF families identified in root tissues in Pokkali and Bengal, respectively (Table 4). The general view of the regulation pattern of TFs indicated that most of the TFs were downregulated in low N conditions while they are upregulated in early N recovery stage (Figure S6). An analysis of DEGs from different transcription factor families revealed that members of ethylene responsive factor (ERF), WRKY, basic helix-loop-helix (bHLH), MYB, C2H2 zinc finger, and NAC families were most abundant with 14.9%, 14.4%, 8.9%, 8.9%, 7.9% and 7.9% of TFs, respectively (Figure S7).
2.7. Phytohormones Related to N Stress and Root Growth and Development
Plant hormones play a significant role in the regulation of plant development and environmental responses. Based on Plant Reactome annotations, several DEGs with phytohormone signaling functions in various N conditions were identified, suggesting the important role of plant hormones in processes involved in improving the NUE (Figure 8).
2.8. Differentially Expressed Genes Overlapping the Nitrogen Stress-Related QTL
We found several DEGs co-localized with the selected QTLs for low N responses (Table S5), and only 81 of these were Pokkali-specific (Figure 9). These co-localized DEGs consisted of loci that lie within, and some overlapping with other QTLs. For example, LOC_Os03g20680 lies within NUE and root number (rtnb) QTLs and LOC_Os11g40690 lie within chlorate resistance (kclo3rs), rtnb, and root length (rtlg) QTLs.
2.9. Alternate Splicing Events under Nitrogen Stress
A total of 39,668 alternative splicing (AS) events were identified in Bengal (19,753) and Pokkali (19,916) (Figure 10A). The most predominant AS event was A3SS with nearly 29% of total events in both genotypes followed by intron retention (IR (BG: 23.4%, PK: 23.3%)), exon skipping (ES (BG: 13.0%, PK: 12.9%)), alternative 5′ splice site (A5SS (BG:12.9%, PK:12.8%)), A5SS or alternative 3′ splice site (A3SS (BG: 2.8%, PK:3.0%)), IR1 + IR2 (BG: 2.7%, PK: 2.5%), ES1 + ES2 (BG: 1.7%, PK: 1.6%), A5SS + A3SS (BG: 1.6%, PK: 1.7%), mutually exclusive exon (MXE (BG: 0.7%, PK: 0.6%)), IR1 or IR2 (BG: 0.7%, PK: 0.8%), and A5SS + ES + A3SS (0.2%).The analysis of the combined assembly of genotypes and treatments showed that Bengal had slightly more AS events in FN and LN treatments compared to Pokkali while Pokkali had higher AS in 1 h after recovery than Bengal (Figure 10B). Looking at the genotype per treatment comparison, Bengal shows the highest number of AS events in low nitrogen treatment followed by a decrease in 1-h recovery treatment. Out of the 148 Pokkali-specific differentially expressed genes, 17 exhibited alternative splicing events (Table S6). Overall, there were differences in the types of AS events between different treatments within a genotype or between both genotypes.
2.10. Gene Expression Validation of Selected Genes by Quantitative Reverse Transcription PCR (qRT-PCR)
Eleven genes were selected based on their response to nitrogen stress, growth, yield, and root development to validate the expression pattern obtained from RNA-seq data using qRT-PCR under different treatments (full N, low N, and 1-h recovery) in both genotypes. After normalizing the Ct values using the internal control gene EF1α (Figure S8), a fold change calculation was used for a comparison between genotypes (i.e., full nitrogen treatment as a reference sample for low nitrogen while low nitrogen treatment was used as reference sample for 1-h recovery treatment). Pearson’s correlation coefficient (r = 0.71), based on the log2 fold change expression values of the selected genes, showed a positive and linear relationship between qRT-PCR and RNA-seq data indicating the reliability of RNA-seq results (Figure S9).
3. Discussion
Since a bulk of N from the applied fertilizers is lost to environment, reducing fertilizer use is necessary for sustainable agriculture. This can be accomplished by enhancing crop plants’ ability to efficiently uptake and assimilate N so that yield can be maintained with reduced fertilizer application. Therefore, an improved understanding of plants’ response N stress condition at the molecular level is needed to make such genetic improvements. In rice, RNA-seq studies have been conducted to reveal differentially regulated genes related to tillering [22] and carbon and nitrogen metabolism [20] in response to nitrogen availability. In this study, we used an indica genotype “Pokkali” and a japonica genotype “Bengal” to focus on root characteristics due to its pivotal role in N uptake [27]. We applied a comparative transcriptomic approach to differentiate transcriptional responses in the root tissues of these two genotypes under chronic nitrogen stress, as well as after the transfer to the optimal nitrogen condition and to identify differentially expressed genes which may have the potential to improve NUE. We hypothesize that variations in the nitrogen uptake capacity are due to the differential expression of genes which regulate the root architecture in plants.
Our results on chlorate assay (Figure 1) and seedling parameters (Figure 2 and Figure 3) showed a clear superiority of indica genotypes over japonica genotypes, with respect to nitrogen uptake and utilization thus confirming earlier studies [9,12,28]. These results prompted us to use Pokkali (indica) and Bengal (japonica) for a transcriptomic analysis of root tissues. The chlorate assay was supplemented with the N starvation study that showed a better agro-morphological response in Pokkali compared to Bengal. Pokkali was earlier reported to have a higher NUE [29]. Responses to nitrogen stress brought about changes in the root morphology and expression of genes related to nitrogen uptake and/or assimilation. A previous study on root architectural changes in response to variable N applications in different rice genotypes demonstrated the importance of root length and root number in enhancing the N foraging ability [30]. Hence, the low nitrogen availability induced root growth and development promoting more and longer lateral roots which resulted in a higher total root surface area and volume in Pokkali compared to Bengal. Our observations on the involvement of genes in root growth and development, as well as their overlapping with the QTL regions controlling root architecture, indicated that these genes may be responsible for the efficient uptake of N under the limited nitrogen condition in Pokkali.
Among the DEGs, there were several TFs and signaling genes involved in the root growth and development under low N conditions. Transcription factors belonging to the ERFs, WRKYs, bHLHs, MYBs, CH2H2s, and NACs play a multifaceted and overlapping role in the regulation of cell growth, development and proliferation, phytohormone signaling, biotic and abiotic stimuli/stress responses, regulation, and metabolism [31,32,33,34,35,36]. The genes involved in phytohormone signaling proteins/pathways control physiological and morphological responses under N stress [37,38]. Cytokinins, auxins, and gibberellin were involved in the network of nitrate-signaling pathways in plant developmental processes [39]. In addition, there were proteins such as CEPs and CLEs which were regulated under N stress and recovery. These proteins were involved in long distance shoot-to-root signaling and the systemic responses of root growth and development, specifically during N stress [26,39,40]. Some CEPs regulate the expression of NRT2.1 in Arabidopsis, suggesting interactions with core N uptake apparatus [41]. Therefore, there may be crosstalk between the physiological (nitrate transport) and morphological (root) responses under N stress and recovery involving phytohormone signaling genes, TFs, signal proteins, and transporters [26,39].
Examination of the root DEGs between Pokkali and Bengal showed more genes in Pokkali than in Bengal under N starvation. Likewise, the number of DEGs during early N recovery was also higher in Pokkali than in Bengal roots. It was not just the involvement of more DEGs, but also that the allelic differences resulting from different genetic backgrounds may be responsible for an increased N uptake efficiency in Pokkali compared to Bengal. There is a sizable number of genes overlapping between both genotypes observed in our study compared to an earlier study [21]. These common genes may be involved in the inherent mechanisms associated with uptake regardless of genotypic differences.
Nitrogen stress conditions caused NRT2.1 (LOC_Os02g02170) to be profoundly more modulated in Pokkali than in Bengal roots, while the nitrogen assimilation-mobilization genes, such as ammonia transporters, nitrate reductases, glutamine synthetases, and glutamine oxoglutarate aminotransferases (GOGAT), were either unchanged or downregulated under low N conditions (Table 3). The nitrate transporters in Pokkali played a major role in enhancing the uptake of N, particularly under limited N conditions, compared to Bengal. In addition, the N assimilation-mobilization activity was reduced as the latter processes were highly regulated by feedback inhibition depending on the nitrogen availability of plants. The downregulation of nitrate reductase, ammonia transporters, and glutamine synthetase could be due to the regulation of the substrate for nitrogen assimilation, particularly during chronic N stress.
There were several Pokkali-specific DEGs that might have contributed to Pokkali’s advantage over Bengal under both low and optimal N conditions. These DEGs were involved in stress responses, growth and development, regulation and metabolism. Some of these DEGs co-localized with QTLs related to root growth and development, chlorate-resistance, and NUE. These observations suggested that a combination of physiological and root morphological adaptations are responsible for N uptake by roots in rice. This has been demonstrated in case of the OsNRT2.1 which was involved in NO3− uptake under a low external NO3− concentration [42] and during root development [26,43]. The OsNRT2.1-overexpressing rice plant not only enhanced vegetative growth under low N conditions [42], but also enhanced the yield compared to the wild type under alternating wet and dry conditions [44] and drought stress [45] by increasing nitrogen uptake. This gene also promoted root elongation via auxin transport [43]. Although both genotypes promoted root architecture systems that could have aided in N foraging under low N, the root architecture is more pronounced in Pokkali than in Bengal. Another gene, OsDUR3 (LOC_Os10g42960), a high-affinity urea transporter, was upregulated under low N, confirming earlier studies [21,46]. Since OsDUR3 improved the rice yield under nitrogen deficient conditions by increasing urea uptake [47], the increased expressions of OsDUR3 should help in N foraging. Therefore, the upregulation of both OsDUR3 and OsNRT2.1 may be enhancing N foraging in the form of either nitrate or urea.
Two important Pokkali-specific differentially expressed genes were Nitrate-Inducible, GARP-type Transcriptional Repressor 1 (NIGT1, LOC_Os02g22020) and GA insensitive dwarf1 (GID1, LOC_Os03g57640), which were involved in nitrogen responses in rice through GA signaling [48,49]. NIGT1 overexpression affected some genes including GID1 [48], which overlapped with our DEGs under low N. Since NIGT1 is a negative regulator of NRT2.1 [50], the downregulation of NIGT1 under N stress and upregulation after the transfer to optimal nitrogen conditions was expected due to an increased availability of nitrate. This observation suggests the superiority of the Pokkali allele in enhancing N uptake under N stress. A recent study discovered the role of GID1 and nitrogen-mediated tiller growth response 5 (NGR5 (LOC_Os5g32270)) in gibberellin signaling by nitrogen fertilization to enhance tillering in rice [49]. NGR5, an APETALA2-domain transcription factor responsible for recruiting polycomb repressive complex 2 (PRC2) to regulate the transcription of tillering repressing genes, is dependent on nitrogen application. Semi-dwarf rice strains resist lodging under heavy N fertilization resulting in high yield despite the accumulation of growth repressing Della proteins whose degradation was regulated by GID1. Similarly, NGR5 is also targeted by GID1 for degradation. Although NGR5 was not among the Pokkali-specific DEGs, it was differentially expressed in other pairwise comparisons. Since there is a competition for the involvement of GID1 in both Della and NGR5 degradation, further investigation is needed to clarify these competitive interactions and how to modulate these interactions to enhance the grain yield of rice with reduced N fertilizer applications.
In our study, the alternative 3′ splice site (A3SS) was the major AS category which contrasts with intron retention reported in other studies in rice [51,52]. The number of AS events was more in the N deficient condition compared to the optimal N condition whereas no significant difference between two conditions was reported in an earlier study [19]. The 17 Pokkali-specific DEGs having AS events included transcription factors, stress-related genes, chloroplast precursor genes, and several key enzymes such as BIM2, zinc finger C3H4 type domain containing protein, pentatricopeptide (PPR) repeat containing protein, eukaryotic aspartyl protease domain containing protein, indole-3-glycerol phosphate synthase, isochorismate synthase 1, hydrolases, transferases, dehydrogenases (Table S6). All these results indicate the possible role of AS in regulating gene expressions for stress adaptation [53,54].
BIM2 (LOC_Os09g29930) is a bHLH type transcription factor that controls cell proliferation to cell lineage establishment [33]. Alternate splicing in this gene was only observed in Pokkali. It was downregulated and co-localized with root length QTL (rtlg; QTL ID: AQHE006) [55] and chlorate resistance QTL (kclo3rs; QTL ID: AQF102) [55], indicating its possible role in N uptake through the modulation of root architecture. Its involvement in brassinosteroid-regulated gene expression, that affects growth in Arabidopsis, further supports its role [56]. The late embryogenesis abundant protein, OsLEA3–2 (LOC_Os03g20680) is another Pokkali-specific DEG showing AS, which imparts tolerance to salt and osmotic stress [57]. It is highly upregulated in this study under N stress in Pokkali when compared to Bengal. Its co-localization to the QTL regions of NUE (qNUE-3) [58] and root number (rtnb; QTL ID: AQEL011) [55] suggests the importance of this gene during low N stress conditions.
Cinnamyl alcohol dehydrogenase 4 (LOC_Os11g40690), a key enzyme in the biosynthesis of an important cell wall component lignin [59], was downregulated in low N treatment in Pokkali root tissues, but its expression was higher than it was in Bengal. Pokkali’s superiority in developing an efficient root system architecture under low N conditions was further strengthened due to the co-localization of this gene with the QTLs for root number (QTL ID: CQAW29), root length (QTL ID: AQAR002), and chlorate resistance (QTL ID: AQF105) [55].
The LOC_Os02g09220 (cytochrome P450) showed the same expression pattern as cinnamyl alcohol dehydrogenase 4. It may be involved in maintaining plant growth under N stress due to its essential role in regulating growth and enhancing stress tolerance through the mediation of gibberellin homeostasis [60]. GDSL-like lipases, that are involved in biotic [61] and abiotic stress tolerance [62], were differentially expressed in our study as well as in a previous study [21]. Besides, Pokkali-specific DEGs included many disease resistance genes (LOC_Os08g28790; LOC_Os05g25350, LOC_Os10g04180, LOC_Os12g36860), which supports the notion of crosstalk as a common phenomenon for plants’ adaptation under different biotic and abiotic stress conditions including nitrogen stress.
4. Materials and Methods 4.1. Plant Materials and Cultivation
Plant materials for the preliminary screening of NUE included four japonica (Cocodrie, Bengal, Cheniere, and Cypress), three indica rice genotypes (PSRR-1, Nona Bokra, and Pokkali), and one aus genotype (Dular). All japonica genotypes were high-yielding cultivars released by LSU Agricultural Center. Dular is a lowland drought tolerant cultivar with deep roots [63]; Nona Bokra and Pokkali are salt tolerant land races [64,65]. PSRR-1 is a weedy rice genotype from Louisiana [66] and is closely related to indica subspecies [67].
4.2. Chlorate Assay
A chlorate assay was initially performed to screen for the nitrogen uptake capacity of eight genotypes in hydroponic experiments following Teng et al. [68] with two treatments: (a) control (no chlorate), and (b) chlorate-treatment (0.1% potassium chlorate solution). For plant materials, 100 seeds from each test cultivar were sanitized by a soaking in 5% hypochlorite for 30 min followed by rinsing with sterilized water. Seeds were germinated in petri-plates for 5 days at 28 °C. Ten seedlings with uniform growth per genotype were transferred into the hydroponics set-up; three replications and shoot lengths were measured after 7 days. A one-way ANOVA (alpha = 0.05) was performed to assess the differences in responses among cultivars and treatments.
4.3. Nitrogen Stress Response Experiment A nitrogen stress experiment was performed to measure the agronomic response of Pokkali (indica) and Bengal (japonica). It involved plants grown in a hydroponic set-up with a Yoshida nutrient solution in a semi-controlled environment. Two treatments were set up: (1) low nitrogen (LN), which contains one tenth of the normal amount of nitrogen (4 ppm of N) in Yoshida solution, and (2) full nitrogen (FN), which contains the normal amount of nitrogen (40 ppm of N). Seed sanitization and germination were similar to the chlorate assay mentioned earlier. There were three replicates for each treatment. Plant measurements were performed at 30 days after transferring to hydroponics. Five plants for each replicate were sampled for the measurement of shoot lengths, root-shoot fresh and dry weight, and chlorophyll content/leaf greenness. After taking the fresh weight, shoots and roots were oven-dried at 60 °C for one week to measure the dry weight. Chlorophyll measurements were performed using a SPAD meter (SPAD 502 Chlorophyll Meter, Spectrum Technologies, Inc., Aurora, IL, USA). The percent reduction in the shoot length, biomass, and chlorophyll content in low N compared to full N was calculated. All phenotypic data were subjected to a Student’s T-test (alpha = 0.05) to assess the difference in response to LN and FN between Pokkali and Bengal.
Intact root samples from at LN and FN treatment were collected and scanned using a Xerox®Altalink® C8035 (600 × 600 dpi) scanner. Image analyses were performed using the plug-in SmartRoot [69] in Image J [70]. Root architecture parameters such as the root length, root diameter/thickness, root number, root volume and surface area, were measured. A one-way ANOVA was conducted to identify significant differences between LN and FN treatments, and between Pokkali and Bengal.
4.4. Hydroponic Experiment for RNA-Seq Both parents were grown in a hydroponic setup under FN and LN, like the nitrogen stress experiment, and the seed preparation and cultivation conditions were also identical. Twenty seedlings were allocated for FN and 40 seedlings for LN. Root samples were collected at the 4-leaf stage around 14 days after transferring to hydroponics. Ten plants were collected and pooled for each replicate of the FN and LN treatment per cultivar. The remaining plants from the LN treatment were transferred into the FN treatment and then ten samples were collected and pooled for each replicate at a 1 h time point after transfer (recovery experiment). Collected root tissues were placed in liquid N during collection and stored in a −80 °C freezer until RNA extraction. 4.5. RNA Isolation, Library Preparation, and Sequencing The total RNA from root tissues of three biological replicates per treatment was isolated using RNApro™ Solution (MP Biomedicals, Santa Ana, CA, USA). The quality and quantity of RNA were evaluated using an ND-1000 Spectrophotometer (Thermofisher Scientific, Waltham, MA, USA) and were then treated with PerfeCTa DNase I (Quantabio, Beverly, MA, USA). Purified RNA samples were shipped to Novogene Co. Ltd. (Sacramento, CA, USA) for quality checking using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), library preparation, and sequencing using an Illumina HiSeq platform with a paired-end 150 bp sequencing strategy. 4.6. Reference Genome-Based Mapping
The quality of the raw reads was checked using FastQC (version 0.11.7, Babraham Bioinformatics, Cambridge, UK) with Phred quality score reads of ≥ 30 used in downstream analyses. Subsequent bioinformatics analyses were performed using high performance computing resources provided by Louisiana State University (http://www.hpc.lsu.edu). The processed paired-end reads were mapped to the Nipponbare reference genome (Os-Nipponbare-Reference-IRGSP 1.0) [71]) using HISAT2 (v2.0.1) [72]. All raw read sequences of Bengal and Pokkali were deposited in the NCBI’s sequence read archive (SRA) under the accession number PRJNA644654.
4.7. Transcript Assembly and Analysis of Alternative Splicing Events
The reference-based assembly of transcripts for each sample mapped with HISAT2 was done with Stringtie v.1.3.3 using the default settings [73]. The merging of transcript assemblies for each genotype per treatment was done with Stringtie using the merge option. The resulting GTF files were then processed in AStalavista v4.0 [74] for alternative splicing (AS) event identification. The analysis of these AS events was done following the method of Sammeth et al. [75], which identifies events such as intron retention (IR; 0,1^2-), exon skipping (ES; 0,1-2^), alternative 3′ splice sites (A3SS; 1-,2-), alternative 5′ splice sites (A5SS; 1^,2^), mutually exclusive exon (MXE; 1-2^,3-4^), and complex events (i.e., A5SS + A3SS, A5SS + ES + A3SS, etc.).
4.8. Differential Expression, Genotype and Treatment Comparison, and Ontology
Mapped reads were processed with featureCounts [76] to get the raw count table for genes in each sample. Genes with less than 10 reads across all samples were removed from downstream analyses. Log2 fold changes (LFC) were estimated using the DESeq2 R package [77] wherein genes with LFCs ≥ 2 (upregulated), ≤ 2 (downregulated), and a Benjamin–Hochberg adjusted p-value of <0.01 were considered as differentially expressed genes.
The expression profiles of Pokkali and Bengal in response to variable nitrogen treatments were analyzed: (1) to identify the DEGs responsible for each N treatment per genotype (N treatment effect) and (2) to identify the DEGs unique to the Pokkali per N condition (genotype effect). Analyses were performed for the following pair-wise comparisons: (1) Pokkali-low N (PKLN) vs. Bengal-low N (BGLN), (2) PKLN vs. Pokkali-full N (PKFN), (3) BGLN vs. Bengal-full N (BGFN), (4) Pokkali 1-h after transfer to full N (PK1H) vs. Bengal 1-h after transfer to full N (BG1H), (5) PK1H vs. PKLN, and (6) BG1H vs. BGLN. To identify the DEGs specific to various N treatments on each genotype, PKLN vs. PKFN and BGLN vs. BGFN were used for low nitrogen responses, and PK1H vs. PKLN and BG1H vs. BGLN for early N responses. To obtain the DEGs between Pokkali and Bengal under various LN and 1H treatments, PKLN vs. BGLN and PK1H vs. BG1H were used, respectively.
For gene ontology and plant pathway analyses, sets of DEGs were clustered in the following manner. For low N responses, DEGs from PKLN vs. PKFN, BGLN vs. BGFN, and PKLN vs. BGLN were combined. For early responses or 1-h N recovery periods, DEGs from PK1H vs. PKLN, BG1H vs. BGLN, and PK1H vs. BG1H were pooled. Gene ontology enrichment of the DEGs was performed to determine the biological significance with respect to the biological processes, molecular functions, and cellular localization of their proteins using a singular enrichment analysis (SEA) with agriGO v2 [78]. In addition, plant pathway analyses were performed using Plant Reactome (https://plantreactome.gramene.org) to determine the involvement of DEGs in plant pathways or processes [79].
4.9. Transcription Factors among DEGs
Transcription factors were identified and classified using the PlantTFDB v5.0, a plant transcription factor database [80]. The transcription factors among DEGs from the low N responses in Pokkali (PKLN vs. PKFN) and Bengal (BGLN vs. BGFN), and early N responses in Pokkali (PK1H vs. PKLN) and Bengal (BG1H vs. BGLN) were identified.
4.10. DEGs Co-Localized in Root Development-Related, ChloratE-Resistance, and NUE QTLs
DEGs identified as specific to N responses in Pokkali were co-localized with earlier reported QTLs for chlorate resistance, root development, and NUE. For potassium chlorate resistance (kclo3rs) and root development traits (root length: rtlg, root thickness: rtth, root volume: rtvol, root depth: rtdp and root number: rtnb), the Gramene QTL database [55] was used. For NUE-related traits, percent nitrogen content (qNCP-3-1), nitrogen use efficiency (qNUE-3) [58], and nitrogen response (qNR4 and qNR6) [81] were considered. In some cases, DEG co-localization with QTLs was not possible due to a lack of information on physical chromosomal positions (Table S7). Co-localized DEGs for each QTL were illustrated using MapChart 2.3 [82].
4.11. Validation of Expression of NUE-Related Genes: qRT-PCR
The qRT-PCR was done following the protocol used in our laboratory [83] to validate the gene expression patterns obtained from the RNA-seq analysis. A total of 1 μg DNase-treated RNA from pooled biological replicates was used to synthesize cDNA with an iScript cDNA Synthesis Kit (Bio-Rad Laboratories, Inc., Hercules, CA, USA). The elongation factor 1 alpha (EF1α, LOC_Os03g08010) was used as an internal control for expression normalization in different cDNA samples. The expression analysis was done for the following genes: WRKY55 (LOC_Os03g20550), transcription factor BIM2 (LOC_Os09g29930), ammonium transporter (LOC_Os02g40730), POT family protein/TGF-beta receptor (LOC_Os05g34000), NIGT1 (LOC_Os02g22020), polygalacturonase inhibitor precursor (PGIP1 (LOC_Os05g01380)), ZIP3 (LOC_Os04g52310), NRT2.1 (LOC_Os02g02170), high-affinity urea transporter DUR3 (LOC_Os10g42960), Pectinesterase inhibitor (LOC_Os12g18560), and MYB family transcription factor (LOC_Os05g50340). The reactions were set up in 96-well plates on an Applied Biosystems QuantStudio 3 Real-Time PCR system using iTaq™ Universal SYBR Green Supermix (Bio-Rad Laboratories, Inc) in a total reaction volume of 10 μL. The reactions were performed in three technical replicates using cDNA synthesized from pooled biological replicates and gene-specific primers designed using the PrimerQuest Tool (Integrated DNA Technologies, Inc., Coralville, IA, USA (Table S8)). The qRT-PCR expression levels of genes were determined using the 2–ΔΔCT method as described earlier [84]. The Pearson correlation coefficients were calculated on the expression data of genes obtained from both qRT-PCR and RNA-seq analysis using the R (v 4.0.1) package ggpubr with ggscatter function.
Enlarge this image.Figure 1. Effect of chlorate uptake on plant height in eight rice genotypes (Cocodrie, Bengal, Cheniere, Cypress, PSRR-1, Nona Bokra, Pokkali, and Dular). The mean seedling heights and standard errors were based on 10 control and 0.1% KCl03-treated seedlings per genotype. Groups are determined by a one-way analysis of variance relative to the respective control at p < 0.05. Means with similar letters indicate no significant difference.
Enlarge this image.Figure 2. Biomass (A), shoot length (B), and (C) chlorophyll content (SPAD readings) in Pokkali and Bengal under low nitrogen (4 ppm) and full nitrogen (40 ppm) in a hydroponic experiment. Means with similar letters indicate no significant difference at p < 0.05 (n = 15) determined using a one-way ANOVA.
Enlarge this image.Figure 3. Root architectures of Pokkali and Bengal seedlings under low and full nitrogen treatment.
Enlarge this image.Figure 4. Venn diagrams showing differentially expressed genes (DEGs) in Pokkali and Beng al under low N stress and 1-h recovery; |log2 fold change| ≥ 2 and padj < 0.01). (A) Comparison of low N vs. full N, and (B) comparison of low N vs. 1 h after transfer to full N. Treatment groups are PKFN (Pokkali-full nitrogen), PKLN (Pokkali-low nitrogen), PK1H (Pokkali 1 h after transfer from low N to full N), BGFN (Bengal-full nitrogen), BGLN (Bengal-low nitrogen), and BG1H (Bengal 1 h after transfer from low N to full N).
Enlarge this image.Figure 5. Number of differentially expressed genes in pairwise comparisons of low N and 1-h recovery treatments in Bengal and Pokkali. (A) On the stacked bars, blue represents downregulated genes while red represents upregulated genes. (B) Combined DEGs of all pairwise comparisons of low N and 1-h recovery treatments in both Bengal (BG) and Pokkali (PK). (C) Pokkali-specific DEGs for each treatment. The treatments included in the comparison are PKFN (Pokkali-full nitrogen), PKLN (Pokkali-low nitrogen), PK1H (Pokkali 1-h after transfer from low N to full N), BGFN (Bengal-full nitrogen), BGLN (Bengal-low nitrogen), and BG1H (Bengal 1-h after transfer from low N to full N).
Enlarge this image.Figure 6. The 148 Pokkali-specific root DEGs under low N and early recovery conditions. The expression pattern of these genes is shown for low nitrogen (LN) and 1-h N recovery/early response (ER) in Pokkali as well as compared to Bengal. The upward arrows in blue blocks and downward arrows in red blocks indicate upregulation and downregulation, respectively.
Enlarge this image.Figure 7. Gene ontology of DEGs from low N and 1-h recovery treatments in Pokkali and Bengal.
Enlarge this image.Figure 8. Number of DEGs identified related to phytohormone signaling during various nitrogen conditions. The plant hormone signaling genes identified were obtained from the combined data of low N (PKLN vs. PKFN, PKLN vs. BGLN and BGLN vs. BGFN) and 1-h response/early N response (PK1H vs. PKLN, PK1H vs. BG1H and BG1H vs. BGLN).
Enlarge this image.Figure 9. Pokkali-specific low N response DEGs co-localized on chromosomes with chlorate resistance, root development and nitrogen use efficiency (NUE) quantitative trait loci (QTL). Locus names in red font: downregulated, blue font: upregulated. QTLs included were chlorate resistance (kclo3rs): black solid bar; root number (rtnb): green diagonal bar; root volume (rtvol): red solid bar; root depth (rtdp): blue diagonal bar; root thickness (rtth): blue solid bar and root length (rtlg): red diagonal bar.
Enlarge this image.Figure 10. Annotation of alternate splicing (AS) events in the root transcriptome of rice. (A) Distribution of the various AS events in combined transcript assembly of all treatments in Pokkali (PK) and Bengal (BG). Others refer to complex events. (B) AS events are in the combined transcript assembly of each genotype and treatment combination.
| Genotype | Treatment | Average Lateral Root Length (cm) | Average Lateral Root Diameter (cm) | Total Lateral Roots (range) | Total Surface Area (cm2) | Total Volume (cm3) |
|---|---|---|---|---|---|---|
| Bengal | Full N | 8.20 ± 0.29 a,b | 0.029 ± 0.003 a | 5–8 | 4.19 ± 0.23 a | 0.036 ± 0.003 a |
| Low N | 9.09 ± 0.56 a | 0.032 ± 0.003 a | 7–9 | 7.40 ± 0.13 b | 0.089 ± 0.003 b | |
| Pokkali | Full N | 7.61 ± 0.40 b | 0.039 ± 0.001 b | 6–10 | 5.30 ± 0.33 c | 0.054 ± 0.006 c |
| Low N | 11.60 ± 0.46 c | 0.036 ± 0.001 a,b | 8–13 | 9.67 ± 0.45 d | 0.107 ± 0.006 d |
Means in each column followed by the same letter are not significantly different (p < 0.05).
| Major Pathway | Pathway Name | Low N (%) a | 1-h after Full N (%) b |
|---|---|---|---|
| Cellular process | DNA replication: activation of the pre-replicative complex | - | 1.4 |
| Protein metabolism: translation | 3.0 | 7.1 | |
| Growth and developmental process | Reproductive structure development | 5.0 | 5.7 |
| Vegetative structure development | 3.0 | - | |
| Metabolism and regulation | Amine and polyamine biosynthesis | 1.0 | - |
| Amino acid metabolism | 13.0 | 18.6 | |
| Carbohydrate metabolism | 8.0 | 5.7 | |
| Cofactor biosynthesis | 9.0 | 18.6 | |
| Cytokinin 7-N-glucoside biosynthesis | 2.0 | 1.4 | |
| Cytokinin 9-N-glucoside biosynthesis | 2.0 | 1.4 | |
| Detoxification | 2.0 | - | |
| Fatty acid and lipid metabolism | 2.0 | - | |
| Generation of precursor metabolites and energy | 2.0 | - | |
| Hormone signaling, transport, and metabolism | 25.0 | 21.4 | |
| Inorganic nutrients metabolism | 7.0 | 7.1 | |
| Photorespiration | 2.0 | - | |
| Secondary metabolism | 9.0 | 7.1 | |
| Response to stimuli: abiotic and stimuli and stresses | Response to cold temperature | 2.0 | 2.9 |
| Response to heavy metals | 1.0 | 1.4 | |
| Response to phosphate deficiency | 1.0 | - | |
| Response to stimuli: biotic and stimuli and stresses | Recognition of fungal and bacterial pathogens and immunity response | 1.0 | - |
Distribution shows percentage of differentially expressed genes (DEGs; |log2 fold change| ≥ 2.0 and padj < 0.01). a All low N DEGs (PKLN vs. BGLN, PKLN vs. PKFN, BGLN vs. BGFN); b all 1-h treatment DEGs (PK1H vs. BG1H, PK1H vs. PKLN, BG1H vs. BGLN).
| Function | MSU ID | Description | Low N | Early N Recovery | Low N | Early N Recovery |
|---|---|---|---|---|---|---|
| Pokkali | Bengal | |||||
| Nitrate transporters | LOC_Os10g40600 | Nitrate transporter 1.1B | - | U | - | U |
| LOC_Os02g02170 | High-affinity nitrate transporter 2.1 | U | U | - | U | |
| LOC_Os02g02190 | High-affinity nitrate transporter 2.2 | U | U | U | U | |
| LOC_Os01g50820 | High-affinity nitrate transporter 2.3 | U | - | U | - | |
| Ammonium transporter | LOC_Os02g40730 | Ammonium transporter 1 member 2 | D | U | - | U |
| Nitrate reductase | LOC_Os08g36480 | Nitrate reductase 1 | - | - | - | U |
| LOC_Os08g36500 | Nitrate reductase 2 | - | - | - | U | |
| LOC_Os02g53130 | NADH/NADPH-dependent NO3-reductase 2 | D | U | D | U | |
| Glutamine synthetase | LOC_Os03g12290 | Glutamine synthetase 1;2 | D | U | D | U |
| GOGAT | LOC_Os01g48960 | NADH-GOGAT | - | U | - | U |
U: upregulated; D: downregulated; ‘-‘: not differentially expressed.
| TF Family | Description | Low N | Early N Recovery | Low N | Early N Recovery |
|---|---|---|---|---|---|
| Pokkali | Bengal | ||||
| AP2 | AP2 family protein | 1 | 0 | 0 | 0 |
| B3 | B3 family protein | 0 | 0 | 3 | 1 |
| bHLH | basic/helix-loop-helix family proteins | 10 | 5 | 6 | 7 |
| bZIP | bZIP family protein | 3 | 6 | 0 | 6 |
| C2H2 | C2H2 zinc finger domain | 9 | 6 | 5 | 3 |
| C3H | Cys3His -containing zinc finger domain | 0 | 1 | 0 | 0 |
| CO-like | CO (CONSTANS) family protein | 1 | 1 | 0 | 1 |
| Dof | Dof (DNA binding with one finger) family | 2 | 0 | 0 | 0 |
| EIL | Ethylene-insensitive3 (EIN3) and EIN3-like (EIL) family proteins | 0 | 2 | 0 | 1 |
| ERF | ERF family protein | 9 | 18 | 8 | 15 |
| FAR1 | FAR1 family protein | 1 | 0 | 1 | 1 |
| G2-like | G2-like family protein | 5 | 3 | 1 | 0 |
| GATA | GATA family protein | 2 | 0 | 2 | 1 |
| GRAS | GRAS family protein | 1 | 0 | 1 | 1 |
| GRF | GROWTH-REGULATING FACTOR family protein | 0 | 1 | 0 | 1 |
| HD-ZIP | HD-ZIP family protein | 0 | 3 | 3 | 1 |
| HSF | Heat stress transcription factors | 1 | 1 | 0 | 1 |
| LBD | LBD family protein | 2 | 5 | 3 | 4 |
| MIKC-MADS | MIKC-MADS family protein | 1 | 0 | 1 | 0 |
| M-type MADS | M-type MADS family protein | 0 | 0 | 0 | 1 |
| MYB | MYB family protein | 8 | 1 | 4 | 4 |
| MYB-related | MYB-related family protein | 3 | 2 | 2 | 2 |
| NAC | NAM, ATAF, and CUC (NAC) transcription factors | 8 | 3 | 7 | 3 |
| NF-YA | NF-YA family protein | 2 | 0 | 2 | 0 |
| NF-YB | NF-YB family protein | 0 | 1 | 0 | 1 |
| Nin-like | Nin (for nodule inception)-like family protein | 1 | 1 | 1 | 0 |
| RAV | RAV family protein | 0 | 1 | 0 | 1 |
| SBP | SQUAMOSA promoter binding proteins (SBPs) | 1 | 0 | 1 | 0 |
| WOX | WOX family protein | 0 | 1 | 0 | 0 |
| WRKY | WRKY family protein | 17 | 2 | 9 | 9 |
| Total | 88 | 64 | 60 | 65 | |
Supplementary Materials
Supplementary Materials can be found at https://www.mdpi.com/1422-0067/21/16/5759/s1.
Author Contributions
Conceptualization, P.K.S.; methodology and data analysis, R.S.G., S.C., and R.T.; draft preparation, R.S.G. and P.K.S.; review and editing, P.K.S., R.S.G., S.C., and R.T. All authors have read and agreed to the published version of the manuscript.
Funding
This work was partly funded by the National Institute of Food and Agriculture, U.S. Department of Agriculture (Grant No. 2013-67013-21238) and the Louisiana Rice Research Board.
Acknowledgments
This manuscript was approved for publication by the Director of Louisiana Agricultural Experiment Station, USA as manuscript number 2020-306-34814.
Conflicts of Interest
The authors declare no conflict of interest.
Abbreviations
| N | Nitrogen |
| NUE | Nitrogen use efficiency |
| DEG | Differentially expressed genes |
| QTL | Quantitative trait loci |
| TF | Transcription factor |
| PK | Pokkali |
| BG | Bengal |
| PKFN | Pokkali-Full nitrogen |
| PKLN | Pokkali-low nitrogen |
| PK1H | Pokkali-1h after transfer from low N to full N |
| BGFN | Bengal-full nitrogen |
| BGLN | Bengal-low nitrogen |
| BG1H | Bengal-1h after transfer from low N to full N |
| GO | Gene ontology |
| SEA | Singular enrichment analysis |
| CEP | C-terminally encoded peptides |
| CLE | CLV3/ESR-related |
| rtnb | Root number |
| rtlg | Root length |
| rtth | Root thickness |
| rtvol | Root volume |
| rtdp | Root depth |
| kclo3rs | Potassium chlorate resistance |
| AS | Alternate splicing |
| IR | Intron retention |
| A3SS | Alternative 3' splice site |
| ES | Exon skipping |
| A5SS | Alternative 5' splice site |
| MXE | Mutually exclusive exon |
| GOGAT | Glutamine oxoglutarate aminotransferase |
| qRT-PCR | Quantitative reverse transcription PCR |
| NIGT1 | Nitrate-Inducible, GARP-type Transcriptional Repressor 1 |
| GID1 | GA insensitive dwarf1 |
| NGR5 | Nitrogen-mediated tiller growth response 5 |
| PRC2 | Polycomb repressive complex 2 |
1. Galloway, J.N.; Townsend, A.R.; Erisman, J.W.; Bekunda, M.; Cai, Z.; Freney, J.R.; Martinelli, L.A.; Seitzinger, S.P.; Sutton, M.A. Transformation of the nitrogen cycle: Recent trends, questions and potential solutions. Science 2008, 320, 889-892.
2. Boyle, E. Nitrogen pollution knows no bounds. Science 2017, 356, 700-701.
3. Golden, B.R.; Slaton, N.A.; Norman, R.J.; Wilson, C.E., Jr.; DeLong, R.E. Evaluation of polymer-coated urea for direct-seeded, delayed-flood rice production. Soil Sci. Soc. Am. J. 2009, 73, 375-383.
4. Good, A.G.; Shrawat, A.K.; Muench, D.G. Can less yield more? Is reducing nutrient input into the environment compatible with maintaining crop production? Trends Plant Sci. 2004, 9, 597-605.
5. Spiertz, J.H.J. Nitrogen, sustainable agriculture and food security. Agron. Sustain. Dev. 2010, 30, 43-55.
6. Vidal, E.A.; Alvarez, J.M.; Araus, V.; Riveras, E.; Brooks, M.; Krouk, G.; Ruffel, S.; Lejay, L.; Crawford, N.; Coruzzi, G.M.; et al. Nitrate 2020: Thirty years from transport to signaling networks. Plant Cell 2020.
7. International Rice Genome Sequencing Project (IRGSP). The map-based sequence of the rice genome. Nature 2005, 436, 793-800.
8. McCouch, S.; Zhao, K.; Wright, M.; Tung, C.W.; Ebana, K.; Thomson, M.; Reynolds, A.; Wang, D.; DeClerck, G.; Ali, L.; et al. Development of genome-wide SNP assays for rice. Breed. Sci. 2010, 60, 524-535.
9. Hu, B.; Wang, W.; Ou, S.; Tang, J.; Li, H.; Che, R.; Zhang, Z.; Chai, X.; Wang, H.; Wang, Y.; et al. Variation in NRT1.1B contributes to nitrate-use divergence between rice subspecies. Nat. Genet. 2015, 47, 834-838.
10. Ta, T.C.; Ohira, K. Effects of various environmental and medium conditions on the response of indica and japonica rice plants to ammonium and nitrate nitrogen. Soil Sci. Plant Nutr. 1981, 27, 347-355.
11. Ta, T.; Tsutsumi, M.; Kurihara, K. Comparative study on the response of indica and japonica rice plants to ammonium and nitrate nitrogen. Soil Sci. Plant Nutr. 1981, 27, 83-92.
12. Chao, D.Y.; Lin, H.X. Nitrogen-use efficiency: Transport solution in rice variations. Nat. Plants 2015, 1, 15096.
13. Good, A.G.; Johnson, S.J.; De Pauw, M.; Carroll, R.T.; Savidov, N. Engineering nitrogen use efficiency with alanine aminotransferase. Can. J. Bot. 2007, 85, 252-262.
14. Lea, P.J.; Azevedo, R.A. Nitrogen use efficiency. 2. Amino acid metabolism. Annal. Appl. Biol. 2007, 151, 269-275.
15. Pathak, R.R.; Ahmad, A.; Lochab, S.; Raghuram, N. Molecular physiology of plant nitrogen use efficiency and biotechnological options for its enhancement. Curr. Sci. 2008, 94, 1394-1403.
16. Shrawat, A.K.; Carroll, R.T.; DePauw, M.; Taylor, G.J.; Good, A.G. Genetic engineering of improved nitrogen use efficiency in rice by the tissue-specific expression of alanine aminotransferase. Plant Biotech. J. 2008, 6, 722-732.
17. Wang, Z.; Gerstein, M.; Snyder, M. RNA-Seq: A revolutionary tool for transcriptomics. Nat. Rev. Genet. 2009, 10, 57-63.
18. Yang, S.; Hao, D.; Song, Z.; Yang, G.; Wang, L.; Su, Y. RNA-Seq analysis of differentially expressed genes in rice under varied nitrogen supplies. Gene 2015, 555, 305-317.
19. Yang, W.; Yoon, J.; Choi, H.; Fan, Y.; Chen, R.; An, G. Transcriptome analysis of nitrogen-starvation-responsive genes in rice. BMC Plant Biol. 2015, 15, 31.
20. Xin, W.; Zhang, L.; Zhang, W.; Gao, J.; Yi, J.; Zhen, X.; Li, Z.; Zhao, Y.; Peng, C.; Zhao, C. An integrated analysis of the rice transcriptome and metabolome reveals differential regulation of carbon and nitrogen metabolism in response to nitrogen availability. Int. J. Mol. Sci. 2019, 20, 2349.
21. Sinha, S.K.; Sevanthi, V.A.M.; Chaudhary, S.; Tyagi, P.; Venkadesan, S.; Rani, M.; Mandal, P.K. Transcriptome analysis of two rice varieties contrasting for nitrogen use efficiency under chronic N starvation reveals differences in chloroplast and starch metabolism-related genes. Genes 2018, 9, 206.
22. Zhang, X.; Zhou, J.; Huang, N.; Mo, L.; Lv, M.; Gao, Y.; Chen, C.; Yin, S.; Ju, J.; Dong, G.; et al. Transcriptomic and co-expression network profiling of shoot apical meristem reveal contrasting response to nitrogen rate between indica and japonica rice subspecies. Int. J. Mol. Sci. 2019, 20, 5922.
23. Yamaya, T.; Obara, M.; Nakajima, H.; Sasaki, S.; Hayakawa, T.; Sato, T. Genetic manipulation and quantitative-trait loci mapping for nitrogen recycling in rice. J. Exp. Bot. 2002, 53, 917-925.
24. Quraishi, U.M.; Abrouk, M.; Murat, F.; Pont, C.; Foucrier, S.; Desmaizieres, G.; Confolent, C.; Rivière, N.; Charmet, G.; Paux, E.; et al. Cross-genome map based dissection of a nitrogen use efficiency ortho-meta QTL in bread wheat unravels concerted cereal genome evolution. Plant J. 2011, 65, 745-756.
25. Luo, B.; Tang, H.; Liu, H.; Shunzong, S.; Zhang, S.; Wu, L.; Liu, D.; Gao, S. Mining for low-nitrogen tolerance genes by integrating meta-analysis and large-scale gene expression data from maize. Euphytica 2015, 206, 117-131.
26. Liu, Q.; Chen, X.; Wu, K.; Fu, X. Nitrogen signaling and use efficiency in plants: What's new? Curr. Opin. Plant Biol. 2015, 27, 192-198.
27. Forde, B. Nitrogen signaling pathways shaping root system architecture: An update. Curr. Opin. Plant Biol. 2014, 21, 30-36.
28. Gao, Z.; Wang, Y.; Chen, G.; Zhang, A.; Yang, S.; Shang, L.; Ruan, B.; Liu, C.; Jiang, H.; Dong, G.; et al. The indica nitrate reductase gene OsNR2 allele enhances rice yield potential and nitrogen use efficiency. Nat. Commun. 2019, 10, 5207.
29. Manneh, B.; Stam, P.; Struik, P.C.; Bruce-Oliver, S.; van Eeuwijk, F.A. QTL-based analysis of genotype-by-environment interaction for grain yield of rice in stress and non-stress environments. Euphytica 2007, 156, 213-226.
30. Ogawa, S.; Valencia, M.O.; Ishitani, M.; Selvaraj, M.G. Root system architecture variation in response to different NH4+ concentrations and its association with nitrogen-deficient tolerance traits in rice. Acta Physiol. Plant. 2014, 36, 2361-2372.
31. Chrispeels, H.E.; Oettinger, H.; Janvier, N.; Tague, B.W. AtZFP1, encoding Arabidopsis thaliana C2H2 zinc-finger protein 1, is expressed downstream of photomorphogenic activation. Plant Mol. Biol. 2000, 42, 279-290.
32. Stracke, R.; Werber, M.; Weisshaar, B. The R2R3-MYB gene family in Arabidopsis thaliana. Curr. Opin. Plant Biol. 2001, 4, 447-456.
33. Toledo-Ortiz, G.; Huq, E.; Quail, P.H. The Arabidopsis basic/helix-loop-helix transcription factor family. Plant Cell 2003, 15, 1749-1770.
34. Nakano, T.; Suzuki, K.; Fujimura, T.; Shinshi, H. Genome-wide analysis of the ERF gene family in Arabidopsis and rice. Plant Physiol. 2006, 140, 411-432.
35. Fang, Y.; You, J.; Xie, K.; Xie, W.; Xiong, L. Systematic sequence analysis and identification of tissue-specific or stress-responsive genes of NAC transcription factor family in rice. Mol. Genet. Genom. 2008, 280, 547-563.
36. Rushton, P.J.; Somssich, I.E.; Ringler, P.; Shen, Q.J. WRKY transcription factors. Trends Plant Sci. 2010, 15, 247-258.
37. Krouk, G. Hormones and nitrate: A two-way connection. Plant Mol. Biol. 2016, 91, 599-606.
38. O'Brien, J.A.; Vega, A.; Bouguyon, E.; Coruzzi, G.; Gutierrez, R.A. Nitrate transport, sensing and responses in plants. Mol. Plant 2016, 9, 837-856.
39. Fredes, I.; Moreno, S.; Diaz, F.P.; Gutierrez, R.A. Nitrate signaling and the control of Arabidopsis growth and development. Curr. Opin. Plant Biol. 2019, 47, 112-118.
40. Sui, Z.; Wang, T.; Li, H.; Zhang, M.; Li, Y.; Xu, R.; Xing, G.; Ni, Z.; Xin, M. Overexpression of peptide-encoding OsCEP6.1 results in pleiotropic effects on growth in rice (O. sativa). Front. Plant Sci. 2016, 7, 228.
41. Tabata, R.; Sumida, K.; Yoshii, T.; Ohyama, K.; Shinohara, H.; Matsubayashi, Y. Perception of root-derived peptides by shoot LRR-RKs mediates systemic N-demand signaling. Science 2014, 346, 343-346.
42. Katayama, H.; Mori, M.; Kawamura, Y.; Tanaka, T.; Mori, M.; Hasegawa, H. Production and characterization of transgenic rice plants carrying a high-affinity nitrate transporter gene (OsNRT2.1). Breed. Sci. 2009, 59, 237-243.
43. Naz, M.; Luo, B.; Guo, X.; Li, B.; Chen, J.; Fan, X. Overexpression of nitrate transporter OsNRT2.1 enhances nitrate-dependent root elongation. Genes 2019, 10, 290.
44. Luo, B.; Chen, J.; Zhu, L.; Liu, S.; Li, B.; Lu, H.; Ye, G.; Xu, G.; Fan, X. Overexpression of a high-affinity nitrate transporter OsNRT2.1 increases yield and manganese accumulation in rice under alternating wet and dry condition. Front. Plant Sci. 2018, 9, 1192.
45. Chen, J.; Qi, T.; Hu, Z.; Fan, X.; Zhu, L.; Iqbal, M.F.; Yin, X.; Xu, G.; Fan, X. OsNAR2.1 positively regulates drought tolerance and grain yield under drought stress conditions in rice. Front. Plant Sci. 2019, 10, 197.
46. Wang, W.H.; Kohler, B.; Cao, F.Q.; Liu, G.W.; Gong, Y.Y.; Sheng, S.; Song, Q.C.; Cheng, X.Y.; Garnett, T.; Okamoto, M.; et al. Rice DUR3 mediates high-affinity urea transport and plays an effective role in improvement of urea acquisition and utilization when expressed in Arabidopsis. New Phytol. 2012, 193, 432-444.
47. Beier, M.P.; Fujita, T.; Sasaki, K.; Kanno, K.; Ohashi, M.; Tamura, W.; Konishi, N.; Saito, M.; Imagawa, F.; Ishiyama, K.; et al. The urea transporter DUR3 contributes to rice production under nitrogen-deficient and field conditions. Physiol. Plant. 2019, 167, 75-89.
48. Sawaki, N.; Tsujimoto, R.; Shigyo, M.; Konishi, M.; Toki, S.; Fujiwara, T.; Yanagisawa, S. A nitrate-inducible GARP family gene encodes an auto-repressible transcriptional repressor in rice. Plant Cell Physiol. 2013, 54, 506-517.
49. Wu, K.; Wang, S.; Song, W.; Zhang, J.; Wang, Y.; Liu, Q.; Yu, J.; Ye, Y.; Li, S.; Chen, J.; et al. Enhanced sustainable green revolution yield via nitrogen-responsive chromatin modulation in rice. Science 2020, 367, eaaz2046.
50. Maeda, Y.; Konishi, M.; Kiba, T.; Sakuraba, Y.; Sawaki, N. Molecular interactions of a soluble gibberellin receptor, GID1, with a rice DELLA protein, SLR1, and Gibberellin. Nat Commun. 2018, 9, 1376.
51. Shankar, R.; Bhattacharjee, A.; Jain, M. Transcriptome analysis in different rice cultivars provides novel insights into desiccation and salinity stress responses. Sci. Rep. 2016, 6, 23719.
52. Mansuri, R.M.; Shobbar, Z.S.; Jelodar, N.B.; Ghaffari, M.R.; Nematzadeh, G.A.; Asari, S. Dissecting molecular mechanisms underlying salt tolerance in rice: A comparative transcriptional profiling of the contrasting genotypes. Rice 2019, 12, 13.
53. Ashwini, N.; Sajeevan, R.S.; Udayakumar, M.; Nataraja, K.N. Identification of splice variant of OsGBF1 in Oryza sativa ssp. indica genotypes under salinity stress. 3 Biotech 2018, 8, 345.
54. Sellamuthu, G.; Jegadeeson, V.; Sajeevan, R.; Rajakani, R.; Parthasarathy, P.; Raju, K.; Shabala, L.; Chen, Z.; Zhou, M.; Sowdhamini, R.; et al. Distinct evolutionary origins of intron retention splicing events in NHX1 antiporter transcripts relate to sequence specific distinctions in Oryza Species. Front. Plant Sci. 2020, 11, 267.
55. Ni, J.; Pujar, A.; Youens-Clark, K.; Yap, I.; Jaiswal, P.; Tecle, I.; Tung, C.W.; Ren, L.; Spooner, W.; Wei, X.; et al. Gramene QTL database: Development, content and applications. Database (Oxford) 2009.
56. Yin, Y.; Vafeados, D.; Tao, Y.; Yoshida, S.; Asami, T.; Chory, J. A new class of transcription factors mediates brassinosteroid-regulated gene expression in Arabidopsis. Cell 2005, 20, 249-259.
57. Duan, J.; Cai, W. OsLEA3-2, an abiotic stress induced gene of rice plays a key role in salt and drought tolerance. PLoS ONE 2012, 7, e45117.
58. Senthilvel, S.; Vinod, K.K.; Malarvizhi, M.; Maheswaran, M. QTL and QTL × environment effects on agronomic and nitrogen acquisition traits in rice. J. Integr. Plant Biol. 2008, 50, 1108-1117.
59. Tobias, C.M.; Chow, E.K. Structure of the cinnamyl-alcohol dehydrogenase gene family in rice and promoter activity of a member associated with lignification. Planta 2005, 220, 678-688.
60. Zhou, J.; Li, Z.; Xiao, G.; Zhai, M.; Pan, X.; Huang, R.; Zhang, H. CYP71D8L is a key regulator involved in growth and stress responses by mediating gibberellin homeostasis in rice. J. Exp. Bot. 2020, 71, 1160-1170.
61. Lee, D.S.; Kim, B.K.; Kwon, S.J.; Jin, H.C.; Park, O.K. Arabidopsis GDSL lipase 2 plays a role in pathogen defense via negative regulation of auxin signaling. Biochem. Biophys. Res. Commun. 2009, 379, 1038-1042.
62. Naranjo, M.Á.; Forment, J.; Roldán, M.; Serrano, R.; Vicente, O. Overexpression of Arabidopsis thaliana LTL1, a salt-induced gene encoding a GDSL-motif lipase, increases salt tolerance in yeast and transgenic plants. Plant Cell Environ. 2006, 29, 1890-1900.
63. De Datta, S.K.; Chang, T.T.; Yoshida, S. Drought tolerance in upland rice. In Major Research in Upland Rice; International Rice Research Institute: Los Baños, Philippines, 1975; pp. 101-116.
64. Senadhira, D.; Neue, H.U.; Akbar, M. Development of improved donors for salinity tolerance in rice through somaclonal variation. SABRAO J. 1994, 26, 19-25.
65. Lin, H.X.; Zhu, M.Z.; Yano, M.; Gao, J.P.; Liang, Z.W.; Su, W.A.; Hu, H.; Ren, Z.H.; Chao, D.Y. QTLs for Na+ and K+ uptake of the shoots and roots controlling rice salt tolerance. Theor. Appl. Genet. 2004, 108, 253-260.
66. Subudhi, P.K.; Parco, A.; Singh, P.K.; Deleon, T.; Karan, R.; Biradar, H.; Cohn, M.A.; Brar, D.S.; Sasaki, T. Genetic architecture of seed dormancy in US weedy rice in different genetic backgrounds. Crop Sci. 2012, 52, 2564-2575.
67. Chai, C.; Shankar, R.; Jain, M.; Subudhi, P.K. Genome-wide discovery of DNA polymorphisms by whole genome sequencing differentiates weedy and cultivated rice. Sci. Rep. 2018, 8, 14218.
68. Teng, S.; Tian, C.; Chen, M.; Zeng, D.; Guo, L.; Zhu, L.; Han, B.; Qian, Q. QTLs and candidate genes for chlorate resistance in rice (Oryza sativa L). Euphytica 2006, 152, 141-148.
69. Lobet, G.; Pagès, L.; Draye, X. A novel image-analysis toolbox enabling quantitative analysis of root system architecture. Plant Physiol. 2011.
70. Schneider, C.A.; Rasband, W.S.; Eliceiri, K.W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 2012, 9, 671-675.
71. Kawahara, Y.; de la Bastide, M.; Hamilton, J.P.; Kanamori, H.; McCombie, W.R.; Ouyang, S.; Schwartz, D.C.; Tanaka, T.; Wu, J.; Zhou, S.; et al. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice 2013, 6, 4.
72. Kim, D.; Landmead, B.; Salzberg, S.L. HISAT: A fast-spliced aligner with 192low memory requirements. Nat. Methods 2015, 12, 357-360.
73. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290-295.
74. Foissac, S.; Sammeth, M. ASTALAVISTA: Dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucl. Acids Res. 2007, 35, W297.
75. Sammeth, M.; Foissac, S.; Guigo, R. A general definition and nomenclature for alternative splicing events. PLoS Comput. Biol. 2008, 4, e1000147.
76. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923-930.
77. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550.
78. Tian, T.; Liu, Y.; Yan, H.; You, Q.; Yi, X.; Du, Z.; Xu, W.; Su, Z. AgriGO v2.0: A GO analysis toolkit for the agricultural community, 2017 update. Nucl. Acids Res. 2017, 45, W122-W129.
79. Naithani, S.; Gupta, P.; Preece, J.; D'Eustachio, P.; Elser, J.; Garg, P.; Dikeman, D.A.; Kiff, J.; Cook, J.; Olson, A.; et al. Plant Reactome: A knowledgebase and resource for comparative pathway analysis. Nucl. Acids Res. 2019.
80. Jin, J.P.; Tian, F.; Yang, D.C.; Meng, Y.Q.; Kong, L.; Luo, J.C.; Gao, G. PlantTFDB 4.0: Toward a central hub for transcription factors and regulatory interactions in plants. Nucl. Acids Res. 2017, 45, D1040-D1045.
81. Wei, D.; Cui, K.; Ye, G.; Pan, J.; Xiang, J.; Huang, J.; Nie, L. QTL mapping for nitrogen-use efficiency and nitrogen-deficiency tolerance traits in rice. Plant Soil 2012, 359, 281-295.
82. Voorrips, R.E. MapChart: Software for the graphical presentation of linkage maps and QTLs. J. Hered. 2002, 93, 77-78.
83. Karan, R.; Subudhi, P.K. A stress inducible SUMO conjugating enzyme gene of a grass halophyte Spartina alterniflora (SaSce9) enhances salinity and drought stress tolerance in Arabidopsis. BMC Plant Biol. 2012, 12, 187.
84. Livak, K.; Schmittgen, T. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods 2001, 25, 402-408.
Prasanta K. Subudhi*, Richard S. Garcia, Sapphire Coronejo and Ronald Tapia†
School of Plant, Environmental, and Soil Sciences, Louisiana State University Agricultural Center, Baton Rouge, LA 70803, USA
*Author to whom correspondence should be addressed.
†Current address: Department of Horticultural Sciences, University of Florida, IFAS Gulf Coast Research and Education Center, Wimauma, FL 33598, USA.
© 2020. This work is licensed under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Garcia, Richard S; Sapphire Coronejo
; Tapia, Ronald