Transcriptomic and physiological analysis of common duckweed Lemna minor responses to NH4+ toxicity

Plants can suffer ammonium (NH4+) toxicity, particularly when NH4+ is supplied as the sole nitrogen source. However, our knowledge about the underlying mechanisms of NH4+ toxicity is still largely unknown. Lemna minor, a model duckweed species, can grow well in high NH4+ environment but to some extent can also suffer toxic effects. The transcriptomic and physiological analysis of L. minor responding to high NH4+ may provide us some interesting and useful information not only in toxic processes, but also in tolerance mechanisms. The L. minor cultured in the Hoagland solution were used as the control (NC), and in two NH4+ concentrations (NH4+ was the sole nitrogen source), 84 mg/L (A84) and 840 mg/L (A840) were used as stress treatments. The NH4+ toxicity could inhibit the growth of L. minor. Reactive oxygen species (ROS) and cell death were studied using stained fronds under toxic levels of NH4+. The malondialdehyde content and the activities of superoxide dismutase and peroxidase increased from NC to A840, rather than catalase and ascorbate peroxidase. A total of 6.62G nucleotides were generated from the three distinct libraries. A total of 14,207 differentially expressed genes (DEGs) among 70,728 unigenes were obtained. All the DEGs could be clustered into 7 profiles. Most DEGs were down-regulated under NH4+ toxicity. The genes required for lignin biosynthesis in phenylpropanoid biosynthesis pathway were up-regulated. ROS oxidative-related genes and programmed cell death (PCD)-related genes were also analyzed and indicated oxidative damage and PCD occurring under NH4+ toxicity. The first large transcriptome study in L. minor responses to NH4+ toxicity was reported in this work. NH4+ toxicity could induce ROS accumulation that causes oxidative damage and thus induce cell death in L. minor. The antioxidant enzyme system was activated under NH4+ toxicity for ROS scavenging. The phenylpropanoid pathway was stimulated under NH4+ toxicity. The increased lignin biosynthesis might play an important role in NH4+ toxicity resistance.


Background
Ammonium (NH 4 + ) and nitrate (NO 3 − ) are the two inorganic nitrogen (N) forms that can be directly absorbed by plants [1]. Compared to NO 3 − , NH 4 + is more easily absorbed by plants as its assimilation requires less energy. But in fact, only few plants are known to be NH 4 + specialists, most of high plants are usually sensitive to NH 4 + [2,3]. Non-specialists could display toxicity symptoms such as leaf chlorosis, growth suppression, yield depressions, and even mortality in high NH 4 + conditions, particularly when NH 4 + is supplied as the sole N source [3]. At the ecosystem level, some studies have even shown that increased NH 4 + in soil and water environment was associated with reduced crop yield, and decline of forest and macrophyte abundances [4][5][6].
NH 4 + toxicity is not only a significant ecological issue, but also an important plant physiological process [7]. Plant scientists have been trying to reveal its occurrence, signal transmission and physiological targets in plants [3,7,8]. Usually, for most plants, the root bears the brunt of NH 4 + toxicity [7,[9][10][11]. The root often accumulates high levels of NH 4 + in high NH 4 + condition, and the root cells could experience a futile transmembrane NH 4 + cycling that could carry high energetic cost resulting to decline plant growth [9,[12][13][14]. Associated enzymes involved in the regulation of NH 4 + influx, a signaling pathway model under NH 4 + toxicity in Arabidopsis thaliana has been described [7,8].
Additionally, recent studies also revealed that the NH 4 + toxicity could break the intracellular pH balance and C/N balance [3,7,15], and cause oxidative damage [16,17]. The knowledge on NH 4 + toxicity has greatly expanded in recent years, but the underlying mechanism are still largely unclear, further researches, especially in the subcellular level, using more advanced -omics approaches to follow up NH 4 + -induced global changes in plants are also required [8,18]. Transcriptome analysis is an effective method for global expression profiling of genes involved in stresses and toxicity in living organisms [19,20]. For example, transcriptomic profiling using microarrays have been used in Arabidopsis to identify molecular changes involved in NH 4 + toxicity [21]. With the rapid development of high-throughput sequencing, the next-generation transcriptome profiling approach or RNA sequencing (RNA-seq) has been gaining wide attention and use. RNA-seq could provide more information at a more affordable cost compared with the microarray and now an emerging powerful tool for transcriptome analysis [22].
Duckweeds are simple floating aquatic plants composed by frond and root. It has been considered to be a model species for aquatic plants and has been greatly used previously especially in the fields of toxicity studies, phytoremediation and biofuels production [23]. Lemna minor L. is one of the most widely distributed duckweed species and gains increasing interests due to its better adaptability to varying environmental conditions including high NH 4 + stress [24,25]. L. minor could grow well in high NH 4 + environment and has been even recognized as 'NH 4 + specialist' , but has been shown to still suffer toxicity in very high NH 4 + levels [15]. On the other hand, mechanisms and processes of toxicity in duckweeds however are a bit different from the terrestrial plant. Such as in Arabidopsis, most of the NH 4 + contact happens mainly in roots, thus the roots firstly suffer NH 4 + toxicity [7,26]. While for the floating duckweeds, the frond and root are all directly exposed to the toxic environment. This may lead to some different responses from the terrestrial plant. In this study, we use RNA-seq to investigate the global changes in common duckweed Lemna minor under NH 4 + toxicity. Those transcriptome analyses may provide some interesting insights and useful information not only in intoxication processes, but also on its potential tolerance mechanisms.

Sample preparation
Samples were prepared as described in Wang et al. [15]. L. minor was collected from a eutrophic pond in Chengdu, Sichuan, China (location: 30°38.86′N, 104°1 8.01′ E; elevation 499 m), and no specific permissions were required for specimen collection. To guarantee genetic uniformity, all of the L. minor materials originated from single colony and cultivated in Hoagland solution with 84 mg/L NO 3 − . The L. minor cultured in the Hoagland solution were used as the control (NC). For the treatments, cultures were grown in two NH 4 + concentrations, 84 mg/L (A84) and 840 mg/L (A840) in improved Hoagland solution, in which NH 4 Cl was used to provide NH 4 + , and KCl and CaCl 2 were used to replace KNO 3 and Ca(NO 3 ) 2 to avoid the impact of nitrate. All the solutions used in this study were adjusted to pH 5.5 with 1 M HCl.
Before inoculation, the fronds collected from Hoagland were washed five times with deionized water. Then, 0.2 g (fresh weight) of plant materials was cultivated in plastic basins with water depth of 2 cm. The plants were grown for one week in incubator at 23 ± 1°C with a photon flux density of 50-60 μmol · m −2 · s −1 provided by cool white fluorescent bulbs in a 16 h light/8 h dark cycle. The medium in each container was replaced every day.

Growth and physiological analysis
The relative growth rate (RGR) based on fronds number was used to evaluate the duckweed growth in different treatments as previously described in Wang et al. [15]. A total of 0.5 g fronds homogenized in 5 ml 0.1 % trichloroacetic acid was used for malondialdehyde (MDA) estimation by the thiobarbituric reaction following Dhindsa and Matowe [27]. Superoxide dismutase (SOD) was measured using a kit from Nanjing Jiancheng Bioengineering Institute (Jiangsu, China). Peroxidase (POD) and catalase (CAT) were measured by absorption photometry using a spectrophotometer as described by Bestwick et al. and Aebi [28,29], respectively. Ascorbate peroxidase (APX) activity assays were according to the method of Chen and Asada where the extinction coefficient of ascorbate at 290 nm was used for calculating APX enzyme activity [30].
Fronds of L. minor from the three treatments were stained by 3,3′-diaminobenzidine (DAB) or nitroblue tetrazolium (NBT) for measuring H 2 O 2 or O 2 − level, respectively [31]. Cell death was examined by Evans blue staining as described by Kim et al. [32].

RNA extraction, cDNA library preparation and sequencing
The whole plants with fronds and roots were ground in liquid nitrogen and total RNA was extracted using RNeasy® Plant Mini Kit (Qiagen) as per manufacturer's protocol. The integrity of RNA was assessed by formaldehyde agarose gel electrophoresis. A total of 30 μg mixed RNA from three biological replicates detected by 2100 Bioanalyzer (Agilent, USA) was digested with DNase I (TAKARA), and then purified by Dynabeads® Oligo (dT)25 (Life, USA). 100 ng derived mRNAs were fragmented and reverse transcribed into first-strand cDNAs with random hexamer and then the secondstrand cDNAs were synthesized by using a NEBNext® Ultra™ RNA Library Prep Kit for Illumina (NEB). The double-stranded cDNAs were purified and ligated to adaptors for Illumina paired-end sequencing. The cDNA library was sequenced using the Illumina HiSeq2500 system by Shanghai Hanyu Biotech lab (Shanghai, China).

De novo assembly of RNA-seq reads and quantifying gene expression
For the assembly library, raw reads were filtered using the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) to remove adapters and low-quality reads (base quality < 20, read length < 40 bp). The obtained quality-filtered reads were de novo assembled into contigs by the Trinity Program [33]. Unigenes were defined after removing redundancy and short contigs from the assembly. The unigenes were predicted by "GetORF" in the EMBOSS package [34] and aligned to the protein sequence database NCBI NR (non-redundant protein database), Swiss-Prot (Annotated protein sequence database), KEGG (Kyoto encyclopedia of genes and genomes) and COG (Clusters of orthologous groups of protein) by Blastp with an E-value threshold of 1 × 10 −5 .
The number of unique-match reads was calculated and normalized to RPKM (reads per kb per million reads) for gene expression analysis. Comparison of unigene expression between treatments was according to DESeq as described by Abders and Huber [35]. The differentially expressed genes (DEGs) between NC and A84, or between NC and A840, or between A84 and A840 were restricted with FDR (false discovery rate) ≤ 0.001 and the absolute value of log2 Ratio ≥1.
To examine the expression profile of DEGs, the expression data υ (from NC, A84 and A840 treatment) were normalized to 0, log2 (υ A84 /υ NC ), log2 (υ A840 /υ NC ), and then clustered by Short Time-series Expression Miner software (STEM) [36]. The clustered profiles with p-value ≤ 0.05 were considered as significantly expressed. Then the DEGs in all or in each profile were subjected to gene ontology (GO) classifications using WEGO [37], and KEGG pathway enrichment analysis.

Validation of differential expression using qRT-PCR
The cDNA was generated from 1 μg total RNA isolated from the fronds using a Prime-Script™ 1st Strand cDNA Synthesis Kit (TAKARA, Japan). Primers for quantitative real time PCR (qRT-PCR) were designed using Primer Premier 5.0 software (Premier, Canada) and synthesized by Sangon Biotech (Shanghai) Co., Ltd. The 18S (Gen-Bank accession number: KJ400889) was selected as reference. All the primers are shown in Additional file 1: Table S1. qRT-PCR was performed on a Bio-Rad iQ5 Optical System Real Time PCR System (Bio-Rad, USA). Each reaction mixture was 20 μL containing 10 μL of SYBR Green PCR Master Mix (TaKaRa, Japan), 250 nM of each primer, and 6 μL of diluted first-strand cDNAs. The qRT-PCRs were run as follows: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s in 96-well optical reaction plates. The Ct values were determined for three biological replicates, with three technical replicates for each value. Expression levels of the tested reference genes were determined by Ct values and calculated by 2 -△△Ct .

Statistical analysis
All data were statistically analyzed by means of the SPSS with LSD to identify differences. Significant differences (P < 0.05) between treatments are indicated by different letters.

Results
Phenotypic and physiological responses to NH 4 + toxicity in L. minor Figure 1 a-c shows changes in the appearance of L. minor fronds at the end of experiment. The fronds in NC looked green and healthy, as well as in A84. But in A840, some mother fronds looked greensick ( Fig. 1 c, shown by arrow). The RGR based on fronds number showed a downward trend from NC to A840 ( Fig. 1 e). This could indicate that the NH 4 + concentrations of 84 mg/L affected the propagation of L. minor, and the much higher concentration of 840 mg/L significantly inhibited the growth and could cause some damage.
Evans blue was used to determine the high NH 4 + -stress induced cell death (Fig. 1d). Almost no dead cell was stained in the fronds cultured in NC. Dead cells were however detected in both mother and newborn fronds in the plants grown in both tested NH 4 + concentrations, especially in A840.
Fronds of L. minor were stained with DAB or NBT to reveal in situ accumulation of two main reactive oxygen species (ROS), H 2 O 2 and O 2 − , respectively (Fig. 1d). Histochemically stained cells showed that the H 2 O 2 and O 2 − significantly accumulated in both the mother and newborn fronds in A840 after seven days. The fronds in A84 were also found to have some ROS accumulation. For the fronds in NC, the ROS was just slightly accumulated in some mother fronds that might be induced by the normal ageing.
The MDA content was used to detect the lipid peroxidation and membrane damage induced by oxidative stress. The contents of MDA in L. minor in A84 and A840 were higher than in NC, and the highest MDA content reached 4.4 nmol/g in A840 (Fig. 1f ). The activity of the antioxidant defense system was also analyzed ( Fig. 1g-j). Like the change of MDA, the activities of SOD and POD all increased from NC to A840, and the values all increased almost doubled. In contrast, the CAT decreased from NC (6.5 u/g) to A840 (3.8 u/g). For the APX, the highest activity was in A84 (45.28 u/g) and the lowest one was in A840 (8.96 u/g).

Overview of three libraries data sets by RNA-seq
As shown in Table 1, a total of 6.62G nucleotides, equivalent to 33,136,337 raw reads and 32,403,455 quality filtered (clean) reads were generated from the three separate libraries from NC, A84, and A840. The RNAseq generated clean reads ranged from 10.4 to 11.1 million on each sample. The Q20 percentages of the three libraries were from 97.21 to 97.44 %, and the GC contents ranged from 50.68 to 51.69 %. All clean reads were pooled together and then de novo assembled by Trinity. Based on chosen criteria, an average of 79.91 % of the clean reads was mapped, with perfect matches were from 47.03 to 47.09 %. In each library, the scales of clean reads uniquely mapped to the database were 76.87, 80.09 and 79.91 %, respectively. There were still approximately 20.09 % of clean reads that cannot be mapped back to any references, which could be due to the limited reference gene database of L. minor.
The final assembly of L. minor had 71,094 contigs with length ≥ 200 bp and after further removal of redundant sequences, 70,728 unigenes were obtained. The length of  unigenes ranged from 201 b to 14,857 b, with a mean size of 620 bp and N50 number of 988 bp (Table 2). In NC library, the number of genes identified increased with the number of reads until above 6 million, but 4 million in the other two libraries (Additional file 2: Figure S1A). The unigene coverage analyzed as a means of evaluating the quality of the RNA-seq data was mostly >50 %. More than half of the sequences have coverage more than 80 % (Additional file 2: Figure S1B). The amino acid sequences predicted by 'Getorf ' were searched by BLASTP. A final number of 29,171, 28,476, 17,209 and 14,172 unigenes (E-value < 1e −5 ) had significant matches in NR, KEGG and COG databases, respectively. As shown in Additional file 3: Figure S2, the sequences matched with the species in NR were determined as following: 27 % Vitis vinifera, 16 % Ricinus communis, 14 % Populus trichocarpa, 8 % Glycine max and 4 % Oryza sativa etc. The COG matched L. minor unigenes dataset were categorized into 25 functional COG clusters (Additional file 4: Figure S3). The five largest functional categories in which sequences were identified to include 1) Posttranslational modification, protein turnover, chaperones; 2) Signal transduction mechanisms; 3) General function prediction only; 4) Translation, ribosomal structure and biogenesis; and 5) Intracellular trafficking, secretion, and vesicular transport.

Identification and overview of the differentially expressed genes
We performed a pairwise comparison using NC as the control, and A84, or A840 as the treatments (Fig. 2a). Most of genes were down-regulated in the two treatments but the up-regulated genes in A840 were slightly higher than the down-regulated genes.
Results from FDR identification showed that 14,207 unigenes were classified as DEGs, which were then used for the subsequent analysis. All the 14,207 DEGs could be clustered into 7 profiles by STEM (Additional file 5: Figure S4; Additional file 6), in which 12,959 DEGs were further clustered into 3 profiles (p-value ≤ 0.05), including two down-regulated patterns (Profile 1 and Profile 0) and one up-regulated pattern (Profile 7) (Fig. 2b-d). Profile 1 and Profile 0 contained 11,625 and 954 DEGs, respectively, while Profile 7 contained 380 DEGs.
Next, the DEGs within the three profiles were subjected to GO-term analysis (Fig. 3). The DEGs were classified into three main categories including cellular component, biological process, and molecular function. Cell and cell parts under cellular component category were the two top abundant subcategories of the two down-regulated patterns (Profile 1 and Profile 0). For the up-regulated pattern of Profile 7, the metabolic process under molecular function was the top subcategories.

Analysis of phenylpropanoid biosynthesis pathway genes from L. minor unigenes
In plants, the phenylpropanoid biosynthesis pathway contributes to multiple biosynthetic branches, such as lignin and flavonoid biosynthesis. The expression of transcripts encoding for key enzymes for lignin and flavonoid biosynthesis were analyzed in this study (Fig. 4). The results showed that most of lignin biosynthesis related genes were up-regulated, but not for the expression of transcripts encoding for key enzymes for flavonoid synthesis.
As shown in Fig. 4, in lignin biosynthesis pathway, the genes of PAL (Phenylalanine ammonia-lyase), 4CL (4hydroxycinnamoyl-CoA ligase), COMT (caffeic acid Omethyltransferase), C3H (p-coumaroyl shikimate 3′-hydroxylase), F5H (coniferaldehyde/ferulate 5-hydroxylase),  This indicated that under NH 4 + stress, shift from caffeoyl CoA to feruloy CoA might be difficult during the biosynthesis of G and S type lignin, which are the main components of monocot lignin [38]. To get feruloy CoA, another way might be potentially utilized which involves the caffeic acid, which can be changed into ferulic acid by COMT, and subsequently changed into feruloy CoA. This mechanism has also been shown to be similar to other monocotyledon like switchgrass [39].

Expression profiles of oxidative-related and PCD-related genes
Expression of the 14 ROS oxidative-related genes including six oxidative marker genes, six ROS-scavenging genes and two ROS-producing genes are summarized in Fig. 5. The oxidative marker genes included a trypsin/  chymotrypsin inhibitor, a DNAJ heat shock protein, a FAD-binding protein and three cytochrome P450 genes, which are regarded as hallmarks for the general oxidative stress response [40,41]. The ROS-scavenging genes consisted of genes of CAT, SOD and POD. NADPH oxidases play a key role in generating ROS [42] and have been shown that RbohD and RbohF genes in Arabidopsis and RbohA of Hordeum vulgare could indeed induce ROS production [43,44]. In this study, two Rboh-like protein genes, the RbohA and RbohD in L. minor were up-regulated under high NH 4 + stress according to the RNA-seq and qRT-PCR results.
On the other hand, three genes indirectly related to oxidative stress were also detected. Under high NH 4 + stress, the nitrate reductase (NR) and non-symbiotic hemoglobin 1 (NSHB1) genes were all down-regulated, but the alternative oxidase 1 (AOX1) gene was upregulated (Fig. 5).
Metacaspases act as initiators and regulators for programmed cell death (PCD) in plants [45]. In DEGs, a metacaspase gene of MAC4 (metacaspase 4) was significantly up-regulated under NH 4 + stress. Conversely, two inhibitor of PCD, namely BAX inhibitor 1 and DAD1 (defender against cell death 1), were significantly down-regulated.  L. minor were currently deposited in public databases such as Genbank (accessed 12/10/2015). It is a challenge to do analysis and characterization of L. minor RNA-seq dataset without a sequenced genome, and in fact a lack of a sequenced genome in the Lemna. Until now, only the chloroplast genome of L. minor has been generated [46]. Despite such limitations, careful curation of the sequences and assembly using the robust Trinity Program allowed us to still identify 70,728 unigenes in L. minor, and more than 40 % of them were annotated. The differential expression analysis of RNA-seq data revealed that most of the DEGs were down-regulated under NH 4 + stress. GO and KEGG enrichment analysis revealed that the down-regulated genes (profile 1 and profile 0) mainly categorized into cellular structure and function, metabolic process and gene transcription indicating that some important physiological or "housekeeping" functions might have been inhibited. The same response has been observed in Arabidopsis under NH 4 + stress [21]. On the other hand, the up-regulated genes were mainly associated with the metabolic processes, especially the secondary metabolism, such as the Phenylpropanoid biosynthetic pathway (Table 2). It has been suggested that some secondary metabolites play an important role in defenses of abiotic stress [47]. For example, some flavonoids and lignin precursors have been reported to accumulate in response to various abiotic stresses [48].

NH 4 + toxicity caused oxidative stress and induced cell
death ROS is usually detected in overproduction under abiotic stresses, where it can cause some damages which ultimately results to oxidative stress [49]. Previous studies have shown that L. minor could suffer NH 4 + toxicity in high NH 4 + concentrations [16], and that NH 4 + concentration of 56 mg/L could induce oxidative stress [50]. In this study, we further explored the effects of two higher NH 4 + concentrations (84 and 840 mg/L) on oxidative damage of L. minor. Like other aquatic plants [51], the accumulated ROS and increased MDA content in L. minor in the two NH 4 + treatments indicated that oxidative damage occurring. This is because the increased ROS could induce oxidative stress that contributes to lipid peroxidation and membrane damage, and the MDA has been considered as the indicator of the damage [52]. In addition to the physical evidence, some molecular evidences on NH 4 + toxicity induced damage were also presented produced from the RNA-seq and qRT-PCR analysis (Fig. 5). The expression of some oxidative marker genes like trypsin/chymotrypsin inhibitor, DNAJ heat shock protein and cytochrome P450 genes were enhanced under high NH 4 + stress. Furthermore, the up-regulated ROS-producing genes of RbohA and RbohD also indicated oxidative stress occurring in L. minor under NH 4 + toxicity. The ROS scavenging enzymes play an important role in the plant's defense system in response to the generation of ROS. Among the enzymatic antioxidants, the enzyme SOD represents the first line of antioxidant defense by transforming O 2 − into H 2 O 2 , and then the APX, POD, and CAT subsequently metabolize H 2 O 2 [53]. Under NH 4 + toxicity, the SOD activity of L. minor increased but not the gene expression of the three types of SOD, indicating a lag from the gene transcription to enzyme action in this 7-day stress treatment. According to Huang's observations, the activity of SOD decreased until the 14th day under NH 4 + toxicity [50]. The activated SOD could transform O 2 − into H 2 O 2 but its removal only relies on POD for the gene expression and activities of CAT and APX which all decreased in A840.
Like ROS, nitric oxide (NO) also plays an important role in plant responses to environmental stress, and there are complex networks of interactions between ROS and NO when plants suffer oxidative stress [54]. In plants, NR could reduce nitrate to produce nitrite, as well as reduce nitrite to produce NO, which possesses antioxidant properties and likely to act as a signal in activating ROSscavenging enzyme activities under oxidative stresses [55]. The non-symbiotic hemoglobin could scavenge NO, thus building a futile cycle with NR [56]. The alternative oxidase can scavenge NO with ROS as the substrates, as well as prevent the production of excess ROS by stabilizing the redox state of the mitochondrial ubiquinone pool [56,57]. In this study, the NH 4 + is the sole N source in the two stress treatments, the down-regulated L. minor NR gene in A84 and A840 indicated the gene might not be activated without nitrate. The NR-mediated NO production might also be suppressed, even though the NSHB1 gene was also down-regulated. The slightly up-regulated AOX1 gene in L. minor may be involved in preventing ROS excessive increase under high NH 4 + stress.
ROS is one of the key regulators of PCD that is an active and genetically controlled form of cell death [58]. In this study, except for the ROS, the cell death was also detected in L. minor suffering from NH 4 + toxicity by staining. RNA-seq results further showed that a metacaspase gene, MAC4, was significantly up-regulated in the two NH 4 + treatments. In plants, the metacaspase is a discovered gene family that has distant caspase homologs closely related to PCD [59]. The MAC4 of Arabidopsis plays a positive regulatory role in abiotic stress-induced PCD [60]. In addition, our results also showed that the PCD inhibiters, like BAX inhibitor 1 [61] and DAD1 [62], significantly decreased in their gene expression. Thus, we can speculate that the NH 4 + toxicity induced PCD of L. minor, and that the ROS might play as an intermediate signaling molecule.
Lignin biosynthesis plays an important role in NH 4 + toxicity resistance Lignin is the major components of cell wall and the main structure in plant mechanical support and defense system [63]. There are two pathways for lignin biosynthesis in plants, namely of monolignol and phenylpropanoid pathways [64]. And the stimulation of the phenylpropanoid pathway has been considered as a common feature of some abiotic stress response such as drought, salinity, ozone intoxication and heavy metals [63,65]. Previous studies also showed that both nitrogen deficiency and fertilization (NH 4 NO 3 ) could induce a set of genes required for phenylpropanoid metabolism [66,67]. In this study, the RNA-seq and qRT-PCR results also showed enhanced expression of some key enzyme genes in phenylpropanoid pathway under high NH 4 + stress in L. minor. In addition, all up-regulated genes were lignin biosynthesis-related, rather than flavonoid synthesis, which could be due to the antagonistic relationships of the two biosynthetic pathways [68]. However, it could still be suggested that the NH 4 + toxicity could stimulate the phenylpropanoid pathway of L. minor, and lead to a shift of metabolism towards lignin.
G and S type lignin are the main components of monocot lignin [39]. In L. minor, a series of genes required for the biosynthesis of two types of lignin were up-regulated, including the rate-limiting enzymatic genes in lignin biosynthesis, like PAL [69] and F5H [70]. The increased lignin synthesis would result to higher lignin content, which together with other antioxidants, could play an important role in limiting ROS production in the apoplast [63]. This mechanism could be one of the reasons why L. minor could resist high NH 4 + stress.

Conclusions
In this study we report the first large transcriptome study carried out in L. minor where we have compared physiological and transcriptional responses to NH 4 + toxicity. Evidence from physiological observations, transcriptome and qRT-PCR analysis indicated that NH 4 + toxicity could induce ROS accumulation which in turn results to oxidative damage and later to cell death in L. minor. The antioxidant enzyme system was activated under NH 4 + toxicity for ROS scavenging. We also identified that the phenylpropanoid pathway was stimulated under NH 4 + toxicity, and the lignin biosynthesis was also up-regulated in this metabolic pathway. The increased lignin biosynthesis might play an important role in NH 4 + toxicity resistance.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data
The dataset is available from the NCBI Sequence Read Archive (SRA). The BioProject and SRA accession are PRJNA302233 and SRP066224, respectively.