Abstract Background Plants can suffer ammonium (NH[4]^+) toxicity, particularly when NH[4]^+ is supplied as the sole nitrogen source. However, our knowledge about the underlying mechanisms of NH[4]^+ toxicity is still largely unknown. Lemna minor, a model duckweed species, can grow well in high NH[4]^+ environment but to some extent can also suffer toxic effects. The transcriptomic and physiological analysis of L. minor responding to high NH[4]^+ may provide us some interesting and useful information not only in toxic processes, but also in tolerance mechanisms. Results The L. minor cultured in the Hoagland solution were used as the control (NC), and in two NH[4]^+ concentrations (NH[4]^+ was the sole nitrogen source), 84 mg/L (A84) and 840 mg/L (A840) were used as stress treatments. The NH[4]^+ toxicity could inhibit the growth of L. minor. Reactive oxygen species (ROS) and cell death were studied using stained fronds under toxic levels of NH[4]^+. 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 NH[4]^+ 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 NH[4]^+ toxicity. Conclusions The first large transcriptome study in L. minor responses to NH[4]^+ toxicity was reported in this work. NH[4]^+ toxicity could induce ROS accumulation that causes oxidative damage and thus induce cell death in L. minor. The antioxidant enzyme system was activated under NH[4]^+ toxicity for ROS scavenging. The phenylpropanoid pathway was stimulated under NH[4]^+ toxicity. The increased lignin biosynthesis might play an important role in NH[4]^+ toxicity resistance. Electronic supplementary material The online version of this article (doi:10.1186/s12870-016-0774-8) contains supplementary material, which is available to authorized users. Keywords: Lemna minor, NH[4]^+ toxicity, RNA-seq, Transcriptome, Oxidative damage, Lignin biosynthesis, Phenylpropanoid pathway, Oxidative damage, Programmed cell death Background Ammonium (NH[4]^+) and nitrate (NO[3]^−) are the two inorganic nitrogen (N) forms that can be directly absorbed by plants [[33]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]^+ [[34]2, [35]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 [[36]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 [[37]4–[38]6]. NH[4]^+ toxicity is not only a significant ecological issue, but also an important plant physiological process [[39]7]. Plant scientists have been trying to reveal its occurrence, signal transmission and physiological targets in plants [[40]3, [41]7, [42]8]. Usually, for most plants, the root bears the brunt of NH[4]^+ toxicity [[43]7, [44]9–[45]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 [[46]9, [47]12–[48]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 [[49]7, [50]8]. Additionally, recent studies also revealed that the NH[4]^+ toxicity could break the intracellular pH balance and C/N balance [[51]3, [52]7, [53]15], and cause oxidative damage [[54]16, [55]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 [[56]8, [57]18]. Transcriptome analysis is an effective method for global expression profiling of genes involved in stresses and toxicity in living organisms [[58]19, [59]20]. For example, transcriptomic profiling using microarrays have been used in Arabidopsis to identify molecular changes involved in NH[4]^+ toxicity [[60]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 [[61]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 [[62]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 [[63]24, [64]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 [[65]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 [[66]7, [67]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. Methods Sample preparation Samples were prepared as described in Wang et al. [[68]15]. L. minor was collected from a eutrophic pond in Chengdu, Sichuan, China (location: 30° 38.86′N, 104° 18.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. [[69]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 [[70]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 [[71]28, [72]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 [[73]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 [[74]31]. Cell death was examined by Evans blue staining as described by Kim et al. [[75]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 second-strand 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 ([76]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 [[77]33]. Unigenes were defined after removing redundancy and short contigs from the assembly. The unigenes were predicted by “GetORF” in the EMBOSS package [[78]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 [[79]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) [[80]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 [[81]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 (GenBank accession number: [82]KJ400889) was selected as reference. All the primers are shown in Additional file [83]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 [84]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. [85]1 c, shown by arrow). The RGR based on fronds number showed a downward trend from NC to A840 (Fig. [86]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. Fig. 1. Fig. 1 [87]Open in a new tab Phenotypic and physiological responses of Lemna minor in NC, A84 and A840. a-c, the appearance of L. minor in NC, A84 and A840, respectively, red arrows showed the greensick fronds, scale bar 5 mm; d Histochemically staining of cell death, O[2] ^− and H[2]O[2] by Evans blue, nitroblue tetrazolium (NBT) and 3,3′-diaminobenzidine (DAB), respectively; e relative growth rate (RGR) based on fronds number; f MDA contents; g-j, enzyme activity of catalase (CAT), superoxide dismutase (SOD), peroxidase (POD) and ascorbate peroxidase (APX), respectively Evans blue was used to determine the high NH[4]^+-stress induced cell death (Fig. [88]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. [89]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. [90]1f). The activity of the antioxidant defense system was also analyzed (Fig. [91]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 [92]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 RNA-seq 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,