Transcriptomic and physiological analyses of Medicago sativa L. roots in response to lead stress

Lead (Pb) is one of the nonessential and toxic metals that threaten the environment and human health. Medicago sativa L. is a legume with high salt tolerance and high biomass production. It is not only a globally important forage crop but is also an ideal plant for phytoremediation. However, the biological and molecular mechanisms that respond to heavy metals are still not well defined in M. sativa. In this study, de novo and strand-specific RNA-sequencing was performed to identify genes involved in the Pb stress response in M. sativa roots. A total of 415,350 unigenes were obtained from the assembled cDNA libraries, among which 5,416 were identified as significantly differentially expressed genes (DEGs) (false discovery rate < 0.005) between cDNA libraries from control and Pb-treated plants. Gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEGs showed they mainly clustered with terms associated with binding, transport, membranes, and the pathways related to signal and energy metabolism. Moreover, a number of candidate genes included antioxidant enzymes, metal transporters, and transcription factors involved in heavy metal response were upregulated under Pb stress. Quantitative real-time PCR(qRT-PCR) validation of the expression patterns of 10 randomly selected candidate DEGs were consistent with the transcriptome analysis results. Thus, this study offers new information towards the investigation of biological changes and molecular mechanisms related to Pb stress response in plants.


Introduction
Metal contaminants such as lead (Pb), cadmium (Cd), and mercury (Hg) are nonessential elements released into the environment by anthropogenic activities and are toxic to plants and can threaten human health through the food chain [1,2]. Pb is a highly toxic heavy metal that can be absorbed by plants; inhibit root and shoot growth; cause leaf chlorosis; and disturb other physiological processes such as mitosis, transpiration, and DNA synthesis [3,4]. Physiological disorders and ultimately death can result from the production of large quantities of PLOS ONE | https://doi.org/10.1371/journal.pone.0175307 April 7, 2017 1 / 16 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The de novo assembly of all the clean reads generated 508,781 contigs with an average length of 473.95 bp and N50 length of 519 bp ( Table 2). The minimum and maximum lengths were 201 bp and 14,904 bp, respectively. In total, 415,350 unigenes were obtained with an average length of 426.68 bp and N50 length of 436 bp. RNA-seq data quality analysis showed more than 80% of unigenes had coverage above 50% for each library (S1 Fig).

Identification of differentially expressed genes (DEGs)
To identify unigenes induced by Pb stress, the five libraries were first divided into a control and a treatment group, C0-C48-C96 and C0-Pb48-Pb96, and a comparison of the two groups identified 54,173 unigenes that were specifically expressed in the C0-Pb48-Pb96 group. Next, a pairwise comparison was performed between libraries C48 and Pb48, and C96 and Pb96, to identify DEGs at specific times after Pb exposure. A total of 31,725 DEGs were identified between Pbtreated and control libraries (Fig 2A). Finally, a comparison between the 54,173 unigenes specifically expressed in the Pb-treated library and the 31,725 DEGs identified 5,416 upregulated DEGs under Pb stress (Fig 2B, S2 Table). Among these 5,416 DEGs, 3,609 (66.64%) were annotated; of these DEGs, 1,749 unigenes were mapped to the gene ontology (GO) database and 792 of these unigenes mapped to 196 KEGG pathways (S2 Table). Of the 5,416 DEGs specific to Pb treatment, 3,215 were specifically expressed in Pb48, 1,265 were specifically expressed in Pb96, and the last 936 were expressed in both Pb48 and Pb96 libraries (S2 Table). The large number of DEGs indicates that Pb treatment alters gene expression across the Medicago genome, which may result in specific stress responses to Pb exposure. In addition, the differences of DEGs specifically expressed in Pb48 or Pb96 suggest extensive biological changes in Medicago after Pb treatment.

Candidate genes related to Pb stress
Pb stress leads to an increased level of ROS in plants, and thus the activities of antioxidative enzymes such as CAT, POD, and SOD increase in plants exposed to Pb stress [28,29]. Here, eleven DEGs coding for POD and four DEGs coding for SOD were predicted to be involved in Pb stress (S5 Table). In addition to the antioxidant enzymes, the biosynthesis of antioxidants was also elevated in response to Pb stress in M. sativa. Glutathione is one of the most important plant antioxidants and is also a substrate for phytochelatins, which are metal-chelating compounds that enable metal tolerance [4,27]. We identified 15 DEGs coding for glutathione Stransferase (GST) in the "glutathione metabolism" (ko00480) pathway (S5 Table). Production of secondary metabolites like flavonoids and isoflavonoids is an effective mechanism adopted by plants against heavy metal toxicity [4]. In this study, 14 DEGs involved in flavonoid and isoflavonoid synthesis, including three chalcone synthase (CHS) genes that are involved in the first step of flavonoid biosynthesis, were identified in the flavonoid-related pathways (ko00940, ko00941, ko00942, ko00943, and ko00944) (S5 Table).
Metal transporters and stress response proteins play vital roles in heavy metal tolerance in plants [19]. In this study, we identified a number of DEGs coding for metal transporters, including 16 DEGs for ABC transporters, three DEGs for IRT, three DEGs for CDFs, and two DEGs for heavy metal transport (S5 Table). Signal transduction and expression of TFs are also responsive to stress in plants. We identified 13 DEGs for calcium-binding proteins or MAPKs. Additionally, 24 DEGs for TFs in the WRKY family, MYB ethylene-responsive transcription factor (ERF), and bZIP were identified in Pb-treated M. sativa (S5 Table).
Previous studies have shown that photosynthesis is sensitive to heavy metal stress [27]. In this study, the KEGG pathways for "photosynthesis-antenna proteins" (ko00196) were significantly enriched in Pb-stress cDNA libraries (Table 3), in which the expression levels of 10 DEGs coding for chlorophyll a/b-binding protein (light-harvesting complex I protein) were significantly upregulated (S5 Table). The DEGs that were upregulated in the "photosynthesis-antenna proteins" pathway are shown in Fig 4.

qRT-PCR validation of candidate DEGs
To evaluate the quality of the transcriptome data generated by the Illumina high-throughput sequencing platform, ten candidate DEGs that were predicted to be involved in the Pb stress response were selected for validation using qRT-PCR analysis (S1 Table). For these genes in control and Pb-treated plants, the expression levels calculated by qRT-PCR were consistent with the expression levels from the transcriptome sequencing data (Fig 5). This indicates that our transcriptome sequencing data are reliable.

Discussion
Pb is now a widespread soil contaminant all over the world and poses a great threat to food safety and human health [30,31]. Different physical, chemical, and biological methods have been employed to address heavy metal pollution in the environment, among which phytoremediation is thought to be an efficient, cost-effective, and environmentally friendly method [32,33]. However, the genes involved in protection against heavy metal toxicity and the molecular mechanisms underlying this protection are still undefined in plants. In this study, we identified a number of genes and pathways related to Pb stress responses in M. sativa using transcriptome sequencing, which improves our current knowledge of heavy metal responses and tolerance in plants.
De novo transcriptome sequencing and assembly had been widely used to evaluate gene expression levels, discover novel genes and alternative splicing in non-model plants. Here [37]. N50 value is an indicator of the transcript contiguity of de novo assembly. It is affected by several factors, such as species, assembly software, and assembly parameters, especially k-mer value [36,37]. Thus, the different N50 values among studies may result from different software and parameters. Although the mean length and N50 value of our study was relatively lower than many other studies, this did not affect the accuracy of the transcriptome, as 71.86% of the unigenes match to the genome of M. truncatula, the closely related species of M. sativa and the model leguminous plant.
The transcript levels of genes are used to reflect the biological differences between tissues and plant developmental stages. DEG identification and the subsequent GO and KEGG enrichment analyses of the cDNA libraries from control and Pb-treated plants can be directly related to the global biological changes and molecular mechanisms that respond to Pb in M. sativa roots. In this study, we identified 5,416 DEGs that were specifically expressed in the Pbtreated roots. The GO term and KEGG enrichment analyses showed that terms for transport (GO:0006810) and membrane (GO:0016020, GO:0031224) (S3 Table), and pathways related to signal (ko03320, ko04075) and energy metabolism (ko00040, ko00020, ko00196, ko00010) ( Table 3), were significantly enriched in the identified DEGs. These results were in agreement with previous studies that reported that heavy metals damaged plant membrane systems, inhibited photosynthesis, and delayed plant growth [27,[38][39][40].
The plasma membrane is the primary environmental barrier for a plant cell and mediates the exchange of information and materials between the cell interior and the extracellular environment [41,42]. Thus, the enzymes and transporters located on the plasma membrane play essential roles in heavy metal tolerance and detoxification. The membrane proteins in the natural resistance-associated macrophage proteins and ABC transporter families, and ATPases, are involved in metal tolerance by transporting metals to the vacuole or pumping them back out of the cells [39,[43][44][45]. Additionally, membrane-associated signaling proteins, such as MAPKs and calcium-binding related protein kinases, are involved in several stress responses and are activated under heavy metal stress [20,46]. In this study, a number of DEGs coding membrane proteins associated with metal transport and signaling were identified, including 16 DEGs for metal transporters (ABC transporters, CDFs, and zinc-regulated transporter, IRT-like proteins), seven DEGs coding for chlorophyll a/b-binding proteins (light-harvesting complex), and 13 DEGs coding for signaling proteins (MAPKs and calcium-binding protein) (S5 Table). This indicates that the membrane and its associated transport proteins are important in Pb response and tolerance in M. sativa.
Photosynthesis is the most important biological process for plants, and heavy metals like Pb inhibit photosynthesis by decreasing the chlorophyll content and number of chloroplasts, and disrupting chloroplast structure [27,47,48]. The chloroplast chlorophyll a/b-binding protein is an essential factor in photosynthesis that participates in the two key processes of light absorption and electron transfer [27,49]. In this study, the DEGs related to "photosynthesisantenna proteins" (ko00196) in our pathway prediction analysis were significantly enriched. Among the 10 DEGs coding for the chlorophyll a/b-binding protein, nine were upregulated after 48 h treatment with Pb (S5 Table). However, they were downregulated after 96 h, which may indicate that Pb stress stimulates the expression of chlorophyll a/b-binding protein, while the excess accumulation of Pb in plant cells can lead to an inhibition of chlorophyll a/b-binding protein expression and affect photosynthesis.
Antioxidant synthesis is a main mechanism for protection against heavy metal toxicity in plants through their ROS-scavenging activity [2,50]. Antioxidant enzymes such as SOD, POD, and CAT play essential roles in controlling the ROS level in cells [51]. Studies of heavy metal tolerance in different plant species show that the activities of SOD, POD, and CAT significantly increase after treatment with toxic metals [28,29,52,53]. In this study, we demonstrated that the activities of SOD and POD significantly increased in the roots of Pb-treated M. sativa (Fig 1). In addition, the transcriptome data showed that 11 DEGs for POD and four DEGs for SOD were upregulated in Pb-treated roots (Fig 3, S5 Table). These results show that antioxidant enzymes act to reduce ROS damage in M. sativa after heavy metal exposure.
In addition to antioxidant enzymes, many kinds of metabolic products are efficient antioxidants in plant cells. Flavonoids and glutathione are well-documented antioxidants that protect against heavy metal toxicity [4,19,50,54]. Flavonoids act as ROS scavengers and metal chelators, and isoflavonoids act as phytoalexins in response to heavy metals [4]. Glutathione plays an important role in heavy metal detoxification [19,50] and is not only an important antioxidant but also participates in metal transport processes [4,45,55]. In this study, both transcriptome data and qRT-PCR analysis indicated that DEGs that encode flavonoid synthase genes were significantly upregulated under Pb stress (Fig 5, Table 3, S5 Table). Moreover, in accordance with previous reports, the DEGs for enzymes involved in the "glutathione metabolism" (ko00480) pathway were significantly upregulated in plants under Pb stress. Taken together, our results indicate that the antioxidant reactions are significantly activated to protect against Pb toxicity in M. sativa roots.
In summary, the de novo transcriptome sequencing of M. sativa roots obtained 415,350 unigenes, among which 5,416 were defined as DEGs that were specifically expressed in Pb-treated cDNA libraries. GO and KEGG pathway enrichment analyses of these DEGs showed that these genes mainly clustered with terms for binding, transport, and membrane, and the pathways were related to signaling and energy metabolism. Furthermore, a number of genes involved in heavy metal response were successfully identified. This study adds to the foundational knowledge of the biological changes and molecular mechanisms that respond to Pb stress in M. sativa roots. In addition to the previous studies on heavy metal response in different plant species, this will promote the use of phytoremediation of heavy metal-contaminated soils.

Plant material, growth conditions, and treatments
Alfalfa seeds (M. sativa L. 'Longdong') were surface-sterilized, rinsed, and germinated in the dark at 26˚C. Seedling were then grown in soil and supplied with 1/2 strength modified Hoagland nutrient solution [52]. Seedlings were grown under a 14 h light/10 h dark photoperiod with 70% humidity at 26˚C. After 30 d, the seedlings were divided into two groups and treated with 0 or 200 mg/L Pb 2+ as Pb(NO 3 ) 2 in 1/2 strength modified Hoagland nutrient solution for 0, 24, 48, 72, 96, or 120 h. The roots of seedlings were then harvested, surface cleaned, and frozen immediately in liquid nitrogen before being stored at -80˚C for future analysis.

Assay of soluble protein content and enzyme activity
The soluble protein content and antioxidant enzyme activities of SOD, POD, and CAT in roots of the Pb-treated and control groups were analyzed at each time point. Frozen root tissues were homogenized in 0.05 M phosphate buffer (pH = 7.8) and the homogenate was centrifuged at 15,000 ×g for 20 min at 4˚C. The supernatant was used as a crude sample for enzymatic activity assays.
The soluble protein content was determined as descripted by Bradford [56], using Coomassie brilliant blue G-250 as a dye and albumin as a standard. The activity of SOD and CAT were assayed according to Li et al. [57]. One unit of SOD was defined as the enzyme amount that caused 50% inhibition reduction in nitro-blue tetrazolium. The enzyme activity of CAT was expressed as micrograms of H 2 O 2 destroyed per minute per gram fresh weight (FW) root material. The activity of POD was assayed according to Fang and Kao [58]. A unit of POD specific activity was defined as enzyme units per gram FW root material.

Library construction and sequencing
Total RNA from roots of Pb-treated and control plants at 0, 48, and 96 h (notated as Pb48, Pb96, C0, C48, and C96, respectively) were isolated using an RNAprep Pure Plant kit (Bio-Teke, China). The concentration and quality of RNA were assessed using Qubit 3.0 (Thermo Scientific, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, USA), respectively. For stranded RNA-seq, cDNA libraries were prepared using an NEBNext 1 Ultra™ Directional RNA Library Prep Kit (NEB, UK) according to the manufacturer's instructions. Briefly, the mRNA was isolated and fragmented into 250-450 bp before synthesis of the first strand cDNA. The second strand of cDNA was then synthesized by adding dUTP as marker. Finally, the double strand cDNA was digested with Uracil-DNA Glycocasylase (UDG) before PCR reaction and sequencing. Thus, only the first strands of cDNA were kept and sequenced. The libraries were sequenced on a Hiseq 4000 (Illumina) using a paired-end run (2 × 150 bp).
The raw transcriptome data has been deposited at the NCBI database with Short Read Archive (SRA) accession numbers of SRR5279707-SRR5279711.

Data processing and assembly
The raw reads were first quality-filtered by removing adapter sequences and reads with a quality under Q20 by using the FastQC tool (http://www.bioinformatics.bbsrc.ac.uk/projects/ fastqc/). The total clean reads from the five libraries were then processed and de novo assembled with Trinity (Version 2.2.0) [59] using the default parameters. Trinity assembled the PE reads as contigs with a fixed k-mer value of 25. Contigs overlapped and reads that astride contigs were assigned into a same group, and then transcriptome was assembled by using de Bruijn graph strategy. The longest transcript extracted from each de Bruijn graph was defined as a Unigene, and these transcripts were used as reference sequences for subsequent analyses. The coding and protein sequences of unigenes were predicted by TransDecoder, a plug-in for Trinity (http://transdecoder.github.io/).

Functional annotation of unigenes
The CDS and amino acid sequences of unigenes were aligned to the Gene Ontology (GO) database, the Cluster of Orthologous Groups of proteins (COG) database, the NCBI nonredundant (NR) protein database, the Swiss-Prot protein database, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, respectively, by using BlastX program to identify functional genes (E-value < le -5 ).

Identification and functional enrichment of DEGs
Gene expression levels were calculated and normalized using the FPKM method. DEGs between different libraries were identified using a Student's t-test, and the P-value was adjusted using multiple testing procedures according to Benjamini and Yekutieli [60], by controlling the FDR. Finally, the significant DEGs were defined with the threshold of FDR < 0.005 and absolute value of log 2 (fold-change) ! 1.
The GO and KEGG enrichments analyses for functional significance was performed using an ultra-geometric test with Benjamini-Hochberg correction [60]. GO terms for significant enrichment of DEGs were defined as corrected P-value < 0.05 when compared to the genome background. Pathways for significant enrichment of DEGs were defined as Q value < 0.05.

Quantitative RT-PCR analysis
For quantitative RT-PCR (qRT-PCR) analysis, 1 μg of root total RNA from the control plants at 0 h 48 h, and 96 h and the Pb-treated M. sativa at 48 h and 96 h were treated with RNasefree DNase I before being used as a template for reverse transcription with a RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA). Quantitative RT-PCR analysis was performed using a SYBR 1 Green master mixture (BioRad, USA) in a LightCycler 1 96 (Roche, USA). Genes and primers are listed in supplementary S1 Table. Supporting information S1 Fig