Abstract

Background

   Nitrogen (N) is the most common limiting factor for crop productivity
   worldwide. An effective approach to solve N deficiency is to develop
   low N (LN) tolerant crop cultivars. Tibetan annual wild barley is
   well-known for its wide genetic diversity and high tolerance to poor
   soil fertility. Up to date, no study has been done to illustrate the
   mechanism of LN tolerance underlying the wild barley at transcriptional
   level.

Results

   In this study, we employed Illumina RNA-Sequencing to determine the
   genotypic difference in transcriptome profile using two Tibetan wild
   barley genotypes differing in LN tolerance (XZ149, tolerant and XZ56,
   sensitive). A total of 1469 differentially expressed genes (DEGs) were
   identified in the two genotypes at 6 h and 48 h after LN treatment.
   Genetic difference existed in DEGs between XZ149 and XZ56, including
   transporters, transcription factors (TFs), kinases, antioxidant stress
   and hormone signaling related genes. Meanwhile, 695 LN
   tolerance-associated DEGs were mainly mapped to amino acid metabolism,
   starch and sucrose metabolism and secondary metabolism, and involved in
   transporter activity, antioxidant activities, and other gene ontology
   (GO). XZ149 had a higher capability of N absorption and use efficiency
   under LN stress than XZ56. The higher expression of nitrate
   transporters and energy-saving assimilation pattern could be attributed
   to its more N uptake and higher LN tolerance. In addition, auxin (IAA)
   and ethylene (ETH) response pathways may be also related to the
   genotypic difference in LN tolerance.

Conclusion

   The responses of XZ149 and XZ56 to LN stress differed dramatically at
   transcriptional level. The identified candidate genes related to LN
   tolerance may provide new insights into comprehensive understanding of
   the genotypic difference in N utilization and LN tolerance.

Electronic supplementary material

   The online version of this article (doi:10.1186/s12870-016-0721-8)
   contains supplementary material, which is available to authorized
   users.

   Keywords: Barley, Low N tolerance, RNA-Seq, Genotypes, Differentially
   expressed genes

Background

   Nitrogen (N) is an essential mineral nutrient element for plant growth
   and development [[37]1], which serves as a constituent of many
   important macro-molecules, including proteins, enzymes, and several
   plant hormones [[38]2–[39]4]. On the whole, it is a major limiting
   factor for crop production in the world. Commonly, plants are often
   exposed to N deficient conditions with the situation becoming more
   severe due to increasingly-declined soil fertility and widely planted
   high-yield crop cultivars. Hence, a great amount of N fertilizer is
   applied to meet N requirement by crops [[40]3], which not only
   increases the cost for farmers but also brings the environmental
   problems. In addition it is well documented that only less than half of
   the applied N fertilizer is used by crops [[41]5], with the surplus
   contributing to severe environmental pollution. Therefore, it is
   extremely imperative and important to develop crop cultivars with high
   LN tolerance or N use efficiency (NUE), as it is a basic and also the
   most efficient approach for coping with low N availability in the soil
   and insufficient N fertilizer supply.

   Actually, plants have evolved many adaptive responses for coping with
   LN condition. Moreover, it was found that such N limitation
   adaptability in crops is closely associated with their yield
   performance [[42]6, [43]7]. Thus, genetic improvement of LN tolerance
   in crops would be of significance for developing sustainable
   agriculture. On the other hand, it has been well documented that NUE is
   a genetically controlled trait, differing dramatically among genotypes,
   such as in Arabidopsis [[44]8], as well as in crops including wheat,
   rice, maize and barley [[45]9–[46]12]. However, narrower genetic
   diversity in cultivated barley has become a bottleneck for genetic
   improvement [[47]13]. On the other hand, the Tibetan annual wild
   barley, growing in the Tibetan Plateau, which has been proved to be one
   of the centers of cultivated barley [[48]14], contains accessions with
   high tolerance to abiotic stresses, such as drought, salinity and
   potassium deficiency [[49]15–[50]17]. Meanwhile, the wild barley shows
   generally better adaption to poor soil fertility, and we also have
   identified some wild barley genotypes with high LN tolerance in a
   previous study [[51]18]. Therefore, it may be assumed that the wild
   barley contains some unique mechanisms of LN stress tolerance.

   RNA-Seq, one of next-generation high-throughput sequencing
   technologies, has been widely used recently, due to low background
   noise, high sensitivity and reproducibility, great dynamic range of
   expression and base-pair resolution for transcription profiling
   [[52]19]. Using this technique, transcriptomic profiles of many plants
   have been dissected under any given conditions, including biotic and
   abiotic stresses, such as heat stress [[53]20] and drought stress
   [[54]21]. Recently, this method has been also performed on some crops,
   viz rice, sorghum and cucumber [[55]22–[56]24], to investigate the
   mechanism of LN stress tolerance.

   In our previous research, we found a large genetic variation in LN
   tolerance among the wild barley accessions [[57]18]. However, up to
   date, no study has been done to analyze the genetic difference of
   transcriptomic profiles in response to LN stress. It is imperative for
   us to reveal the underlying mechanism or to explore the relevant genes
   of LN tolerance in wild barley. In this study, a comprehensive
   transcriptome analysis was conducted on two wild barley accessions
   (XZ149, LN-tolerant and XZ56, LN-sensitive), according to a previous
   study [[58]18]. The main objectives of the current study were to
   understand the mechanisms of LN tolerance existed in wild barley; and
   to determine the signaling pathways and regulatory networks of LN
   tolerance.

Methods

Plant materials and treatments

   The experiment was carried out in a greenhouse at Zijingang Campus,
   Zhejiang University, China. Healthy seeds of two wild barley genotypes
   (XZ149, LN-tolerant and XZ56, LN-sensitive) were disinfected with 2 %
   H[2]O[2] for 30 min, rinsed thoroughly with distilled water 5 times,
   and then soaked for 6 h at room temperature. Then the seeds were moved
   onto moist filter papers in germination boxes, which were placed in a
   growth chamber (22/18 °C, day/night) at dark for 3 days, and then
   incubated for another 7 days with light. At the second leaf stage
   (10-day old), the similar seedlings were transplanted into black
   plastic pots (5 L). Endosperm was removed away from the seedlings when
   they were transplanted, so as to eliminate nutrient supply from seeds.
   The used full-strength nutrient solution was the same as that reported
   by Yang et al. [[59]18]. The pH of the solution was adjusted to
   5.8 ± 0.1 with NaOH or HCl as required, and was continuously aerated
   with pumps and renewed every five days. Low N treatment was initiated
   on the three-leaf-stage seedlings, with 0.2 mM N as LN treatment and
   2 mM N as normal level (control).

   For biomass and N content determination, the seedlings were harvested
   and separated into shoots and roots, at 14 d after LN treatment. All
   the plant samples were heated at 105 °C for 30 min, dried at 80 °C
   until their weight remained constant, and then dry weight was recorded.
   N content was determined using Foss Kjeltec 8400.

   In order to know the time course of gene HvNRT2.1 expression under LN
   stress, the roots of XZ149 were sampled with three biological
   replicates at 3 h, 6 h, 12 h, 24 h, 48 h, 4 d, and 8 d after LN
   treatment, frozen in liquid nitrogen immediately, and stored at −80 °C
   for use in RNA extraction.

RNA-Seq sampling, RNA extraction and quality control

   For RNA-Seq sampling, the seeds of XZ149 and XZ56 were incubated at the
   same condition as mentioned above. The samples were taken at 6 h and
   48 h after exposure to LN stress (0.2 mM) and control (2 mM),
   respectively. Roots of four seedlings for each treatment were collected
   and mixed together at each time point. All samples [totally 16, 2
   genotypes (XZ149, LN tolerant and XZ56, LN sensitive) × 2 treatments
   (LN stress and control) × 2 time points (6 h and 48 h) × 2 biological
   replications] were prepared for further RNA-Seq analysis.

   Total RNA was extracted according to the instructions of miRNeasy mini
   kit (QIAGEN, Germany). RNA degradation and contamination were monitored
   on 1 % agarose gels. RNA purity was examined using the NanoPhotometer®
   spectrophotometer (IMPLEN, CA, USA). RNA concentration was determined
   using the Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life
   Technologies, CA, USA). RNA integrity was analyzed using the RNA 6000
   Pico Assay Kit of the Agilent 2100 Bioanalyzer (Agilent Tecnologies,
   Santa Clara, CA, USA).

RNA-Seq library construction, sequencing and reads mapping

   Sequencing libraries were generated using the Illumina TruSeq™ RNA
   Sample Preparation Kit (Illumina, San Diego, CA, USA) following the
   manufacturer`s instructions. Initially, mRNA was purified from the
   total RNA using poly-T oligo-attached magnetic beads. Then the purified
   mRNA was fragmented using divalent cations under elevated temperature
   in Illumina proprietary fragmentation buffer, and reversely transcribed
   into cDNA using SuperScript II (Invitrogen, Carlsbad, CA, USA).
   Illumina paired-end sequencing adapters were ligated for preparation of
   hybridization after making adenylation of the 3’ ends of DNA fragments.
   In order to preferentially select cDNA fragments of 150 bp, the library
   fragments were purified using an AMPure XP system (Beckman Coulter,
   Beverly, USA). With ligated adapters on both ends, DNA fragments were
   selectively amplified and enriched. Then PCR products were purified
   again using an AMPure XP system and quantified by Agilent Bioanalyzer
   2100 system. Finally, each final cDNA library was applied on one lane
   of the Illumina paired-end flow cell for the cluster generation
   processing, and subsequently sequenced on an Illumina NextSeq 500
   platform, thus generating 2 × 75 bp pair-ends reads.

   The raw reads were generated through the Illumina data processing
   pipeline (version 1.8). For further analysis, the clean data were
   obtained by removing low quality bases, empty reads and adaptor
   sequences at the 3’ end from the raw reads. Meanwhile, the Q20, Q30, GC
   contents, and sequence duplication level of the clean data were
   calculated. We downloaded the barley genome sequence and annotation
   data, and used the ultra high-throughput short read aligner to align
   RNA-Seq reads to the barley reference genomes on TopHat
   ([60]http://tophat.cbcb.umd.edu/), and then identified splice junctions
   between exons by analyzing the mapping results.

Identification of the DEGs and validation of RNA-Seq by quantitative RT-PCR

   For gene expression analysis, FPKM (fragments per kilo-base of exon per
   million fragments mapped reads) was calculated at expression level
   [[61]25]. The difference in expression between control and treatment
   (two biological replicates per time point) was analyzed using the DESeq
   R package (1.10.1) [[62]26]. A FDR (false discovery rate) was set as
   0.05 for the threshold of DEGs [[63]27].

   To confirm the validation of the RNA-Seq results, 1 μg total RNA of
   each sample for RNA-Seq was used for real-time quantitative PCR assays.
   After eliminating the genomic DNA contamination, first strand cDNA was
   synthesized with oligo dT primer and Random 6 mers in a 20 μl reaction
   (Takara, Japan). All quantitative RT-PCR (qRT-PCR) analyses were done
   in two biological replicates and three technical replicates using a
   CFX96 system (Bio-Rad, USA). The PCR profiles were as follows: 30 s at
   95 °C for pre-denaturation, 40 cycles of 5 s at 95 °C for denaturation
   and 30 s at 60 °C for annealing, followed by Melt-Curve analysis (60 °C
   - 95 °C, 0.5 °C increment for 5 s per step) to test the amplicon
   specificity. For relative quantification, the comparative CT method was
   used [[64]28]. The amplification of HvGAPDH sequence was used as an
   endogenous reference to normalize all the data. The gene-specific
   primers were designed using primer-blast
   (http:/[65]www.ncbi.nlm.nih.gov/tools/primer-blast/), and all the
   primers were presented in Additional file [66]1: Table S1.

Analysis of GO enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG)

   GO annotation and KEGG analysis for the DEGs were performed using the
   Blast2GO program [[67]29] and the similar steps as reported by Zeng et
   al. [[68]30]. The GOs distribution associated with DEGs were then
   obtained from three levels: biological process (BP), molecular function
   (MF) and cellular component (CC) [[69]30]. The KEGG maps, which
   contained the EC numbers and enzymatic functions in many metabolic
   pathways were available in a variety of formats [[70]31].

Statistical analysis

   Significant differences of physiological traits and gene expression
   among treatments and genotypes were tested using the Duncan’s Multiple
   Range Test (DMRT) on data processing system (DPS) statistical software,
   and the difference at P < 0.05 and P < 0.01 was considered as
   significant and highly significant, respectively.

Results

Effect of N level on growth performances of two wild barley accessions

   XZ149 and XZ56 were considered as LN tolerant and sensitive,
   respectively, according to a previous study [[71]18]. Although LN
   treatment caused a significant reduction of shoots dry weight for the
   two accessions, XZ149 was much less affected than XZ56, with XZ149 and
   XZ56 showing 7.0 % and 28 % reduction, respectively (Table [72]1).
   Meanwhile, LN enhanced root growth increase by 35 % and 28 % for XZ149
   and XZ56, respectively (Table [73]1). The less reduction of shoot dry
   weight and the much increase in root dry weight both contributed to the
   higher relative total plant dry weight in XZ149. On the other hand,
   although there was little difference in shoot N concentration between
   the two accessions under normal N, XZ149 had significantly higher shoot
   N concentration than XZ56 under LN (Table [74]1). As a result, N
   accumulation of XZ149 was 1.58 times larger than that of XZ56 under LN
   (Table [75]1). Obviously the current results proved that XZ149 is more
   LN tolerant than XZ56.

Table 1.

   Growth performances of two wild barley genotypes XZ149 (Low-N-tolerant)
   and XZ56 (Low-N-sensitive) at 14 d after low N treatment
   Trait XZ149 XZ56
   CK LN Relative CK LN Relative
   Dry weight (mg plant^−1 DW) Shoot 208.00a 192.83b 0.93 216.00a 134.00c
   0.62
   Root 52.17b 70.33a 1.35 56.58b 72.25a 1.28
   Total 260.17a 263.17a 1.01 272.58a 206.25b 0.76
   N concentration (%) Shoot 5.41a 4.67b 0.860 5.34a 4.27c 0.80
   N accumulation (mgplant^−1DW) Shoot 11.25a 9.01b 0.800 11.53a 5.72c
   0.50
   [76]Open in a new tab

   CK: Normal N level (2 mM N); LN: Low N level (0.2 mM N); Relative:
   LN/CK. For each line, different lowercase letters indicate significant
   differences (P, 0.05) among the treatments and genotypes, n = 6

Identification of DEGs

   In order to investigate an appropriate time of sampling for RNA-Seq
   analysis, we observed the response kinetic of the N starvation
   responsive gene, HvNRT2.1. Relative expression of HvNRT2.1 was
   calculated using pair-wise comparison between LN and control at 3 h,
   6 h, 12 h, 24 h, 48 h, 4 d, 6 d and 8 d after LN treatment,
   respectively (Additional file [77]2: Figure S1). The transcript level
   of HvNRT2.1 was slightly increased at 3 h and 6 h under LN stress,
   after then peaked abruptly at 12 h followed by another little peak at
   48 h, and thereafter decreased slowly but still remained on a little
   higher level at 4 d and 8 d (Additional file [78]2: Figure S1). The
   results demonstrated that the roots were capable of sensing LN signal
   and activating relevant signal transduction as early as 3 h after
   treatment, resulting in differential expression of the relevant genes,
   which in turn showed highly significant differences at 12 h and 48 h
   between the N levels. Accordingly, in view of sampling convenience in
   time arrangement, we then took the samples at 6 h and 48 h for RNA-Seq
   analysis.

   To obtain an overall view of the early LN responsive transcriptome in
   the two accessions, RNA samples were prepared from the roots of both
   accessions at 6 h and 48 h after LN treatment. Gene expression profiles
   of the wild barley roots under both control and LN conditions were
   analyzed. For each sample, two biological replicates were performed in
   sequencing. In total, 265597518 clean reads were obtained in the tested
   samples. For the most samples, 70 % of the sequenced reads could be
   uniquely mapped (Additional file [79]3: Table S2).

   The transcriptional levels were normalized using the FPKM method.
   Meanwhile, FDR < 0.05 was used as screening thresholds to test the
   significance of difference in transcript abundance. Consequently, 1469
   DEGs under LN stress were identified using pair-wise comparison of each
   accession between normal and low N conditions at each time point
   (Additional file [80]4: Table S3). These included both up-regulated
   (782) and down-regulated (728) genes (Fig. [81]1). Interestingly, DEGs
   in the tolerant genotype XZ149 (1203) were nearly twice as much as
   those in the sensitive genotype XZ56 (524) (Additional file [82]5:
   Table S4, Additional file [83]6: Table S5). However, there were 258
   DEGs commonly found in both XZ149 and XZ56. The two accessions
   displayed dissimilar expression patterns, in which the amount of
   up-regulated DEGs in XZ56 at 48 h decreased to one third of that at
   6 h, while XZ149 maintained little change between 6 h and 48 h
   (Fig. [84]1).

Fig. 1.

   Fig. 1
   [85]Open in a new tab

   A Venn diagram describing overlaps among differentially expressed genes
   (DEGs) in XZ149 and XZ56. a Up-regulated genes at 6 h and 48 h under
   low-N treatment. b Down-regulated genes at 6 h and 48 h under low-N
   treatment

   To confirm the validation of the RNA-Seq data, 15 responsive genes were
   randomly selected for quantitative RT-PCR analysis. The results from
   both qRT-PCR and RNA-Seq analysis showed that expressions of these
   genes were highly consistent, thus validating the RNA-Seq data
   (Fig. [86]2).

Fig. 2.

   Fig. 2
   [87]Open in a new tab

   Quantitative real-time PCR validation of 15 differentially expressed
   genes (DEGs). a Transcript levels of 15 DEGs and the corresponding
   expression data of RNA-Seq. The bars represent SE (n = 6). The columns
   represent relative expression obtained by qRT-PCR, and solid lines
   represent relative expression obtained by RNA-Seq. b Comparison between
   the relative expression obtained from RNA-Seq data and qRT-PCR. The
   RNA-Seq log[2] value of the relative expression (y-axis) has been
   plotted against the developmental stages (x-axis)

DEGs involved in nitrogen metabolism

   Many genes involved in nitrogen absorption and assimilation were
   differentially expressed under LN stress relative to the normal
   condition. In the current study, 12 DEGs encoding nitrate transporters
   were detected (Table [88]2). The abundance of these gene transcripts
   were increased under LN stress. Whereas, the expression patterns of
   these DEGs in XZ149 differed from those in XZ56, in which transcript
   levels of most DEGs were increased both in XZ149 and XZ56 at 6 h, and
   remained being enhanced in XZ149 but declined in XZ56 at 48 h. In
   addition, three DEGs (MLOC_60308, MLOC_56891, MLOC_51737) were
   up-regulated only in XZ149, and one DEG (MLOC_65110) was down-regulated
   only in XZ56. Two LN responsive genes encoding ammonium transporters
   (MLOC_33834, MLOC_35002) were identified (Table [89]2), in which,
   MLOC_35002 was responsive to LN stress only in XZ149. Moreover, seven,
   one and six DEGs encoding amino acid, lysine histidine and oligopeptide
   transporters, respectively, were found (Table [90]2). There were also
   some DEGs encoding the key enzymes in nitrate assimilation, including
   two nitrate reductases (NR), one nitrite reductase (NiR) and one
   glutamate synthase 1 (GOGAT 1) (Table [91]2). Interestingly, all of
   these DEGs were down-regulated in XZ149 but remained little changed in
   XZ56.

Table 2.

   Genes encoding protein transporters and nitrate assimilation enzymes
   showing genotypic difference expression in response to low N stress.
   Blank presented in the table means without significant difference in
   gene expression
   Group Gene id Log2(Fold change) Description
   XZ149 XZ56
   6 h 48 h 6 h 48 h
   Nitrate transporter MLOC_75087 1.23 1.19 1.51 High affinity nitrate
   transporter
   MLOC_73802 2.23 0.47 0.72 High-affinity nitrate transporter -like
   MLOC_3053 0.66 1.47 0.70 High-affinity nitrate transporter -like
   MLOC_60308 0.49 0.69 High-affinity nitrate transporter -like
   MLOC_14298 1.25 0.71 1.43 Nitrate transporter
   MLOC_1673 2.14 2.26 2.20 Nitrate transporter
   MLOC_58437 1.22 2.33 0.97 Nitrate transporter
   MLOC_52621 0.95 4.25 1.52 1.10 Nitrate transporter
   MLOC_58438 2.40 0.97 Nitrate transporter
   MLOC_56891 0.61 Peptide nitrate transporter
   MLOC_51737 0.85 Nitrate transporter
   MLOC_65110 −0.80 Nitrate transporter
   Ammonium transporter MLOC_33834 0.61 0.53 −0.56 Ammonium transporter
   MLOC_35002 0.58 Ammonium transporter 3 member 1-like
   Amino acid transporter MLOC_67247 0.57 Uncharacterized amino-acid
   permease
   MLOC_18832 0.97 Cationic amino acid transporter 5
   MLOC_62449 −0.62 Probable amino acid permease 7 isoform x1
   MLOC_29817 −0.74 Amino acid permease 6-like
   MLOC_36386 −0.69 −1.22 −0.67 −0.69 Amino acid-polyamine transporter
   MLOC_40001 −0.75 −0.98 Sodium-coupled neutral amino acid transporter
   2-like
   MLOC_47904 0.57 Amino acid permease 3-like
   Lysine histidine transporter MLOC_54046 0.75 Lysine histidine
   transporter 1-like
   Oligopeptide transporter MLOC_71333 2.10 −0.78 2.37 Oligopeptide
   transporter
   MLOC_76910 1.23 −0.68 1.02 Oligopeptide transporter 3-like
   MLOC_16637 0.54 0.76 Oligopeptide transporter 7
   MLOC_16638 1.14 Oligopeptide transporter 7
   MLOC_64771 0.63 0.59 Oligopeptide transporter 7-like
   MLOC_51375 1.04 1.02 Peptide transporter ptr2
   Phosphate transporter MLOC_12153 −1.18 −0.49 Phosphate transporter
   pho1-2
   MLOC_56639 −0.74 Phosphate transporter pho1-2
   MLOC_6492 −1.02 −1.03 Phosphate transporter pho1-3
   MLOC_28370 1.91 1.24 Phosphate transporter protein
   MLOC_48765 1.38 0.59 1.22 Phosphate transporter protein
   Potassium transporter MLOC_17989 0.75 High-affinity potassium
   transporter
   MLOC_74879 0.48 Potassium channel akt1-like
   MLOC_13576 −0.67 Potassium channel kor2-like
   MLOC_74565 −0.86 −0.60 Potassium transporter
   Sulfate transporter MLOC_61788 −1.84 −0.91 Probable sulfate transporter
   MLOC_44387 0.52 1.34 Sulphate transporter
   Zinc transporter MLOC_70898 0.86 Zinc transporter 4
   Iron transporter MLOC_57969 1.04 0.84 0.92 Fe(2+) transport protein
   1-like
   MLOC_16370 −0.69 Vacuolar iron transporter homolog 5-like
   MLOC_38253 −1.34 −0.56 Vacuolar iron transporter homolog 5-like
   MLOC_36909 −0.52 Vacuolar iron transporter -like
   MLOC_60633 1.39 1.46 Iron-phytosiderophore transporter ysl15-like
   MLOC_5701 0.61 0.48 Metal transporter nramp5-like
   MLOC_61170 0.63 0.91 Metal-nicotianamine transporter ysl3-like
   MLOC_67385 2.36 Metal-nicotianamine transporter ysl6
   MLOC_37227 0.50 Probable metal-nicotianamine transporter ysl12-like
   Molybdate transporter MLOC_61355 0.72 Molybdate transporter 1-like
   Nitrate reductase MLOC_5716 −1.11 −0.70 Nitrate reductase
   MLOC_60358 −0.69 Nitrate reductase
   Nitrite reductase MLOC_43860 −0.80 Ferredoxin--nitrite chloroplastic
   Glutamate synthase 1 MLOC_13604 −1.05 Glutamate synthase 1
   [92]Open in a new tab

   Meanwhile, expression of many genes associated with absorption or
   translocation of other nutrients changed under LN stress, such as
   phosphate (5), potassium (4), sulfate (2), Zinc (1), iron (9) and
   molybdenum (1), indicating that uptake of these nutrients in barley
   plants is affected by N metabolism under cross-talking regulation.

Transcription factors and protein kinases

   In total 89 DEGs encoding transcription factors (TFs) were identified
   in this study, and they belonged to different families, such as Zinc
   finger (30), bHLH (17), MYB (10), bZIP (10), ERF (7), NAC(7), WRKY (5)
   and HSF (3) (Fig. [93]3). XZ149 (75) had more than two fold TFs DEGs
   than XZ56 (34). Moreover, we found that the proteins with zinc finger
   domains were the most enriched among the TFs, accounting for 34 % of
   all these DEGs (Fig. [94]3).

Fig. 3.

   Fig. 3
   [95]Open in a new tab

   Average linkage hierarchical cluster analysis of transcription factors
   (TFs) identified in differentially expressed genes (DEGs). The samples
   and treatments are displayed above each column. Genes are displayed by
   different colors. Relative levels of expression are showed by a color
   gradient from low (blue) to high (red)

   Kinases play important roles in the development of eukaryotic cells
   [[96]32], and adaptation to abiotic stresses. Because some of their
   targets are transcription factors, they also have the functions of
   regulating transcription [[97]33]. In the current study, 85 DEGs
   encoding kinases of different groups were identified, which include
   serine threonine-protein kinase (STK), leucine-rich repeat
   receptor-like kinase (LRR), lectin protein kinase (LPK), cysteine-rich
   receptor-like kinase (CRK), CBL-interacting protein kinase (CIPK),
   wall-associated receptor kinase and calmodulin-binding receptor-like
   cytoplasmic kinase (CRCK) (Additional file [98]7: Table S6).

DEGs related to hormone signaling

   In addition to the basic roles in growth and development, phytohormones
   are also involved in various environmental responses, such as light,
   salt and drought [[99]34, [100]35]. It has been proposed that some
   hormones coordinate demand and acquisition of nitrogen [[101]36]. A
   total of 47 hormone signaling-related DEGs were found in this study,
   including gibberellin (GA, 6), cytokinin (CTK, 7), auxin (IAA, 18),
   ethylene (ETH, 7), abscisic acid (ABA, 8), jasmonic acid (JA, 7), and
   brassinosteroid (BR, 5) (Additional file [102]8: Figure S2). Heat map
   clustering analysis was performed to detect these DEGs involved in
   hormone signaling. It was found that cytokinin dehydrogenase DEG
   (MLOC_58639) was up-regulated and remained unchanged in XZ149 and XZ56
   under LN stress, respectively (Additional file [103]8: Figure S2). In
   addition, three aminocyclopropane-1-carboxylate (ACC) oxidase (ACO)
   DEGs were up-regulated in XZ149 but remained little changed in XZ56
   (Additional file [104]5: Table S4, Additional file [105]6: Table S5).

GO analysis and pathway enrichment analysis of LN tolerance related DEGs

   Of the total 1469 DEGs, 695 DEGs, showing significant up-regulation in
   XZ149, but down-regulation or unchanged in XZ56, or little change in
   XZ149 but down-regulation in XZ56, were selected for further
   investigation. Based on hierarchical clustering analysis, these DEGs
   could be mainly grouped into four classes. GO functional enrichment
   analysis were done to classify these DEGs into their corresponding
   biological process (BP), molecular function (MF) and cellular component
   (CC) (Fig. [106]4). The DEGs with known annotation could be categorized
   into 37 functional groups. In the biological process ontology, GO terms
   associated with ‘metabolic process’, ‘single-organism process’, and
   ‘cellular process’, occupied the majority (Fig. [107]4). The genes
   associated with catalytic activity and binding were the most enriched,
   accounting for 74.3 % of molecular function ontology (Fig. [108]4).
   Meanwhile, DEGs related to LN tolerance also acted as diverse cellular
   components (Fig. [109]4).

Fig. 4.

   Fig. 4
   [110]Open in a new tab

   Hierarchical cluster and gene ontology (GO) categories analysis of low
   N tolerance related DEGs. A total of 695 low N tolerance related 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). Gene ontology categories
   from three levels: (b) Molecular function; (c) Biological process; (d)
   Cellular component

   In addition to GO analysis, 695 DEGs were mapped to terms in KEGG
   pathway enrichment, and the encoded enzymes were assigned to 72 KEGG
   pathways (Fig. [111]5), including amino acid, nucleotide, lipid,
   carbohydrate, energy and other metabolisms. Among these pathways, DEGs
   involved in phenylpropanoid biosynthesis (34), phenylalanine metabolism
   (33), and starch and sucrose metabolism (13) were the most abundant.
   DEGs (MLOC_4686, MLOC_14829) (Additional file [112]4: Table S3)
   encoding the two key enzymes phenylalanine ammonia-lyase (PAL) and
   cinnamate 4-hydroxylase (C4H), respectively, involved in phenylalanine
   metabolism, as well as DEG (MLOC_6177) encoding flavonoid
   3'5'-hydroxylase (F3'5'H) homolog, were up-regulated only in XZ149 but
   not changed in XZ56 under LN stress. In addition, the DEG (MLOC_68184)
   encoding chalcone isomerse (CHI) was down-regulated only in XZ56 and
   had little change in XZ149. Ten enzymes encoded by 13 LN
   tolerance-related DEGs were found to be associated with starch and
   sucrose metabolism.

Fig. 5.

   Fig. 5
   [113]Open in a new tab

   KEGG overview of low-N tolerance related DEGs under low-N stress.
   X-axis represents the number of DEGs involving in each pathway; Y- axis
   depicts the different pathway

Discussion

   Nitrogen is an essential mineral nutrient required for plant growth and
   development. However, N deficiency in soil is becoming a big issue in
   crop production worldwide. It is imperative to develop the crop
   cultivars with high LN tolerance or NUE for coping with the issue.
   There is a large difference among species or genotypes within a species
   in LN tolerance. Relative biomass or dry weight is often used as an
   indicator of plant tolerance to low nutrition stress [[114]37,
   [115]38]. In the present study, the differences in growth performance
   between the two wild barley genotypes proved our previous result that
   XZ149 is more tolerant to LN stress than XZ56, and the results also
   indicated that XZ149 had the higher capability of N absorption and
   translocation than XZ56. Meanwhile, we used the RNA-Seq to reveal the
   differences of the transcriptome profiling between the two wild barley
   genotypes under low N stress. Clearly, there is a distinct difference
   in the transcriptional level between the two genotypes in their
   responses to LN.

Nitrogen metabolism genes responsive to LN stress

   Nitrate transporters are responsible for nitrate absorption from soils.
   On the whole, N deficiency enhances the expression of high affinity
   transport systems for nitrate [[116]39]. It was reported that high
   affinity nitrate transporters AtNRT2.1, AtNRT2.4 and AtNRT2.5 were
   induced in N-starved Arabidopsis roots [[117]40, [118]41]. Similarly,
   the current study showed that the nitrate transporter DEGs were
   up-regulated under LN stress (Table [119]2). Moreover, three nitrate
   transporter DEGs were up-regulated only in XZ149, and in general, the
   relative increase in abundance of most nitrate transporter DEG
   transcripts last longer in XZ149 than in XZ56. This unique and higher
   expression of nitrate transporters in XZ149 may contribute to higher
   nitrate uptake efficiency, producing more N-containing metabolites
   required for its survival under LN stress. Thus, it could be assumed
   that high LN stress tolerance observed in XZ149 is described to its
   more N uptake and accumulation in plants.

   Nitrogen assimilation is another fundamental biological process in
   plants, which is very energy consuming [[120]42]. The energy cost is
   particularly larger when nitrate is used as a major N source. When
   inorganic N is assimilated in roots, energy and C skeletons are
   provided through respiration, in which sucrose must be supplied from
   source leaves [[121]43]. Thus, the location of nitrate assimilation may
   affect energy budget in plants. Against the backdrop, it is reasonable
   that N assimilation related DEGs (Table [122]2) in the roots of the
   tolerant genotype XZ149 were depressed under LN condition, so as to
   keep energy supply more efficient than transporting sucrose from source
   leaves to roots for generating ATP and NAD(P)H used in nitrate
   assimilation. Alternatively, its nitrate and ammonium may be
   assimilated in the shoots after their transportation out of the roots
   via the xylem. This less energy-consuming counter measure may also
   contribute to XZ149’s tolerance to LN stress.

Carbon metabolism and reducing equivalents

   In order to survive under LN stress, some genes related to alleviation
   of the detrimental effect are abundantly expressed, thus resulting in
   enhanced stress tolerance. Currently, the transcription of soluble acid
   invertase (MLOC_60412) was increased in XZ149 but not in XZ56
   (Additional file [123]5: Table S4), indicating that sucrose degradation
   was enhanced in the tolerant genotype. The similar findings were
   observed in drought and LN stresses in other studies [[124]23,
   [125]44]. Therefore, we may assume that the enhanced invertase
   expression in the roots of the tolerant genotype XZ149 may stimulate
   the cycling of sucrose and carbon partitioning in favor of sucrose
   accumulation for counteracting the LN stress [[126]4].

   6-Phosphogluconate dehydrogenase (6PGDH) is a key enzyme involved in
   the pentose phosphate pathway, where it catalyzes the conversion of
   6-phosphogluconate to ribulose-5-phosphate, with the generation of
   NADPH. In the current study, the DEGs (MLOC_80338) encoding this enzyme
   were up-regulated in XZ149 (Additional file [127]5: Table S4), but
   unchanged in XZ56. As one of the primary end products of the pentose
   phosphate pathway, NADPH is necessary for fatty acid synthesis, and is
   needed in response to oxidative stress. Apart from that, it also helps
   to maintain the reduced state of glutathione (GSH) by serving as
   co-substrate for glutathione reductases that reduce oxidized
   glutathione.

Transcription factors

   Several TFs have been described in plants exposed to limited N. For
   example, Peng et al. found that NAC29 showed elevated expression at N
   deficiency [[128]45]. An R2R3-type MYB TF, CmMYB1, enhances the
   expression of CmNRT, CmNAR, CmNIR, CmAMT, and CmGS under N-starvation
   [[129]46]. It was reported that in fungi some GATA factors were
   involved in regulation of N metabolism and may activate expression of N
   catabolic enzymes during N-deficiency [[130]47]. However, little has
   been known about their functions in plants. Here, we found one GATA
   type zinc finger transcription factor family protein DEG (MLOC_53547)
   to be up-regulated only in XZ149 (Fig. [131]3). In addition, three MYB
   and three NAC transcription factors were found only up-regulated in
   XZ149 or only down-regulated in XZ56 (Fig. [132]3). Obviously, it will
   be quite imperative to determine the possible roles of these TFs in LN
   stress tolerance in the future.

Hormone signaling

   As a hormone, CTK (cytokinin) acts as a critical signalling molecule in
   communicating N availability in plants [[133]48, [134]49]. It has been
   well documented that CKT is involved in the repression of high-affinity
   NO[3]^− transporter genes and root-shoot-root signaling in Arabidopsis
   [[135]50, [136]51]. In this study, cytokinin dehydrogenase DEGs
   (MLOC_58639) involved in CTK degradation pathway was up-regulated and
   not changed in XZ149 and XZ56 under LN stress, respectively (Additional
   file [137]8: Figure S2). Obviously, the regulation of CTK homeostasis
   in both accessions is consistent with the expression profile of
   high-affinity NO[3]^− transporter genes.

   Interestingly, we found that three ACO homologs were increased in XZ149
   but remained little changed in XZ56 (Additional file [138]5: Table S4,
   Additional file [139]6: Table S5). It is well known that ACO, a
   rate-limiting enzyme in ethylene synthesis, can oxidize ACC (the direct
   precursor of ethylene synthesis) to ethylene,and ACO4 and another ACO
   homolog had responses to N deficiency [[140]45, [141]52]. Recent
   investigations on Arabidopsis detected a negative feedback loop between
   NRT2.1 expression and ethylene biosynthesis under low nitrate level
   [[142]53]. Thus we deduced that it may be the finely tuning between
   up-regulation of NRT2.1 expression and ethylene biosynthesis in XZ149
   should be beneficial for its LN tolerance.

Antioxidant stress

   Recent study visualized reactive oxygen species (ROS) accumulation when
   Arabidopsis roots were subjected to N deprivation [[143]54], due to the
   reduction in the frequency of electron carriers of electron transport
   systems [[144]55]. Excessive accumulation of ROS in plant cells can
   result in oxidative stress, a major cause of cellular damage and cell
   death [[145]56]. To cope with this, plants have developed numerous
   strategies for the detoxification of ROS. Glutathione S-transferases
   (GSTs) can directly scavenge peroxides with the help of glutathione as
   electron acceptor [[146]57]. In this study, three DEGs encoding
   putative GSTs were up-regulated only in XZ149 (Additional file [147]9:
   Table S7). In addition, Cytochrome P450s (CYPs), which are involved in
   biosynthesis and detoxification of a wide variety of molecules
   [[148]58], had higher expression in rice [[149]59] and sorghum
   [[150]23] under LN stress, indicating their specific effect of
   N-induced CYPs on the antioxidant level. Currently, we found nine CYPs
   DEGs were up-regulated in XZ149, but unchanged in XZ56, and another
   three unchanged in XZ149, but down-regulated in XZ56 under LN stress
   (Additional file [151]9: Table S7). The expression of one gene
   (MLOC_75745) for aldehyde dehydrogenase (ALDH) was up-regulated only in
   XZ149 under LN stress (Additional file [152]9: Table S7). Moreover, a
   large amount of peroxidases, up to 29, were up-regulated only in XZ149,
   apparently contributing to the higher capacity of antioxidant defense
   in XZ149 (Additional file [153]9: Table S7). Clearly, XZ149 may develop
   the higher ability of scavenging excessive ROS through forming a
   stronger antioxidant system under LN stress, facilitating higher
   tolerance to LN stress.

Phenylpropanoid metabolism pathway regulated by PAL under LN stress

   Phenylalanine is a key amino acid at the interphase of primary and
   secondary metabolism, and PAL is an initial rate limiting enzyme in
   phenylpropanoid synthesis [[154]60]. In addition to its important role
   in plant development, PAL also acts as a key enzyme in response to
   stress. Its biosynthesis would be stimulated in plants exposed to
   pathogenic attack, low temperature, salt stress, and low nitrogen,
   phosphate, or iron [[155]60]. Therefore, PAL has generally been
   considered as one of the main markers of environmental stress
   [[156]61]. Possibly, the response to N deficiency could be partly
   altered through regulation of PAL. In our study, DEG (MLOC_4684)
   encoding PAL was found only in XZ149 under LN, suggesting that the
   phenylpropanoid metabolism pathway mediated by PAL regulation may
   confer to genotypic difference in LN tolerance. Furthermore, lots of
   DEGs, up to 34, related to LN tolerance, were assigned to
   phenylpropanoid metabolism pathway according to KEGG metabolic pathway
   enrichment (Fig. [157]5). It may be assumed that the enhanced
   phenylpropanoid metabolism observed in XZ149 under LN stress may
   contribute to its high tolerance.

Regulatory network of the flavoniod synthesis pathway under LN stress

   LN stressed plants shows significant reduction of photosynthetic
   capacity and become more susceptible to oxidative damage caused by
   excessive light [[158]62]. As an adaptive strategy, synthesis of
   photo-protective pigments, such as anthocyanins and flavonols could be
   enhanced in plants exposed to LN [[159]63, [160]64], because
   anthocyanins may act as the ‘sunscreen’ and scavenger of ROS [[161]62],
   and flavonols could filter off damaging wavelengths of radiation
   [[162]65].

   Flavonols and anthocyanins synthesis shared the first two phases.
   Structural genes in anthocyanins and flavonols synthesis specially
   displayed higher transcriptional level under LN stress in the tolerant
   genotype XZ149 (Fig. [163]6). Starting from the first shared phase
   phenylpropanoid pathway, the DEGs encoding PAL, C4H and F3’5’H were
   up-regulated only in XZ149 and remained little changed in XZ56 under LN
   stress (Fig. [164]6). Moreover, the DEGs (MLOC_68184) encoding CHI were
   down-regulated only in XZ56 and had little change in XZ149
   (Fig. [165]6).

Fig. 6.

   Fig. 6
   [166]Open in a new tab

   The phenylpropanoid pathway for the biosynthesis of anthocyanin and
   flavonol. p-Coumaroyl CoA is at the junction of the metabolic routes
   leading to flavoniods synthesis. Genes are displayed by different
   colors. Relative levels of expression are showed by a color gradient
   from low (blue) to high (red). For each heatmap from left to right:
   XZ149-6 h (first column), XZ149-48 h (second column), XZ56-6 h (third
   column), XZ56-48 h (fourth column)

   In addition to the altered structure genes, the regulatory genes in
   anthocyanins and flavonols synthesis were also affected by LN stress. A
   MYB12 orthologous gene (MLOC_69570) was up-regulated in XZ149 under LN
   stress, while remain unchanged in XZ56. The MYB12 transcription factor,
   known to be a specific activator of flavonoid metabolism, activated the
   expression of chalcone synthase (CHS), CHI, flavanone 3-hydroxylase
   (F3H), and flavonol synthase (FLS) in Arabidopsis [[167]66], and was
   also intersection of signaling pathways for auxin- and
   ethylene-mediated flavonol increases [[168]67], indicating that the
   ethylene and auxin signalling may play an important role in
   low-N-induced flavonoid synthesis in barley roots. Moreover, three ACO
   homologs were up-regulated in XZ149 but unchanged in XZ56. Besides, two
   ethylene-responsive transcription factors (ERFs) (MLOC_18164,
   MLOC_75601) were unchanged/increased in XZ149 while decreased/unchanged
   in XZ56 (Fig. [169]3).

   In short, the current results strongly suggest that higher flavonoid
   accumulation in the tolerant genotype XZ149, especially the
   accumulation of anthocyanins and flavonols, may be attributed to its
   higher LN tolerance.

Conclusion

   Identification of DEG transcripts in plants under LN stress would
   reveal the genetic mechanism of LN tolerance. Here, the results of
   RNA-Seq analysis (all the sequences of the DEGs were available in
   Additional file [170]10: Table S8) demonstrated that there was the
   dramatic difference at transcriptional level between the two Tibetan
   wild barley genotypes in response to LN stress. Accordingly, a
   hypothetical model was developed for the mechanism of LN tolerance in
   XZ149 (Fig. [171]7). We deduced that more N absorption, less
   energy-consuming N assimilation pattern, more energy-producing model,
   which is contributed to high LN tolerance, may explain its well growth
   performance under LN stress. In addition, the current study identified
   some candidate genes related to LN tolerance, and cast a light on
   comprehensive understanding of the genotypic difference in N
   utilization and LN tolerance.

Fig. 7.

   Fig. 7
   [172]Open in a new tab

   A hypothetical model of low N tolerance mechanism underlying in XZ149.
   Genes are shown by different colors and relative expression levels are
   shown by a color gradient from low (blue) to high (red). For each
   heatmap from left to right: XZ149-6 h (first column), XZ149-48 h
   (second column), XZ56-6 h (third column), XZ56-48 h (fourth column)

Availability of data and materials

   The datasets supporting the conclusions of this article are included
   within the article and its additional files.

Acknowledgments