Phenotypic characterization of perennial ryegrass at different concentrations of nitrogen
To describe the growth state of perennial ryegrass under different treatments (Figure 1), the plant height, root length, fresh weight, tiller number and chlorophyll content of all samples were determined (Figure 2 A-D, Figure3). The results showed that the plant heights of treatments N0.5 and N1 were significantly higher than those of the other two treatments (Figure 2A). The fresh weight and tiller number of the N0.5 treatment were significantly higher than those of the other treatments (Figure 2C, 2D). However, the root length of treatment N0 as a control was highest and significantly higher than that in the other treatments (Figure 2B). The results of chlorophyll content determination indicated that there was a significant increase in treatment N0.5 compared to the other three groups; meanwhile, a trend of first up and then down was verified among these four treatments (Figure 3). Based on the above, it could be inferred that ryegrass in treatment N0.5 has a better growth state than the other three groups. Treatments N1 and N10 contained excessive nitrogen, which had a stress effect on plant growth.
Sequencing and transcriptome assembly
To comprehensively identify the transcripts of perennial ryegrass in response to nitrogen, 12 cDNA libraries were constructed from four treatments (N0, N0.5, N1, N10) corresponding to the conditions of no nitrogen, moderate amounts of nitrogen, a certain amount of excessive nitrogen and high concentration excessive nitrogen, each of which had three biological replications. In total, 80.58 GB clean data were released, and the content of each sample varied from 6.02 GB to 7.20 GB(table S1) All clean data were subjected to genomic localization analysis (http://126.96.36.199:5080/ryegrassgenome) by HISAT. A minimum of total mapped clean data occupied 61.84% in the relative example, and the proportion of multiple mapped reads or fragments in each example was less than 0.8% (table S2). Through transcriptome assembly, these transcripts were classified into 101376 unigenes, containing 2185 novel genes(table S3). The raw data was loaded on the SRA database in the NCBI(BioProject ID: PRJNA660099) (table S4).
Analysis of differentially expressed genes (DEGs)
To understand the response of the ryegrass gene to nitrogen, DEGs in three treatments (N0.5, N1 and N10) were screened with treatment N0 as a control. Under these conditions, 883, 1597, 1778 were identified separately as significantly differentially expressed in the treatments N0.5, N1 and N10. In treatment N0.5, 422 genes were upregulated and 461 genes were downregulated compared to the control group. Accordingly, the expression of 546 genes was upregulated and the expression of 1051 genes was downregulated in treatment N1. A total of 642 genes were upregulated and 1166 genes were downregulated in treatment N10 (table S5 A1-2, B1-2, C1-2, Figure 4A-C).
The differential expression of these DEGs was caused by the addition of nitrogen. To identify more genes related specifically to nitrogen utilization, three groups of DEGs were drawn into a Venn diagram (Figure 5a). The results showed that 345 genes are considered to be directly involved in the regulation of nitrogen utilization, among which 322 genes are known and the rest are novel genes(table S6).
In plant growth status identification, the conditions of treatment N0.5 were more suitable for ryegrass growth than the other three groups, and treatment N1 and N10 produced excessive nitrogen stress on the plants. Based on this, we further attempted to screen genes that specifically respond to excessive nitrogen stress. DEGs in treatments N1 and N10 were identified using treatment N0.5 as a control (table S5 D1-2, E1-2, Figure 4 D-E). There were 456 differentially expressed genes (169 upregulated, 287 downregulated) and 494 differentially expressed genes (154 upregulated, 340 downregulated) in treatments N1 and N10. Furthermore, 82 known genes and 22 novel genes were identified in the common genes in treatment N1 and N10 compared to treatment N0.5 (table S7, Figure 5b). These common genes may play an important role in the response of perennial ryegrass to excessive nitrogen stress.
Clustering of unigenes related to nitrogen utilization and excessive nitrogen stress
To gain insights into the expression pattern of those screened unigenes related to nitrogen utilization and excessive nitrogen stress, their expression levels in the different treatments containing different concentrations of nitrogen were clustered separately (Figure 6). Clustering analysis of unigenes related to nitrogen utilization revealed that the expressions of those genes in treatment N0.5 were most similar to those in treatment N1, followed by those in treatment N0, and then those in treatment N10, in general (Figure 6A). Considering the nitrogen content in different treatments, presumably, we could tell that these genes have analogous expression in response to moderate amounts of nitrogen and a certain amount of excessive nitrogen but distinct expression in response to high concentrations of excessive nitrogen. It is worth mentioning that most of the genes related to nitrogen utilization had relatively high expression in treatment N10, whereas they had relatively low expression in the other three groups. This result suggested that most of them possess similar expression patterns.
On the other hand, there were similar expressions of the genes related to nitrogen excessive nitrogen stress between treatment N0.5 and treatment N0, and something similar happened in treatment N1 and treatment N10 (Figure 6B). It could be inferred that those genes related to nitrogen excessive nitrogen stress are expressed homoplastically under the condition of no nitrogen and moderate amounts of nitrogen, in view of the different nitrogen contents in each treatment. Similarly, these genes had analogous expression in response to excessive nitrogen ranging from a certain concentration to high concentration.
GO and KEGG enrichment analysis of DEGs
To further understand the functional classification and metabolic pathways involved in the response of perennial ryegrass to nitrogen, GO and KEGG enrichment of DEGs was analyzed. We selected 30 GO terms with the most significant enrichment and showed them in the histogram (Figure 7). GO enrichment analysis of DEGs in treatment N0.5 compared to treatment N0 as a control revealed that all selected DEGs have significantly enriched GO terms, which are classified into “biological process” and “molecular function”, respectively. The “oxidation-reduction process” had the maximum enrichment in all “biological processes”. The number of DEGs enriched in “catalytic activity”, “transferase activity”, and “oxidoreductase activity” was also higher than the others in “molecular function” (Figure 7A). This result suggested that the genes of plant response to nitrogen were mainly enriched in “oxidation-reduction process”, “catalytic activity”, “transferase activity” and “oxidoreductase activity” under the conditions of relatively suitable nitrogen. Nevertheless, DEGs in treatment N1 compared to treatment N0.5 as a control only had five significantly enriched GO terms, which were mainly classified into “molecular function” (Figure 7B). We could deduce that more GO terms in “molecular function” are involved in the response mechanism of plants to nitrogen under a certain amount of nitrogen stress. Based on the GO enrichment analysis of DEGs in treatment N10 compared to treatment N0.5 (Figure 7C), there were no significantly enriched GO terms.
Furthermore, DEGs were separately enriched in 60, 37 and 38 pathways (Table S8, S9 and S10). Twenty pathway items with the most significant enrichment were selected for display in Figure 8. DEGs in treatments N0.5 and N1 (Figure 6A, 6B) displayed a certain degree of enrichment in the “nitrogen metabolism” pathway. The results in treatment N0.5 compared to N0 showed that the top enriched pathway of DEGs was “Photosynthesis-antenna proteins”, except “nitrogen metabolism” (Figure 8A). This result suggested that “photosynthesis-antenna proteins” may also respond positively to nitrogen under appropriate nitrogen conditions. Similarly, “steroid biosynthesis” and “carotenoid biosynthesis” were identified as the top two significantly enriched pathways in the KEGG enrichment analysis of DEGs between treatments N1 and N0.5 (Figure 8B). In addition, there were two relatively enriched pathways (“phenylpropanoid biosynthesis”, “phenylalanine metabolism”) with lower Q values (Figure 8B). This means that genes in these special pathways play an important role in the response mechanism of plants to a certain amount of nitrogen stress. Regarding the pathways enriched in treatment N10 compared to N0.5, “C5-branched dibasic acid metabolism” (only Novel00490|osa:4328147 enriched) and “steriod biosynthesis” pathways showed the top significant enrichment (Figure 8C). It can be indicated that “steriod biosynthesis” may perform some unknown and irreplaceable function under excessive nitrogen stress.
Photosynthesis-antenna proteins, steroid biosynthesis and carotenoid biosynthesis pathways
In KEGG enrichment analysis of DEGs, four special pathways with significantly enriched DEGs were found in three DEG groups caused by the different nitrogen contents. To further explore the response of these pathways to different nitrogen conditions, the enriched DEGs were mapped to those filtered pathways, and their expression was analyzed by clustering severely (Figure 9). There were 9, 2, 2, and 3 DEGs with significant enrichment classified into photosynthesis-antenna proteins, steroid biosynthesis and carotenoid biosynthesis pathways. These DEGs enriched in three special pathways displayed relatively broad expression patterns. Notably, the different DEGs in treatment N1 and N10 compared to treatment N0.5 were significantly enriched in the steroid biosynthesis pathway simultaneously. It could be indicated that the diverse parts of the steroid biosynthesis pathway function in the perennial ryegrass response to excessive nitrogen ranging from a certain concentration to a high concentration.
Nitrogen metabolism pathway
By KEGG enrichment analysis, we found 15 DEGs related to the right amount of nitrogen in ryegrass, which had homology to the previously identified genes in nitrogen metabolism pathway(table S8). 12 of them were successfully matched in the set of genes(table S6) we screened as genes related to nitrogen utilization. This was due to the different methods of statistical analysis of data. In the same way, 1 DEG and 0 DEG related to a certain amount of nitrogen stress and high concentration of nitrogen stress, respectively, performed homology to the previously identified genes in nitrogen metabolism pathway(table S9 and S10). This DEG was not included in the gene collection we screened as genes related to excessive nitrogen stress(table S7). Considering that we screened 345 and 104 DEGs as candidate genes related to nitrogen utilization and excessive nitrogen stress, respectively(table S6 and S7), it was reasonable to assume that the remaining 437 DEGs were genes related to the nitrogen response specifically in ryegrass, excluding the 12 DEGs enriched in the nitrogen metabolism pathway.
The expression analysis of these 12 DEGs showed that they performed a similar expression pattern, which was up-regulated in N0 treatment and down-regulated in other treatments, except
ref0006917-exonerate_est2genome-gene-0.0 with the reverse trend(Figure 10A). This indicated that ref0006917-exonerate_est2genome-gene-0.0 might play a positive role in nitrogen utilization, while the remaining genes performed some important functions in the absence of nitrogen.
Nitrate transporter family(NRT)
A total of 209 members of this transcriptome that were homologous to the NRTs in Arabidopsis were extracted, and 104 of them were differentially expressed in each treatment(table S11). The expression patterns of these differentially expressed genes varied, but they were more closely expressed in N0.5, N1, and N10 treatments than in N0 treatments(Figure 10B). This suggested that these NRTs might perform different functions in the absence of nitrogen than in the presence of nitrogen.
Identification of transcription factors
There is accumulating evidence that transcription factors (TFs) play key roles in various regulatory mechanisms(Hu, et al. 2017). Therefore, transcription factors from DEGs were identified for the purpose of further studying the molecular regulatory network of perennial ryegrass in response to nitrogen (Figure 11A). In total, 185 TFs, among which there were 14 novel genes, were distinguished in all differentially expressed genes. These TFs belonged to 37 transcription factor families, including the top three TF families (AP2-EREBP, MYB, NAC) with quantities greater than ten.
The expression patterns of these TFs changed in a relatively large range. Based on this, TFs related to nitrogen utilization and excessive nitrogen stress were characterized (Figure 11B, 11C). In total, 21 TFs related to nitrogen utilization, among which AP2-EREBP and MYB TF families occupied the largest number, were classified into 10 transcription factor families (Figure 10B). These TFs exhibited semblable expression in treatments N0.5, N1 and N10, and most of them possessed similar expression patterns in all treatments. Accordingly, 4 TFs related to excessive nitrogen stress were identified, which belonged to 4 transcription factor families, including LOB, NAC, AP2-EREBP and HB (Figure 11C). The expression of these TFs was relatively similar in treatments N1 and N10, and 2 TFs in the HB transcription factor family had similar expression patterns. Remarkably, ref0030344-processed-gene-0.2 from the AP2-EREBP transcription factor family concurrently exhibited a positive response to nitrogen absorption and excessive nitrogen stress.
Expression validation by qRT-PCR
To verify the accuracy of expression from RNA-Seq, 22 unigenes with significantly different expression levels were randomly selected for qRT-PCR analysis. These unigenes, among which 12 genes from DEGs related to nitrogen utilization and 10 from excessive nitrogen stress, contained 8 transcription factors and 14 functional genes. Almost all expression patterns of these genes were accordant in qRT-PCR analysis compared to those in the transcriptome profiling (Figure 12). Thus, the results of qRT-PCR with good coherence to those of RNA-Seq contributed to reflecting the high credibility and accuracy of transcriptome data.