For comparative transcriptome analysis, we utilized the Tibetan hulless barley cultivar Zangqing 320, which shows a tolerance to N stress. The seeds were sterilized in 10% (v/v) H₂O₂ for 15 min, rinsed with distilled water, and germinated in a plant growth chamber under a 14 h day/10 h night cycle (28°C/21°C day/night temperature cycle) with 60–70% relative humidity. The uniform seedlings with three leaves were transferred into Hoagland solution. After 2 weeks, seedlings of uniform size and growth were picked randomly and planted into four treatment groups including optimal-nitrate (ON) (4 mM/L nitrate), free-nitrate (FN) (0 mM/L nitrate), low-nitrate (LN) (0.04 mM/L nitrate), and high-nitrate (HN) (40 mM/L nitrate) treatment groups. Ca(NO₃)₂·4H₂O and NH₄NO₃ were used as the N source and the consequent Ca²⁺ deficiency was supplemented with CaCl₂. The other components of the nutrient solution were as described previously (Xu et al., 2011). After being subjected to N stress for 2 days, seedling leaves were harvested and two independent replicates were collected for each sample. Plant materials were frozen in liquid nitrogen immediately and stored at −80°C until subsequent analyses.

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocols. Library construction was performed by staff at the Beijing Genome Institute (BGI, Shenzhen, China) comprising the following steps: enrichment of mRNA, fragment interruption, addition of adapters, size selection, and polymerase chain reaction (PCR) amplification. In total, eight paired-end libraries were constructed and 90 bp paired-end reads were generated using Illumina HiSeq™ 2000.

To ensure the high quality of sequencing data, clean reads with high quality were obtained by removing low-quality sequencing reads. Gene expression quantification was conducted using RSEM software (Li and Dewey, 2011). Fragments per kilobase pair of exon model per million fragments mapped (FPKM) was used to normalize gene expression values. The pairwise comparisons between N-treated and control samples were executed using the NOIseq method (Tarazona et al., 2012). DEGs were obtained based on a threshold of significance as P≥0.8. Cluster analysis of expression patterns was conducted using R language and Mev v4.7.4 software (Saeed et al., 2003). Gene ontology (GO) annotation and enrichment analyses were performed based on the GO Database (http://www.geneontology.org/) involving three hierarchies: biological process, molecular function, and cellular component. Pathway enrichment analysis of DEGs was performed utilizing the KEGG database (http://www.genome.jp/kegg/).

To validate the findings of the RNA-Seq assay, 20 DEGs were randomly chosen and their relative expression confirmed by qPCR analysis using the fluorescent intercalating dye SYBRGreen in the Opticon 2 detection system (MJ Research, Waltham, MA, USA). Details of the selected genes and the respective primers are listed in Table S1. Gene expression levels were normalized against the hulless barley gene HvADP (Ferdous et al., 2015) and calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). Three technical replicates were generated for each biological sample.

After removing sequencing adaptors and low quality data, we obtained 50,781,646 and 51,951,054 paired-end 90 bp reads from N treated and control samples, respectively, corresponding to approximately 4.57 and 4.68 Gb data (Table S2). The results of data quality assessment showed that clean reads exhibited good quality scores with Q20 percentages of all samples over 95%, whereas the uncalled base (“N”) percentages were lower than 0.01% (Table S2). The GC contents were almost identical for all eight hulless barley leaf tissues, ranging from 51.43 to 54.41%. These results indicated that the read number and quality were sufficient for further analysis.

On average, we mapped 77.61% of the clean reads to the full gene set of Tibetan hulless barley (Table S3). A total of 28,481 known expressed genes were detected in all samples, out of which 21,547 genes were expressed in all samples (Table S3). Gene expression data showed a Pearson's correlation between biological replicates of over 99.08% for all samples analyzed, indicating high correlation between biological replicates (Figure S1).

Compared with the ON group, 2428, 1274, and 1861 genes were classified as DEGs in the FN, LN, and HN groups, respectively (Table S4). Under FN stress, 837 genes showed increased expression and 1591 genes showed decreased expression (Table S4). In comparison, relatively fewer DEGs were identified under LN stress, wherein 675 transcripts were up-regulated and 599 were down-regulated (Table S4). In the HN group, 729 genes were induced and 1132 genes were suppressed (Table S4). Venn-diagram analysis revealed that 809 DEGs displayed differential expression under all three N stress treatments (Figure S2). In addition, 879, 138, and 407 DEGs were specifically detected in the FN, LN, and HN groups, indicating different transcriptional changes under varied N stress. To validate the accuracy and reproducibility of the RNA-Seq results, 20 DEGs were randomly chosen and the expression profiles were evaluated using qPCR. The qPCR and RNA-Seq analyses showed a positive correlation coefficient (R² > 0.80), suggesting the reliability of the RNA-Seq results (Figure S3).

Next, the three DEG sets were assigned to 51 GO classes by GO annotation analysis (Figure 1A). In the molecular function classification, “binding” and “catalytic activity” were dominant among the GO terms (Figure 1A). Under the classification of biological processes, “metabolic process,” “cellular process,” and “biological regulation” were prominently represented (Figure 1A). It is noteworthy that 303, 167, and 247 DEGs in the FN, LN, and HN groups, respectively, were separately annotated to “response to stimulus” (Figure 1A). The GO enrichment analyses further identified common and specific GO enrichment terms in the three DEG sets (Table S5). The DEGs in FN were significantly overrepresented in “photosystem I” (GO:0009522), “photosynthesis” (GO:0015979), and “nitrogen compound metabolic process” (GO:0006807). The DEGs in LN were found to be enriched in “positive regulation of catalytic activity” (GO:0043085), “cell killing” (GO:0001906), and “terpenoid biosynthetic process” (GO:0016114) (Table S5). In addition, the DEGs in HN were specifically enriched in “branched-chain amino acid metabolic process” (GO:0009081), “regulation of transport” (GO:0051049), and “oxidoreductase activity, acting on other nitrogenous compounds as donors” (GO:0016661) (Table S5).

KEGG pathway enrichment analysis also revealed that the three DEG sets were commonly enriched in nitrogen metabolism, N assimilation, and related metabolic pathways such as “biosynthesis of amino acids,” “carbon metabolism,” “carbon fixation in photosynthetic organisms,” “photosynthesis,” and “valine, leucine, and isoleucine degradation” (Figure 1B, Table S6).

Overall, the differences of enrichment degree and specific enrichment pathways implied that responsive differences existed between the three N stress groups. Thus, these results suggested that these pathways and processes might participate in N stress response.

To investigate the impact of stress degree differences on gene expression, we further compared the differences of DGEs between two pairwise comparisons (FN vs. ON and LN vs. ON). A total of 1010 DEGs overlapped, out of which 473 were upregulated and 368 were down-regulated (Figure S4). In addition, 254 DEGs appeared only in LN vs. ON comparison, while 1408 DEGs were in FN vs. ON comparison (Figure S4).

Hierarchical cluster analysis was then carried out based on expression profile of these common DEGs. On this basis, GO enrichment analysis of each cluster was performed and illustrated in Figure S5A. Enriched GO terms of gradually up-regulated genes under N deficiency included nitrogen compound metabolic process (GO:0006807), photosynthesis (GO:0015979), and carbohydrate catabolic process (GO:0016052). The common DEGs with down-regulated expression trend were significantly enriched in cellular amino acid catabolic process (GO:0009063), alpha-amino acid catabolic process (GO:1901606), and branched-chain amino acid catabolic process (GO:0009083). KEGG pathway enrichment revealed that numerous DEGs participated in biosynthesis of amino acids, photosynthesis, and carbon metabolism (Figure S5B). Among them, a high proportion of DEGs showed up-regulated expression profiles. The down-regulated DEGs mainly related to valine, leucine and isoleucine degradation, fatty acid degradation, and alanine, aspartate and glutamate metabolism.

We then focused on specific DEGs in FN vs. ON comparison, and identified 4 clusters based on dynamic expression profiles using a k-means clustering approach (Li et al., 2010). Cluster 1 and 3 displayed dramatically changes under free nitrogen stress (Figure 2A). The DEGs in cluster 1 significantly overrepresented in negative regulation of endopeptidase activity (GO:0010951), negative regulation of hydrolase activity (GO:0051346), negative regulation of protein maturation (GO:1903318,). Cluster 3 were highly enriched in response to stress (GO:0006950), defense response (GO:0006952), and regulation of gene expression (GO:0010468). KEGG pathway analysis further demonstrated distinct functional enrichments in biological process among four clusters (Figure 2B). Interestingly, there were 11 DEGs in cluster 4 significantly induced expressions under free nitrogen stress, involving in photosynthesis pathway. Carbon fixation in photosynthetic organisms and carbon metabolism also belonged to enrichment pathway for cluster 4.

To explore the genes associates with low- and high-N stress response in Tibetan hulless barley, we performed comprehensive comparisons of the three N stress groups. We identified 809 common DEGs among the three N stress-treated groups as compared to the ON group that were associated with multiple biological processes including “cellular aldehyde metabolic process,” “cellular amino acid metabolic process,” and “photosynthesis” (Figure 3). The common DEGs were further subjected to KEGG pathway enrichment analysis, which revealed that these common DEGs were significantly enriched in “nitrogen metabolism,” “carbon metabolism,” “photosynthesis,” “biosynthesis of amino acids,” and “starch and sucrose metabolism” (Figure S6). In contrast, the 407 specific DEGs detected in the HN vs. ON comparison were not significantly enriched in any biological processes (Figure S6). However, KEGG analysis found that these specific DEGs were primarily involved in “phenylalanine metabolism,” “phenylpropanoid biosynthesis,” “cyanoamino acid metabolism,” and “mRNA surveillance pathway” (Figure S6). These results suggested that Tibetan hulless barley had the ability to respond to N stress by regulating multiple common biological processes and pathways. In addition, a large number of specific DEGs involved in diverse biological processes were detected in the HN vs. ON comparison, indicating that the regulatory patterns of low- and high-N stress response also exhibited differences.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.