Transcriptomic analysis reveals adaptive strategies to chronic low nitrogen in Tibetan wild barley

Background Development of crop cultivars with high low nitrogen (LN) tolerance or nitrogen use efficiency (NUE) is imperative for sustainable agriculture development. Tibetan wild barley is rich in genetic diversity and may provide elite genes for LN tolerance improvement. Little has been known about transcriptional responses of the wild barley to chronic LN stress. Results In this study, transcriptomic profiling of two Tibetan wild barley genotypes, LN- tolerant XZ149 and LN-sensitive XZ56 has been conducted using RNA-Seq to reveal the genotypic difference in response to chronic LN stress. A total of 520 differentially expressed genes (DEGs) were identified in the two genotypes at 12 d after LN stress, and these DEGs could be mainly mapped to 49 metabolism pathways. Chronic LN stress lead to genotype-dependent responses, and the responsive pattern in favor of root growth and stress tolerance may be the possible mechanisms of the higher chronic LN tolerance in XZ149. Conclusion There was a distinct difference in transcriptional profiling between the two wild barley genotypes in response to chronic LN stress. The identified new candidate genes related to LN tolerance may cast a light on the development of cultivars with LN tolerance. Electronic supplementary material The online version of this article (10.1186/s12870-019-1668-3) contains supplementary material, which is available to authorized users.


Background
As an essential component of key macromolecules, nitrogen (N) is quantitatively the most important mineral nutrient for plants [1]. In soils, N is often the most important factor limiting plant growth, and plants frequently encounter N deficiency in their natural habitats. In the past several decades, the increasing use of N fertilizers in crop production has played a major role in increasing yields [2]. However, the main problem is the fact that crop plants only use less than half of the applied N [3], with the remaining N causing severe environmental pollution. There is thus an impending need to realize high productivity while decreasing the rate of N application. This necessitates a comprehensive understanding of molecular mechanisms underlying morphological and physiological adaptation to low N (LN) stress in crops.
Under LN stress, the sessile plants have evolved many adaptive responses. Plants improves the efficiency of N uptake by modifying root architecture, enlarging root system and enhancing the expression of high-affinity transport systems for nitrate and ammonium [4][5][6][7][8][9]. Meanwhile, N utilization efficiency could be also improved in the plants [10]. In addition, the remobilization of N from source organs might be stimulated when plants are subjected to N limitation [11], resulting in the enhancement of N re-assimilation, and maintenance of N economy in plants [12,13]. Furthermore, the expression levels of a number of genes associated with N metabolism in plants were altered to ensure the survival or complete their life cycle [14,15].
Tibetan annual wild barley is regarded as one of the progenitors of cultivated barley [16,17], possessing wider genetic diversity and generally better adaption to N deficiency in comparison with cultivated barley [18]. Some wild barley genotypes with high LN tolerance were identified, providing the elite genetic materials for improving LN tolerance of barley as well as other cereal crops [19]. Two wild barley genotypes differing dramatically in LN tolerance were used for transcriptome analysis at early stage of LN stress (6 h and 48 h) in our previous study [20]. However, plants respond to nutrient deficiency by inducing or repressing different sets of genes at special time [21]. Obviously, the mechanisms of LN tolerance in wild barley still remain to be revealed.
Although understanding of immediate responses is necessary for revealing molecular mechanisms of LN tolerance in plants, it still needs knowing the responses to long-term LN stress. In this study, we investigated genotypic differences in transcriptomic responses to chronic LN stress using the two wild barley genotypes (XZ149 and XZ56) differing in LN tolerance. The major objective of this study was to identify genes associated with the LN tolerance and to understand the molecular mechanisms underlying LN tolerance in wild barley.

Results
Effect of LN stress on growth performances of two wild barley genotypes XZ149 and XZ56 used in this study exhibited a distinct difference in plant growth at 12 d under LN stress (Fig. 1). Shoot dry weight of the two genotypes was significantly reduced under LN stress, with 15 and 37% reduction in comparison with control for XZ149 and XZ56, respectively (Fig.  1a). On the other hand, root dry weight was increased by 37 and 22% under LN stress for XZ149 and XZ56, respectively (Fig. 1b). As a result, total dry weight per plant under LN stress was only reduced by 7% in XZ149, but as much as 29% in XZ56 (Fig. 1c). LN stress also caused a significant reduction in SPAD value (chlorophyll content), N concentration and accumulation of the two genotypes, with XZ149 being reduction of 13, 11, 25% and XZ56 of 27, 19, 49%, respectively (Fig. 1d, e and f). The less reduction of shoot dry weight, N concentration and accumulation, and SPAD value, as well as the more increase in root dry weight in XZ149 than in XZ56 proved previous finding that XZ149 is more LN tolerant than XZ56 [18].

DEGs at the two N levels
The responding kinetics of HvNRT2.1 to LN stress was investigated in both XZ149 and XZ56. The relative expression was obtained using pair-wise comparison between normal and LN conditions at 1 d, 2 d, 4 d, 8 d and 12 d after LN treatment, respectively (Additional file 1: Figure S1). The transcript level of HvNRT2.1 was the maximum at 12 d after LN stress in both XZ149 and XZ56 (Additional file 1: Figure S1). In our previous study, 12 d was the time for sampling to determine the long-term changes of metabolites under LN stress [22]. Accordingly, in this study we took the samples at 12 d after LN stress for RNA-Seq analysis, to reveal the long-term responses to LN treatment.
In this study, eight cDNA libraries were constructed, and a total of 132,152,200 clean reads were obtained in the all samples. For most samples, 70% of the sequenced reads could be uniquely mapped (Additional file 2: Table S2). The results from real time PCR of 15 responsive genes highly validated the RNA-Seq data (Fig. 2).
The transcriptional levels were analyzed by calculating the FPKM. Meanwhile, FPKM ≥1 at least in one of the samples and FDR < 0.05 was set as screening thresholds for DEGs. Consequently, 520 DEGs were identified using pair-wise comparison for each accession (Additional file 3: Table S3, all the sequences of the DEGs were listed in Additional file 4: Table S4), including both up-regulated (229) and down-regulated (296) genes (Fig. 3). Notably, DEGs in XZ149 (278) were almost equal to those in XZ56 (290) (Fig. 3). Meanwhile, there were 48 DEGs shared commonly in XZ149 and XZ56 (Additional file 5: Table S5). The two accessions showed the different expression patterns, with XZ149 being basically equal in up-regulated and down-regulated DEGs and XZ56 having less up-regulated DEGs than down-regulated ones (Fig. 3).

DEGs over the different time under LN stress
Combined with the previous data on the samples taken in 6 h and 48 h after treatment, the expression patterns of different DEGs were clustered into 13 classes (Additional file 6: Figure S2). Across three time points examined, most of DGEs were only transiently down-or up-expressed at a specific time point in each genotype (Additional file 6: Figure S2). Only 29 DEGs in XZ149 and 11 in XZ56 were consistently expressed at three time points, with only one DEG being consistently up-regulated in XZ149, which encoded aminocyclopropane-1-carboxylate oxidase (ACO) (Additional file 7: Table S6). No gene was consistently down-regulated and only two DEGs, encoding nitrate transporter, were consistently up-regulated in both XZ149 and XZ56 at three time points (Additional file 7: Table S6). The results suggest that plants respond to LN stress by taking different countermeasures, inducing or repressing different sets of genes at certain time. In other words, the molecular mechanisms of LN tolerance could vary with growth stages.

Functional classification of DEGs
By hierarchical clustering analysis, 520 DEGs could be clarified into four classes (Fig. 4). Furthermore, these DEGs were enriched in 10 functional groups through GO enrichment analysis (Fig. 4b, Additional file 8: Figure S3). Molecular function Ontology included four GO terms 'catalytic activity' , 'antioxidant activity' , 'molecular function regulator' and 'binding' , while cellular component Ontology only consisted of term 'extracellular region'. The DEGs associated with 'metabolic process' were the most enriched and accounted for 35% of the biological process Ontology (Fig. 4b).
KEGG pathway analysis showed that 53 encoded enzymes by the 520 DEGs were mapped to 49 KEGG pathways (Fig.  4c), such as amino acid, carbohydrate, energy and secondary metabolite metabolisms. Among these pathways, 11 DEGs were involved in sucrose and starch metabolism, and three of them were up-regulated only in XZ149 and one down-regulated only in XZ56 (Additional file 9: Table S7).

DEGs involved in nitrogen/carbon metabolism and other nutrient uptake
Many genes involved in N absorption were differentially expressed, with ten DEGs encoding nitrate transporters, which were found at 12 d after LN stress ( Table 1). Most of these DEGs were up-regulated and their relative expression was higher in XZ149 than in XZ56 except MLOC_65110 and MLOC_70747 (Table 1). In addition, asparagine synthetase (ASN) DEG (MLOC_63089) involved in N assimilation was up-regulated in XZ149 but not changed in XZ56 (Table 1), and relative content of asparagine (Asn) was also higher in XZ149 than in XZ56 in both leaves and roots (Additional file 10: Table S8). Furthermore, in comparison with XZ56, XZ149 showed higher relative content of soluble protein and activities of the two key assimilating enzymes in leaves, namely NR and GS (Additional file 10: Table S8).
Carbon and N metabolisms are highly interconnected [23], and LN stress caused alteration in the expression of a b c d e f the genes related to carbon metabolism in barley roots.
In the current study, the DEG encoding soluble acid invertase (MLOC_60412) was up-regulated in XZ149, but not changed in XZ56 under LN stress (Table 1). Two trehalose-6-phosphate phosphatase (TPP) DEGs were identified, with one (MLOC_59991) being up-regulated in XZ149 and unchanged in XZ56, and the other one (MLOC_4115) being down-regulated in XZ56 and not changed in XZ149 (Table 1). Moreover, XZ149 had higher relative contents of sucrose and trehalose in comparison with XZ56 (Additional file 10: Table S8). Futhermore, the DEGs related to cell wall remodeling, including pectinesterase ( significantly changed under LN stress. MLOC_55102 showed elevated expression in XZ149 but not changed in XZ56, and MLOC_65390 and MLOC_59629 showed the reduced expression in XZ56 but not changed in XZ149 (Table 1). In addition to nitrate transporter, many DEGs encoding other transporters were also identified, including phosphate transporters, potassium transporters, amino acid transporters, sugar transporters and ABC transporters (Table 1), indicating that N metabolism could affect uptake of many other nutrients in barley through interactive regulation under long-term LN stress. It was also noted that all of the four DEGs encoding sugar transporters were up-regulated only in XZ149 or down-regulated only in XZ56 (Table 1).

Transcription factors and protein kinases
Totally 37 DEGs encoding transcription factors (TFs) were identified in XZ149 and XZ56 under long-term LN stress, less than half of the early responsive TFs (89).
These TFs belonged to different families, such as Zinc finger (8) (Fig. 5a). Among the 37 TFs, only one DEG was shared by the two genotypes. All other TFs showed different responses to LN stress between the two genotypes, while most DEGs were up-regulated in XZ149 and down-regulated in XZ56. It is noteworthy that only 13 TFs were found in response to long-term LN stress in XZ149, with nine TFs being up-regulated (Fig. 5a).
In the present study, 29 DEGs encoding different groups of protein kinases (PKs) were identified under LN stress (Fig. 5b). The PKs belonging to receptor-like kinase (RLK) family was the most enriched among the PKs. Six of the fifteen RLKs carried leucine-rich repeat (LRR) domains, and three contained lectin domains (Fig. 5b). Among RLKs, a wall associated kinase 3 homolog (MLOC_11098) was up-regulated in XZ149 but not changed in XZ56 (Fig. 5b). Besides, a N-acetyl-l-glutamate kinase (NAGK, MLOC_ 77291) was up-regulated in XZ149 and down-regulated in XZ56, while one NSP-interacting kinase 3-like DEG (MLOC_52843) was down-regulated in XZ56 but not changed in XZ149 (Fig. 5b).

Discussion
Nitrogen deficiency in soil is a major factor limiting crop production worldwide. It is imperative to understand the molecular mechanisms of LN tolerance in plants in order to improve LN tolerance through breeding. It is generally believed that the early responsive DEGs play more important roles for plant survival. Indeed, more DEGs were identified at 6 h (922) and 48 h (842) than at 12 d (520) after LN stress, suggesting the more genes are involved in early response to LN tolerance in plants. Moreover, plants respond to nutrient deficiency by regulating different sets of genes at special time [21]. Actually, the DEGs identified at 12 d of LN stress were different from those identified at the earlier time (Additional file 6: Figure S2). In the present study, we analyzed the transcriptomic responses of the two wild barley genotypes to LN stress by the RNA-Seq. Clearly, there is a distinct difference at the transcriptional level between the two genotypes in their responses to chronic LN stress.
The DEGs associated directly with N metabolism Nitrate metabolism starts from N acquisition from soil by nitrate transporters. N deficiency enhances the expression of nitrate transporters, and high affinity nitrate transporters played critical roles in N acquisition and remobilization in N-starved plants [5,6]. The higher relative expression of most nitrate transporter DEGs in XZ149 shows its higher ability of nitrate uptake and N reallocation compared with XZ56 under LN stress, which is also confirmed by higher relative soluble protein and activity of NR and GS in the leaves of XZ149. Consequently, XZ149 had higher N concentration and Fig. 4 Hierarchical cluster, gene ontology (GO) enrichment and KEGG analysis of 520 DEGs at 12 d after low N stress. A total of 520 DEGs were performed on a Hierarchical cluster analysis. The samples and treatments are displayed above each column. Genes are displayed by different colors and relative levels of expression are showed by a color gradient from low (blue) to high (red). b GO enrichment. c KEGG overview. X-axis depicts the different pathway; Y-axis represents the number of DEGs involving in each pathway more N accumulation in the plant tissues in comparison with XZ56 (Additional file 10: Table S8, Fig. 1e and f ).
It is well documented that N limitation acts as an exogenous senescence-triggering factor and enhances N remobilization [11]. Amino acids from protein degradation or assimilation will be translocated and reallocated as glutamine (Gln) and Asn, the major N carriers in the phloem of many plants [24][25][26]. Asn synthetase gene 1 (ASN1), controlling Asn synthesis will be up-expressed in the plants exposed to LN stress. In addition, over-expression of ASN1 in Arabidopsis enhanced its tolerance to limited N condition [27]. In the current study, an Arabidopsis ASN1 homolog DEG (MLOC_ 63089) was only up-regulated in XZ149, but not changed in XZ56 at 12 d under LN stress. In view of higher relative content of Asn in XZ149 than in XZ56 in both leaves and roots (Additional file 10: Table S8), it may be suggested that higher LN tolerance in XZ149 is partially contributed by its better reallocation of N.

Chronic LN stress alters carbon metabolism
It is well known that carbon and N metabolism are highly interconnected [23,28,29]. Indeed, 11 DEGs involved in sucrose and starch metabolism were found in barley roots under LN stress (Additional file 9: Table  S7). Moreover, four DEGs encoding sugar transporters were up-regulated in XZ149 but down-regulated in XZ56 (Table 1). Meanwhile, more accumulation of sucrose was also observed in the roots of LN-treated plants relative to control (Additional file 10: Table S8), and moreover, XZ149 accumulated more sucrose in its  roots than XZ56 under LN stress. As sucrose plays a role as a stimulant for lateral root formation [30], thus it could be understood that XZ149 showed better growth than XZ56 under LN stress (Fig. 1). The gene encoding soluble acid invertase (MLOC_60412) was up-regulated in XZ149 but not changed in XZ56 after 12 d under LN stress, and the results were consistent with those in the previous study, where the sampling were taken at 6 h and 48 h after LN treatment [20]. As invertase could irreversibly catalyze the hydrolysis of sucrose [31], it may be suggested that sucrose degradation was continuously enhanced in XZ149. In fact, it is well known that invertase is crucial in the carbon supply from sucrose to the non-photosynthetic cells of plants [32]. Thus, the consistently enhanced expression of invertase in the roots of XZ149 may provide more carbon source and energy, accounting for its better root growth under LN stress.
Trehalose-6-phosphate phosphatase (TPP) involves in trehalose biosynthesis, and trehalose plays a protective role against stress in plants [33]. Currently, one DEG encoding TPP was up-regulated in XZ149 but unchanged in XZ56, and another one was down-regulated in XZ56 but not changed in XZ149. The higher relative expression of DEGs encoding TPP in XZ149 may enhance its higher biosynthesis of stress protectant trehalose (Additional file 10: Table S8). It was reported that OsTPP1 was induced by stress [34], and its over-expression conferred stress tolerance in rice [35]. Therefore, it may be assumed that the higher relative expression of TPP and biosynthesis of trehalose in XZ149 may be a mechanism for protecting the plants against LN stress.
The DEGs associated with cell wall remodeling The dynamic and active structure of cell wall is necessary for plants to remodel when responding to external stimuli [36,37]. Cell wall would be remodeled when plants are subjected to nutrition deficiency, in particular N deficiency [36]. Cell wall related genes play a role in inducing cell wall relaxation and maintaining N assimilation under N starvation [38]. In this study, we found that cell wall remodeling DEGs, including PE, PEL, and PAE were more highly expressed in XZ149 than in XZ56 under LN stress. In addition, it should be noted that above-mentioned three enzymes are involved in pectin degradation. Hence, in comparison with XZ56, more pectin breakdown in XZ149 may provide it more materials for use in other biological processes, compensating for the depressed photosynthetic carbon assimilation under LN stress [15]. Thus, higher relative expression of cell wall related DEGs in XZ149 may be beneficial for its coping with LN stress.

Transcription factors and hormone signaling
Expression of plant-specific Dof1 TFs has been proved to improve plant growth under LN condition [39,40]. Several R2R3-type MYB TFs are involved in plant stress responses [41], and over-expression of OsMYB48-1 enhanced drought and salinity tolerance [42]. In this study, we found that one Dof and two MYB were down-regulated in XZ56, while three MYB were up-regulated in XZ149 (Fig. 5a). It may be assumed that these IFs may contribute to the development of complex signaling webs and should be potential candidates for improving LN stress tolerance.
The synthesis of plant hormone ethylene is tightly regulated in response to stresses [43]. Genes encoding ACO, a rate-limiting enzyme in ethylene synthesis, have been reported to be response to N deficiency [44,45]. In our previous study, three ACO homologs were up-regulated only in XZ149 under LN stress, however all of them were transiently expressed at either 6 h or 48 h [20]. Here, another ACO DEG (MLOC_34709) was detected, which was up-regulated in XZ149 at all the three time points, but unchanged in XZ56 at the later two time points (Additional file 7: Table S6). Considering the fine tuning between NRT2.1 expression and ethylene biosynthesis under LN level [46], it may be suggested that the up-regulation of MLOC_34709 in XZ149 contributes to its early and long-term LN tolerance.

Protein kinases
Protein kinases are important for adaptation to abiotic stress in plants [47]. The NSP-interacting kinase (NIK), belonging to the five leucine-rich repeats-containing receptor-like kinase subfamily, is a transducer of plant defense signaling [48]. In this study, we found that NSP-interacting kinase 3-like (MLOC_52843) was down-regulated only in XZ56 (Fig. 5b), suggesting that XZ149 may enhance its LN stress tolerance by maintaining its defense capability. In addition, a wall associated kinase 3 (WAK) homologous DEG (MLOC_11098) was up-regulated in XZ149 but unchanged in XZ56 (Fig. 5b). WAK is involved in root growth under LN stress [49], thus its higher expression in XZ149 could partially account for its better adaptation to LN stress.
N-acetyl-l-glutamate kinase (NAGK), the prototype of the amino acid kinase family, can catalyzes a key step of Arginine (Arg) synthesis. Arg content was increased in plants under N-deficient condition [50]. It was reported that some amino acids play roles in alleviating abiotic stress [51]. Currently, a NAGK homolog (MLOC_77291) was up-regulated in XZ149 and down-regulated in XZ56 (Fig. 5b), indicating that the increase of NAGK in XZ149 could confer its LN tolerance.

Conclusion
The current study identified the DEGs in barley roots, which respond to chronic (12 days) LN stress. There was a dramatic difference between the two Tibetan wild barley genotypes in transcriptomic response to chronic LN stress. The pattern of genetic responses in favor of better root growth and higher stress tolerance may be the possible mechanisms of the higher LN tolerance in XZ149. In addition, some new candidate genes related to LN tolerance were identified, which could be useful for developing barley cultivars with LN tolerance.

Plant materials and LN treatment
The experiment using two Tibetan wild barley accessions XZ149 and XZ56 (LN-tolerant and sensitive genotypes, respectively) was carried out in black plastic pots (5 L) in a greenhouse with natural light. The wild barley accessions were collected from Tibetan area in last century and kindly presented by professor Sun of Huazhong Agricultural University, China. The nutrition solution was used according to Quan et al. [20]. The solution was renewed every five days, continuously aerated with pumps. Treatments were conducted on three-leaf-stage seedlings with two N levels (0.2 mM N as LN treatment, 2 mM N as control).
For biomass determination, the plants were harvested at 12 d after LN treatment and separated into shoots and roots. Dry weight was recorded after the samples were dried at 105°C for 30 min and to constant weight at 80°C. Meanwhile fresh plant tissues were taken for use in determining nitrate reductase (NR) activity, glutamine synthetase (GS) activity and soluble protein content with three biological replications, and the content of metabolites with four biological replications. The roots of the two barley accessions under both N treatments were sampled with three biological replicates at 1 d, 2 d, 4 d, 8 d and 12 d after treatments for the time course analysis of gene HvNRT2.1 expression. For RNA-Seq analysis, the samples were taken at 12 d after treatments. Roots of four seedlings for each treatment were pooled as one biological replication. Totally eight samples [2 genotypes × 2 treatments × 2 biological replications] were taken for analysis.

Physiological measurement
N concentration in plant tissue was determined using Foss Kjeltec 8400. Soluble protein content was measured as described by Andrews et al. [52]. NR and GS activity were determined according to Kaiser et al. [53] and Masclaux-Daubresse et al. [54], respectively. The metabolites were extracted and analyzed using 7890A/5975C GC-MS system (Agilent, USA) and AMDIS 32 software according to Quan et al. [22].
cDNA library construction and sequencing for RNA-seq The total RNA was extracted using miRNeasy mini kit (QIAGEN, Germany) following the manufacturer's specification. RNA degradation, integrity, abundances and purity were assessed for meeting the requirements [20]. cDNA libraries were constructed using the Illumina Tru-Seq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Briefly, mRNA was obtained from the total RNA using magnetic beads with poly-T oligonucleotide. Then the random fragmentation of the purified mRNA was reversely transcribed into cDNA. After ligated with the adapters on both ends, DNA fragments were selectively amplified and enriched. Subsequently, the purified PCR products were quantified using Agilent Bioanalyzer 2100 system. After cluster generation, the final cDNA library was sequenced on an Illumina NextSeq 500 platform.
Raw reads with 75 bp single-end were initially processed to remove adapters sequences, empty sequences and low-quality bases, and then the Q20, Q30, GC contents, and sequence duplication level of the clean data were analyzed. Then the clean reads were mapped against the barley reference genomes using TopHat (http://tophat.cbcb.umd.edu/), and finally the mapping results were analyzed to identify splice junctions.

Identification of the DEGs and validation of RNA-Seq by real time PCR
Gene expression levels were calculated by the FPKM values (fragments per kilobase of exon per million fragments mapped reads) [55]. Fold-changes were defined as normalized read count abundance for the LN-stressed samples divided by that of the control samples. To identify differentially expression genes (DEGs), the difference in expression between control and LN treatment was analyzed using the DESeq R package (1.10.1) [56]. An FDR (false discovery rate) of 0.05 was used for determining significant DEGs [57].
To validate the reliability of the RNA-Seq results, the expression of candidate genes was determined by real time PCR assay using the RNA for RNA-Seq. The first strand cDNA was synthesized using PrimerScript™ RT reagent Kit with gDNA Eraser (Takara, Japan). The gene-specific primers, designed by primer-blast (http:/ www.ncbi.nlm.nih. gov/tools/primer-blast/), were presented in Additional file 11: Table S1. All real time PCR analyses were performed on a CFX96 system (Bio-Rad USA) with two biological replicates and three technical replicates. HvGAPDH was used as an internal control. The relative expression was calculated by the comparative CT method and expressed as the fold change referred to the expression in the control plants [58].