Transcriptome profiling analysis for two Tibetan wild barley genotypes in responses to low nitrogen

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₂O₂ 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. [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 (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 [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) [26]. A FDR (false discovery rate) was set as 0.05 for the threshold of DEGs [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 [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:/www.ncbi.nlm.nih.gov/tools/primer-blast/), and all the primers were presented in Additional file 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 [29] and the similar steps as reported by Zeng et al. [30]. The GOs distribution associated with DEGs were then obtained from three levels: biological process (BP), molecular function (MF) and cellular component (CC) [30]. The KEGG maps, which contained the EC numbers and enzymatic functions in many metabolic pathways were available in a variety of formats [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.

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

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 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 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 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 4: Table S3). These included both up-regulated (782) and down-regulated (728) genes (Fig. 1). Interestingly, DEGs in the tolerant genotype XZ149 (1203) were nearly twice as much as those in the sensitive genotype XZ56 (524) (Additional file 5: Table S4, Additional file 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. 1).

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. 2).

DEGs involved in nitrogen metabolism

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. 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. 3).

Kinases play important roles in the development of eukaryotic cells [32], and adaptation to abiotic stresses. Because some of their targets are transcription factors, they also have the functions of regulating transcription [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 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 [34, 35]. It has been proposed that some hormones coordinate demand and acquisition of nitrogen [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 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 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 5: Table S4, Additional file 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. 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. 4). The genes associated with catalytic activity and binding were the most enriched, accounting for 74.3 % of molecular function ontology (Fig. 4). Meanwhile, DEGs related to LN tolerance also acted as diverse cellular components (Fig. 4).

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. 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 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.

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 [39]. It was reported that high affinity nitrate transporters AtNRT2.1, AtNRT2.4 and AtNRT2.5 were induced in N-starved Arabidopsis roots [40, 41]. Similarly, the current study showed that the nitrate transporter DEGs were up-regulated under LN stress (Table 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 [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 [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 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 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 [23, 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 [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 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 [45]. An R2R3-type MYB TF, CmMYB1, enhances the expression of CmNRT, CmNAR, CmNIR, CmAMT, and CmGS under N-starvation [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 [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. 3). In addition, three MYB and three NAC transcription factors were found only up-regulated in XZ149 or only down-regulated in XZ56 (Fig. 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 [48, 49]. It has been well documented that CKT is involved in the repression of high-affinity NO₃ ⁻ transporter genes and root-shoot-root signaling in Arabidopsis [50, 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 8: Figure S2). Obviously, the regulation of CTK homeostasis in both accessions is consistent with the expression profile of high-affinity NO₃ ⁻ transporter genes.

Interestingly, we found that three ACO homologs were increased in XZ149 but remained little changed in XZ56 (Additional file 5: Table S4, Additional file 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 [45, 52]. Recent investigations on Arabidopsis detected a negative feedback loop between NRT2.1 expression and ethylene biosynthesis under low nitrate level [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 [54], due to the reduction in the frequency of electron carriers of electron transport systems [55]. Excessive accumulation of ROS in plant cells can result in oxidative stress, a major cause of cellular damage and cell death [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 [57]. In this study, three DEGs encoding putative GSTs were up-regulated only in XZ149 (Additional file 9: Table S7). In addition, Cytochrome P450s (CYPs), which are involved in biosynthesis and detoxification of a wide variety of molecules [58], had higher expression in rice [59] and sorghum [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 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 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 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 [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 [60]. Therefore, PAL has generally been considered as one of the main markers of environmental stress [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. 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 [62]. As an adaptive strategy, synthesis of photo-protective pigments, such as anthocyanins and flavonols could be enhanced in plants exposed to LN [63, 64], because anthocyanins may act as the ‘sunscreen’ and scavenger of ROS [62], and flavonols could filter off damaging wavelengths of radiation [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. 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. 6). Moreover, the DEGs (MLOC_68184) encoding CHI were down-regulated only in XZ56 and had little change in XZ149 (Fig. 6).

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 [66], and was also intersection of signaling pathways for auxin- and ethylene-mediated flavonol increases [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. 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.

Availability of data and materials

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