Functional and Integrative Analysis of the Proteomic Profile of Radish Root under Pb Exposure

Lead (Pb) is one of the most abundant heavy metal (HM) pollutants, which can penetrate the plant through the root and then enter the food chain causing potential health risks for human beings. Radish is an important root vegetable crop worldwide. To investigate the mechanism underlying plant response to Pb stress in radish, the protein profile changes of radish roots respectively upon Pb(NO3)2 at 500 mg L−1(Pb500) and 1000 mg L−1(Pb1000), were comprehensively analyzed using iTRAQ (Isobaric Tag for Relative and Absolute Quantification). A total of 3898 protein species were successfully detected and 2141 were quantified. Among them, a subset of 721 protein species were differentially accumulated upon at least one Pb treatment, and 135 ones showed significantly abundance changes under both two Pb-stressed conditions. Many critical protein species related to protein translation, processing, and degradation, reactive oxygen species (ROS) scavenging, photosynthesis, and respiration and carbon metabolism were successfully identified. Gene Ontology (GO) and pathway enrichment analysis of the 135 differential abundance protein species (DAPS) revealed that the overrepresented GO terms included “cell wall,” “apoplast,” “response to metal ion,” “vacuole,” and “peroxidase activity,” and the critical enriched pathways were involved in “citric acid (TCA) cycle and respiratory electron transport,” “pyruvate metabolism,” “phenylalanine metabolism,” “phenylpropanoid biosynthesis,” and “carbon metabolism.” Furthermore, the integrative analysis of transcriptomic, miRNA, degradome, metabolomics and proteomic data provided a strengthened understanding of radish response to Pb stress at multiple levels. Under Pb stress, many key enzymes (i.e., ATP citrate lyase, Isocitrate dehydrogenase, fumarate hydratase and malate dehydrogenase) involved in the glycolysis and TCA cycle were severely affected, which ultimately cause alteration of some metabolites including glucose, citrate and malate. Meanwhile, a series of other defense responses including ascorbate (ASA)–glutathione (GSH) cycle for ROS scavenging and Pb-defense protein species (glutaredoxin, aldose 1-epimerase malate dehydrogenase and thioredoxin), were triggered to cope with Pb-induced injuries. These results would be helpful for further dissecting molecular mechanism underlying plant response to HM stresses, and facilitate effective management of HM contamination in vegetable crops by genetic manipulation.


INTRODUCTION
Heavy metal (HM) contamination from natural, agricultural, and industrial sources has become a worldwide public health concern, which could seriously deteriorate the environment and cause adverse impacts on human health through the food chain (Pourrut et al., 2013;Singh et al., 2016). Lead (Pb) is one of the most abundant HM pollutants with no physiological function, mostly penetrating the plant through the roots and accumulating in different parts (Gupta et al., 2013;Pourrut et al., 2013). Pb can influence various morphological, physiological and biochemical processes leading to a restriction of plant growth, inducing membrane cell damages and, ultimately, to cell death (Sharma and Dubey, 2005;Gupta et al., 2011). Accordingly, plants also have developed diverse defense mechanisms to prevent the toxic effect of heavy metals including generating signal sensing and transduction proteins, compartmenting into the vacuoles and induction of higher levels of metal chelates like a protein complex, organic and inorganic complexes (Gupta et al., 2011(Gupta et al., , 2013Thapa et al., 2012).
Radish (Raphanus sativus L.), a major member of the Brassicaceae family, is an important annual or biennial root vegetable crop worldwide (Wang and He, 2005). Because the root is considered as the vulnerable part which is easily affected by HM (Wang et al., 2015a), it has become of vital importance to investigate the HM-response mechanisms and explore the molecular regulatory network of tolerance and homeostasis in radish. The identification of the HM-responsive genes or proteins is a fundamental step in understanding the molecular mechanism underlying response to HM stress (Ahsan et al., 2009). In our previous study, NGS-based transcriptome, miRNA and degradome analysis were employed to investigate the expression patterns of genes and miRNAs in radish exposure to Pb stress. A lot of Pb-responsive transcripts, miRNA and its targets were detected, which were predominately involved in stressrelated signal sensing and transduction, specific metal uptake and homeostasis, glutathione metabolism-related processes and carbohydrate metabolism-related pathways (Wang et al., , 2015b. Although transcriptomics provides a useful tool for unraveling gene expression networks at the mRNA level which enhanced our understanding the response of radish under Pb stress, proteomics can offer a new platform for investigating complex biological functions involving large numbers of proteins and provide further insight into posttranscriptional modifications thereby complementing genomics analysis (Ralhan et al., 2008;Ahsan et al., 2009;Wang et al., 2014). In the past two decades, classical two-dimensional electrophoresis (2-DE) technology has been widely employed for protein identification and analysis. However, there were some limitations in its applications such as poor reproducibility, weak sensitivity and low automation (Sazuka et al., 2004). Isobaric tags for relative and absolute quantification (iTRAQ) is a robust mass spectrometry of protein quantitative technology, which can perform relative and absolute quantification in up to eight samples in parallel (Bindschedler and Cramer, 2011;Glibert et al., 2014). Recently, iTRAQ has been widely used for large-scale quantitative plant proteomic studies in exploration of various metabolic processes at the post-transcriptional level (Kambiranda et al., 2013;Martínez-Esteso et al., 2013) in response to various stresses Li et al., 2015;Fu et al., 2016). Additionally, iTRAQbased proteomics has been proven as a powerful method for unraveling the molecular regulatory networks involved in interactions between heavy metals and plant species, and a set of HM-responsive candidate proteins have been successfully identified. Alvarez et al. (2009) reported that exposure of Brassica juncea roots to cadmium (Cd) could activated several protein species involved in sulfur assimilation, redox homeostasis and xenobiotic detoxification, and depressed multiple proteins involved in protein synthesis and processing by two quantitative proteomics approaches including fluorescence two-dimensional differential gel electrophoresis (DIGE) and iTRAQ technology. More recently, the protein abundance changes in rice roots in response to Aluminum (Al) at an early phase were conducted with iTRAQ, and a total of 700 distinct protein species with >95% confidence were identified and 106 protein species were differentially accumulated upon Al toxicity in sensitive and tolerant cultivars (Wang et al., 2014). However, investigation of the dynamically protein abundance changes in response to Pb stress in radish has not been reported.
In this study, an iTRAQ-based quantitative proteomics approach was firstly employed to detect the effects of Pb responses at the protein abundance levels in radish roots. The differential abundance protein species (DAPS) involved in Pbresponse of radish were quantified, and the enriched networks for regulating Pb stress at the protein level were acquired. Furthermore, to deeply reveal the integrative molecular network of radish plant response to Pb stress, the proteomic data were integrated with our previous transcriptomic, miRNA, degradome and metabolomic data, which provided a more global view of the molecular and cellular changes elicited by Pb stress in radish. This work would provide valuable information for further functional analyses of the critical Pb-responsive protein species in radish, which will be helpful for effectively facilitating the management of Pb and other HM contaminations in vegetable crops by genetic manipulation.

Plant Material
The variety of high-Pb-accumulation "NAU-RG" was selected for exploring the molecular regulation mechanisms in radish roots responding to the Pb stress. This genotype is an advanced inbred line with a medium size root in globular shape, white skin and flesh. According to the related evaluation methods and previous studies, different concentrations of Pb (NO 3 ) 2 (100, 200, 400, 500, 1000, and 1500 mg·L −1 ) were set to investigate the changes of visible physiological symptoms with different temporal durations in the preliminary experiment. Interestingly, no special obvious morphologic differences were found among individuals when exposed to low dose of heavy metal treatment, while the plants were seriously hampered and grew abnormally when exposed to 1500 mg·L −1 Pb (NO 3 ) 2 . Therefore, the concentrations of Pb (NO 3 ) 2 at 500 mg·L −1 (Pb500) and 1000 mg·L −1 (Pb1000) were selected for the comparative proteomic analysis. Additionally, a control group was defined using non-treated seedlings (CK).The growth conditions of radish plants were conducted according to the reported descriptions . Plants were collected after 72 h with three different concentration treatments including an untreated control (CK) and two Pb-stressed conditions (Pb500 and Pb1000), based on the previous reported study . Each treatment consisted of three biological replicates. Equal amount of radish taproot samples (2.5 g) from three randomly selected individual plants of each replicate were pooled and cut into small pieces, which were rapidly frozen in liquid nitrogen and stored at −80 • C for protein extraction.

Protein Extraction
Total protein was extracted from the radish taproots using a phenol (Phe) extraction procedure according to the reported methods with some modifications (Saravanan and Rose, 2004;Yang et al., 2013). In brief, about 2.5 g frozen taproots were finely powdered in liquid nitrogen adding 0.5 polyvinylpolypyrrolidone (PVPP) and then suspended in 20 mL of precooled extraction buffer containing 500 mM Tris-HCl pH 7.5, 50 mM ethylenediaminetetraacetic acid (EDTA), 100 mM KCl, 2 mM dithiothreitol (DTT), 2 mM phenylmethylsulfonyl (PMSF) and 30% sucrose (w/v). The mixture was extensively homogenized on ice for 2 min, and then an equal volume of Tris-HCl pH 7.5-saturated Phe was added. The mixture was thoroughly vortexed and proteins were collected by centrifuging at 12,000 g for 15 min at 4 • C. The upper Phe phase was removed and re-extracted two or three times with extraction buffer. Proteins were precipitated from the final Phe phase with five volumes of saturated ammonium acetate in methanol at −20 • C overnight. After centrifuging at 12,000 g for 15 min at 4 • C, the pellets were washed twice with 20 mL of 0.1 mol·L −1 cooled acetone, and then dried by lyophilization and finally stored at −80 • C until use. The protein was quantified by a Bradford protein assay kit (Sangon Biotech, China), and an equal amount of protein from three replicates of each treatment (CK, Pb500, and Pb1000) was respectively pooled for iTRAQ analysis.

iTRAQ Labeling and Strong Cation Exchange (SCX) Chromatography
The protein samples were dissolved in 100 mM triethylammonium bicarbonate (TEAB, pH 8.5) containing 1% SDS (w/v), reduced with 10 mM DTT at 56 • C for 1 h, and followed by alkylation with 55 mM iodoacetamide (IAM) for 45 min at room temperature in the dark. Trypsin was then added to a final enzyme/substrate ratio of 1:20 (w/w) for protein digestion, which was incubated at 37 • C for 12 h. The resulting tryptic peptides were vacuum-concentrated and then labeled with iTRAQ reagents (Applied Biosysterms, USA) according to manufacturer's instruction. The control sample (CK) was labeled with 118 iTRAQ reagent, and Pb-stressed samples treated with 500 and 1000 mg·L −1 Pb (NO 3 ) 2 were labeled with 119 and 114, respectively. The labeling reactions were incubated for 2 h at room temperature. Subsequently, all the peptides from three groups were combined and further fractionated using SCX chromatography on a Ultremex SCX column (250 × 4.6 mm, 5 µm particle size, 200 Å pore size, Phenomenex, USA) by high performance liquid chromatography (HPLC) system (Shimadzu LC-20AB, Japan). The HPLC gradient consisted of buffer A (25 mM NaH 2 PO4, 25% ACN, pH 2.7) for 10 min, 5−35% buffer B (25 mM NaH 2 PO 4 , 25% ACN, 1M KCL, pH 2.7) for 11 min and 35−80% buffer B for 1 min, which were eluted at a flow rate of 1 mL·min −1 . The chromatograms were recorded at 214 nm and 20 constituents were obtained. The collected fractions were desalted with StrataX (Phenomenex, USA) and concentrated to dryness using a vacuum centrifuge.

Liquid Chromatography Coupled with Tandem Mass Spectrometry (LC-MS/MS)
The mass spectroscopy analysis was performed using an ABSCIE X TripleTOF 5600 mass spectrometer (AB SCIE X, Framingham, MA, USA), coupled with an online flow HPLC nanoAcquity system (Waters, USA), which comprised two parts of columns including Symmetry C18 column (5 µm particles, 180 um × 20 mm) and BEH130 C18 column specification (1.7 µm particles, 100 µm × 100 mm). Symmetry C18 column was used for adsorption and desalination of peptides, while BEH130 for the separation of peptides. The correction liquid (Thermo Fisher Scientific) was respectively added to liquid A (water:acetonitrile:formic acid = 98:2:0.1) and B (water:acetonitrile:formic acid = 2:98:0.1) with a certain proportion. A 2.25 µg aliquot of supernatant (9 µl) was transferred, and then the samples were eluted with liquid A at 2 µl·min −1 for 15 min for the purpose of adsorption and desalination of peptides. After that, the mobile phase which contain 5% liquid B was used to elute the supernatant (300 nl·min −1 , 1 min). Next, the concentration of liquid B was raised from 5% to 35% in 40 min, from 35% to 80% in 5 min using linear gradients, and then the samples were eluted for 5 min with 80% liquid B. The separation was performed at a constant flow rate of 20 µl·min −1 .

Data Analysis and Protein Identification
Raw data files were converted into MGF files using Proteome Discoverer 1.2. Proteins identification were performed by using Mascot search engine (Matrix Science, London, UK; version 2.3.02) against our sequenced radish root transcriptome database (SRX316199 and http://www.ncbi.nlm.nih.gov/sra/; Wang et al., 2013).The search parameters were set according to the published studies (Qiao et al., 2012;Yang et al., 2013;Wang et al., 2014). For protein identification, a mass tolerance of 0.05 Da was permitted for intact peptide masses and 0.1 Da for fragmented ions, with allowance for one missed cleavages in the trypsin digests. The methionine oxidation (M), protein N-terminal acetylation, deamidation (N, Q), and iTRAQ (Y) were selected as variable modifications, and carbamidomethyl modification of cysteines (C), iTRAQ (Nterminal), iTRAQ (K) were as fixed modifications. The search results were passed through additional filters before exporting the data. For protein identification, the filters were set as follows: significance threshold P, 0.05 (with 95% confidence) and ion score or expected cutoff <0.05 (with 95% confidence). For protein quantitation, it was required that a protein contains at least two unique spectra. The quantitative protein ratios were weighted and normalized by the median ratio in Mascot. The DAPS were defined: Showing a fold-change of greater than 1.2 or less than 0.83 using a greater statistically significant value P < 0.05.

Functional Analysis
The online software BLAST2GO (http://www.blast2go.com/ b2ghome) was used to automatically assign protein description and obtain annotations from homologous sequences of public databases (Conesa et al., 2005;Conesa and Götz, 2008). A metabolic pathway analysis was undertaken based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The Gene Ontology (GO) and metabolic pathway enrichment analysis of the DAPS were conducted through two publically available tools (webservers), namely, DAVID 6.7 and KOBAS 2.0. Furthermore, an integrated expression analysis of Pb-responsive miRNAs, mRNAs, metabolites and proteins were conducted based on our previous studies (Wang et al., , 2015a. The same gene annotations were used to link expression levels across the four technologies.

Protein Identification and Quantification of Radish Root in Response to Pb Stress
A total of 241,107 spectra were obtained from the iTRAQ LC−MS/MS proteomic analysis of these three group samples including an untreated control (CK), Pb500 and Pb1000. After data filtering to eliminate low-scoring spectra, a total of 17,579 unique spectra were matched to 10,896 peptides representing 9725 unique peptides. The criteria that the peptides with ≥ distinct unique peptide were used for protein identification, and a total of 3898 protein species were successfully detected at a 95% confidence limit. The relative molecular mass distribution of the identified protein species showed that the most abundant sequences ranged from 10 to 80 kDa.
Only the proteins with at least two unique identified peptides were retained for further quantification analysis. In total, 2141 protein species were permitted the quantification of their abundance. Among them, 345 protein species were identified as differentially accumulated between Pb500 and CK groups including 177 (51%) up-accumulated, and 168 (49%) downaccumulated (Table S1). Among the 509 DAPS responding to the Pb1000 treatment, 48% (244) showed increased abundance, while 52% (265) showed decreased ( Table S2). The proportion of common DAPS were calculated in response to the two Pb stress conditions. In total, 721 species were differentially accumulated under Pb treatments with 19% (135) common DAPS in both Pb treatments and 81% (586) distinct DAPS, among which 211 unique to Pb500-stress and 375 unique to Pb1000-stress (Figure 1). FIGURE 1 | The Venn diagram of differential abundance protein species (DAPS) involved in radish root under Pb500 and Pb1000 exposure.

Functional Classification and Enrichment Analysis of the DAPS
All of the 721 non-redundant DAPS sequences were functionally annotated using BLAST2GO (Conesa et al., 2005;Conesa and Götz, 2008). GO terms were assigned to query sequences and cataloged groups were produced basing on biological processes, molecular functions and cellular components. In total, 709 DAPS were assigned with 1687 GO terms and could be classified into 43 functional groups at the second level 2 (Figure 2). Among them, biological processes represented the largest category, containing 22 groups with metabolic (78.0%) and cellular process (75.3%) as the two most frequent terms. Within the molecular function category, the predominant groups were assigned to catalytic activity (56.7%) and binding (55.1%). For cellular components, all the DAPS were mostly located in cell (84.2%) and organelle (67.7%).
To systematically understand its biological functions in terms of networks, the 721 DAPS were mapped into the KEGG pathway database. A total of 75 pathways were assigned, which were largely involved in carbohydrate metabolism (such as starch and sucrose metabolism, pyruvate metabolism, glycolysis/gluconeogenesis, pentose phosphate pathway, and amino sugar and nucleotide sugar metabolism), amino acid metabolism (phenylalanine metabolism, alanine, aspartate, and glutamate metabolism, cysteine and methionine metabolism and arginine and proline metabolism) and energy metabolism (carbon fixation in photosynthetic organisms, carbon fixation pathways in prokaryotes, methane metabolism, oxidative phosphorylation and sulfur metabolism) ( Table S3). The pathways with the DAPS number larger than five were shown in Figure 3.

Characterization of the Critical Pb Stress-Responsive Proteins in Radish
Protein modification, the balance between synthesis and degradation, is a critical form of regulation that is coordinated to achieve a unified cellular under the stress of environmental stimuli (Hinkson and Elias, 2011). In this study, 27 protein species implicated in protein translation, processing and degradation were identified to be co-differentially accumulated      The abundance levels for three translation initiation family proteins (CL89.Contig1, Unigene926 and Unigene26860), two RNA binding proteins (Unigene29764 and CL7993.Contig7) and one protein disulfide isomerase (CL1707.Contig1) were observed to be increased under these two Pb exposures. However, one RNA helicase (CL1044.Contig6) was identified to be decreased responding to these two Pb treatments. All three identified ubiquitin-conjugating enzymes (CL4590.Contig3, CL83.Contig2, and Unigene282) and one cyclase family protein (Unigene2385) were coordinately increased at the Pb500 exposure while decreased at the Pb1000 (Table 1). It is known that the antioxidant enzymes could ensure cellular protection from the reactive oxygen species (ROS) mediated damage in plant responding to various environmental stresses. In the present study, the abundance changes for many antioxidant enzymes were detected including two 12-oxophytodienoate reductases (CL2689.Contig1 and CL5307.Contig5), one glutaredoxin-like protein (CL10604.Contig2), one thioredoxin (CL13300.Contig1) and eight peroxidases (CL7174.Contig1, CL1027.Contig4, Unigene8343, Unigene21004, Unigene21694, CL3864.Contig1, Unigene27310, and CL959.Contig2). Among these identified DAPS for antioxidant enzymes, most of them were decreased in radish under the exposure of Pb except one OPR (CL2689.Contig1) and one GR (CL10604.Contig2) up-accumulated at both Pb treatment, and two peroxidases (Unigene 8343 and Unigene 21694) up-accumulated under the Pb500 exposure (Table 1).

Associated Analysis of Pb-Responsive Genes through Proteomic in Combination with Other-Omic Techniques of Radish in Response to Pb Stress
In order to investigate whether the protein level is correlated with the corresponding mRNA level alterations, the proteomic data was compared with our previous transcriptome data between CK vs. Pb1000 . Proteins were considered to be correlated when quantified proteins have expression information at the transcriptome level, and a total of 1667 protein species (77.9% of all the quantified proteins) were detected from the transcriptome data ( Table S4). Out of the 510 protein species found to be significantly accumulated under Pb1000 stress of radish, there were only 57 DAPS could be matched with their cognate differentially expressed genes (DEGs) including 25 regulated in the same trends and 32 found in the opposite direction (Table S5). In addition, there were 343 DAPS exhibiting no change in mRNA expression and 144 DEGs without altered expression at protein level (Table S5). These results indicated that the accumulation of transcripts and proteins occur independently. A generally low congruency of proteomic and transcriptional profiles has been reported in other previous studies (Lan et al., 2011;Zhuang et al., 2013).
Increasing evidences have revealed that miRNA-mediated gene regulation play significant roles in plant response to HM stress, and 19 Pb-responsive miRNAs and their corresponding target mRNAs were identified in radish by siRNA sequencing and degradome analysis (Wang et al., 2015a). By exploring an integrated expression analysis of miRNA, mRNAs and proteins, a total of seven miRNA-mRNA pairs and matching proteins were identified in radish response to Pb stress, which mainly included two miRNA families, miR156 and miR396. As shown in Table 4, almost all Pb-responsive miRNAs and their corresponding mRNA targets had an anti-relationship at protein level. For example, two members of miR156 family (miR156b and miR156c) were found to be up-regulated in response to Pb exposure, while the expressions of their targets (glutaredoxin, aldose 1-epimerase and malate dehydrogenase) were all repressed at protein levels. Two corresponding protein species for the down-regulated miR396b were identified during Pb stress of radish, one (Translation initiation factor eIF-3) was up-accumulated, while the other one (thioredoxin-like 1-2, chloroplastic) was shown to be down-accumulated.
For deeply dissecting the molecular mechanism underlying Pb tolerance and homeostasis in radish, we further performed an association analysis of the seven Pb-responsive miRNA-mRNAproteins pairs with the differentially regulated metabolites identified in radish roots during Pb stress (Wang et al., 2015b). Based on the intersected pathway analysis, three Pb-responsive metabolites including glucose, citrate and malate were found to be linked with the miRNA-mRNA-proteins, which may play a significant role in radish response to Pb stress (Table S6).

DISCUSSION
Proteomic technique provides a powerful tool for the analysis of molecular mechanism of plant response against stresses, and a path for linking gene expression to cell metabolism in the rapidly processing post-genome era (Ahsan et al., 2009;Liang et al., 2013). Increasing evidences have revealed that proteomic studies are playing important roles in the post-genomic era for characterizing the molecular mechanisms underlying plant responses to HM stresses (Fukao et al., 2011;Wang et al., 2013Wang et al., , 2014Sebastiani et al., 2014). Radish is an important root vegetables worldwide, which has been considered as one of the most significant root species for dissecting the molecular regulatory network of HM stress in Brassicaceae crops (Xu et al., 2012;Wang et al., 2014Wang et al., , 2015aLiu et al., 2015;Xie et al., 2015). In the current study, comparative proteomics analysis by iTRAQ together with other-omics techniques reveals complex regulatory network and provides insights into the response of radish root under Pb exposure. To the best of our knowledge, this is the first study to systematically investigate and characterize the protein abundance changes under Pb stress exposure in radish root with iTRAQ-based proteomics analysis.

Proteins of Signal Sensing Mechanisms Involved in Radish Root Response to Pb Stress
It was known that the cell wall could activate a variety of specific stress-responsive signaling proteins when perceive outer stress conditions (Jamet et al., 2006;Day et al., 2013). In previous studies, a set of cell wall-related protein species exhibited dynamic changes in response to HM stress, such as Al exposure in rice and cadmium (Cd) exposure in flax (Hradilová et al., 2010;Wang et al., 2013Wang et al., , 2014. Currently, "cell wall" (GO: 0005618) was identified as the most enriched GO term for all the 135 common DAPS involved in the Pb500 and Pb1000 stress, which indicated that cell wall may play a vital role in response to Pb stress of radish. Two critical inhibiting protein species including a pectinesterase inhibitor (CL2232.Contig2) and a polygalacturonase inhibitor (CL2803.Contig4) were both downaccumulated during the stresses of Pb500 and Pb1000 in radish. The two protein species were verified to modulate the activities of pectinesterases (PE) and polygalacturonases (PG), which could help the plant to fine-tune cell wall remodeling processes when exposure stress conditions (Pelloux et al., 2007;Ferrari et al., 2012). The down-accumulation of the inhibiting proteins maybe result in up-accumulation of PE and PG, indicating their positive role to the cell wall remodeling when encounter the Pb stress of radish. Calreticulin (CRT) is an abundant endo reticulum Ca 2+ binding protein, and plays a critical role in Ca 2+ homeostasis and signaling sensing network (Corbett and Michalak, 2000;Jia et al., 2009). In this study, one CRT (calreticulin-3-like isoform x1, CL10110.Contig2) was found to be induced by the Pb1000, indicating that the altered expression of CRT maybe function as signaling molecular in modulating the radish plant to adapt to the Pb-stressed environments. Additionally, jasmonate (JA) is known to be a vital signaling molecule which can be activated in response to a wide range of environmental cues including heavy metals (Maksymiec et al., 2005;Schommer et al., 2008;Liu et al., 2009). In this study, the expression level of 12oxo-phytodienoic acid reductase (OPR), a key enzyme involved in the JA biosynthesis pathway, showed differential abundance levels during the Pb stress. The OPR-1 (CL2689.Contig1) was up-accumulated in radish roots when subjected to the either Pb500 or Pb1000 stress, and a similar result was found at the proteome level in rice root exposure to As stress (Ahsan et al., 2008;Srivastava et al., 2009). However, the expression level of OPR-3 (CL5307.Contig5) was decreased in radish roots under the conditions of Pb500 and Pb1000 exposure, indicating that different members of the OPR were differentially accumulated under the Pb exposure in plant.

Proteins of Carbohydrate and Energy Metabolism Involved in Radish Root Response to Pb Stress
The accumulation of HM in plants can severely affect the photosynthetic pathway and thus lead to symptoms of toxicity, such as chlorosis and growth reduction (Pourrut et al., 2013). To maintain the normal growth and development, or at least to protect the cells against excess damages, plants need to activate the carbohydrate and energy metabolism related metabolic pathways (Sarry et al., 2006;Thapa et al., 2012). In this study, "citric acid (TCA) cycle and respiratory electron transport" was identified as the most enriched pathway, and a lot of critical metabolic enzymes were shown to exhibit abundance changes under Pb exposure. Beta-glucosidase catalyzes the hydrolysis of the glycosidic and release of glucose, which then entry into the glycolysis process. ATP-citrate lyase (ACL, EC4.1.3.8) is a key enzyme involved in TCA cycle, which catalyzes the cleavage of citrate to yield acetyl CoA, oxaloacetate, ADP, and orthophosphate, and the ACL gene was found to be induced under various stresses including low temperature and drought stimulus (Hu et al., 2015). Pyruvate dehydrogenase could exert a key role in linking the glycolysis to the TCA cycle by catalyzing the formation of an acetyl-CoA from pyruvate (Vuoristo et al., 2015), which is a vital rate-limiting step reaction that determines the rate and efficiency of TCA cycle. In this study, a beta-glucosidase 21 (Unigene24884), ATP citrate lyase (ACL12533.Contig1), and a pyruvate dehydrogenase (CL6783.Contig1) were observed to be both up-accumulated during the Pb500 and Pb1000 stress of radish root, which revealed their positive potential role in producing more reducing power to compensate high-energy demand response to Pb stress. However, the abundances of the some other metabolic enzymes including Isocitrate dehydrogenase (CL4038.Contig1), fumarate hydratase (Unigene27210, Unigene20278), and malate dehydrogenase (Unigene23018) were decreased under the Pb-stresses, indicating their activity suffered repression. These decreased proteins may influence the biosynthesis and accumulation of organic acids such as citrate and malate, which play critical roles in HM tolerance (Ma et al., 2001;Wang et al., 2015b).

Proteins of Antioxidative Defense and Detoxification Involved in Radish Root Response to Pb Stress
The presence of excess metal ions causes ROS in plants, which can irreversibly damage the cells and attack macromolecules (Ahsan et al., 2009). However, the ROS can be scavenged by plant antioxidant defense system consisting of antioxidant compounds and enzymes (Apel and Hirt, 2004;Suzuki et al., 2012). The ascorbate (ASA)-glutathione (GSH) cycle is one of the main antioxidant systems in plants to keep ROS under control being, which involves a lot of critical antioxidant enzymes in a series of cyclic reactions to detoxify H 2 O 2 (de Sousa et al., 2016;Noshi et al., 2016). One of the antioxidant enzymes is peroxidase (POD), which can detoxify H 2 O 2 by oxidizing ascorbate. It was reported that different members of the peroxidase gene family were differentially accumulated under various HM exposures in Arabidopsis (Kumari et al., 2008), wheat (Houde and Diallo, 2008) and aspen (Grisel et al., 2010). In this study, the term of "peroxidase activity" (GO: 0004601) was identified as one of the most over-represented GO terms and six POD protein species were identified as DAPS, indicating the POD may be the critical ROS-scavenging protein in radish root response to Pb stress. The abundance of another antioxidant enzyme involved in ASA-GSH cycle, monodehydroascorbate reductase (MDHAR, CL11999.Contig2) was also found to be altered in response to Pb stress. The main function of MDHAR is responsible for ascorbate regeneration in plant tissues, and the expression level changed during the HM stress in B. juncea root (Alvarez et al., 2009) and A. Halleri shoot (Farinati et al., 2009). Additionally, many thiol-containing antioxidants, peroxiredoxin (PRX33, unigene27310), glutaredoxin-like protein (GRX, CL10604.Contig2) and thioredoxin (TRX, CL13300.Contig1), were also found to be altered during the Pb-stressed condition of radish root. The Prx is a thiol peroxidase with multiple functions, which was found to be induced under various HM stresses such as Cd  and As (Pandey et al., 2012). Glutaredoxin (GRX) and thioredoxin (TRX) could be oxidized by peroxides and regenerate peroxiredoxins, which were verified to play direct roles in the antioxidative system by regenerating peroxiredoxins oxidized by peroxides .

The Genetic Regulatory Network of Radish Root Response to Pb Stress
The mechanism of HM response is a complex process that a variety of genes and response components involved in plants (Fukao et al., 2011;Sebastiani et al., 2014). In the present study, to reveal the molecular mechanism underlying Pb stress response in radish root, a putative schematic network was put forward based on proteomic information of this study in conjunction with integrated analysis of miRNA, transcriptomic and metabolomic data (Figure 4). Once Pb enters the cell of radish plants, the cell wall firstly perceive stress signals and then activate specific stress-responsive signaling molecules including fine-tuning cell wall remodeling processes through special proteins like PE and PG, triggering the HM stress responsive hormones levels such as JA as well as modulating calcium-signaling molecules. With the aid of signaling molecules, the stress signals were transmitted and ultimately give rise to the alterations in gene expression and protein levels. A direct consequence of HM stress was the disturbance in the balance of protein synthesis and degradation, which is essential to both cellular homeostasis and dynamics because almost all biological processes need the involvement of enzymes. Under Pb stress, many key enzymes (i.e., ATP citrate lyase, Isocitrate dehydrogenase, fumarate hydratase and malate dehydrogenase) involved in the glycolysis and TCA cycle were severely affected, which ultimately cause alteration of some metabolites including glucose, citrate and malate. The glucose could act as the osmoprotectants protecting the cell constituents, and the organic acids (such as citrate and malate) played critical roles in chelating toxic HM ions. Meanwhile, a series of other defense responses were triggered to cope with Pb-induced injuries. For example, the ASA-GSH cycle was the main antioxidant systems for scavenging the accumulated ROS to alleviate oxidative damages. Moreover, several Pbdefense protein species (glutaredoxin, aldose 1-epimerase malate dehydrogenase and thioredoxin) encoding genes targeted by miRNA156 and miR396, which were found to play critical FIGURE 4 | The putative genetic regulatory network of Pb stress response in radish root. 12-oxo-phytodienoic acid reductase-1 (OPR-1), 12-oxo-phytodienoic acid reductase-3 (OPR-3), Ascorbate (ASA)-Glutathione (GSH), Calreticulin-3-like isoform x1 (CRT3), Glutaredoxin (GRX), thioredoxin (TRX), Leucine-rich repeat family protein (LRR), Monodehydroascorbate reductase (MDHAR), Pectinesterases (PE), Peroxidase (POD), Peroxiredoxin (PRX), Polygalacturonases (PG), Protein disulfide isomerase (PDI), Translation initiation factor (TIF). roles in defense system when radish root responds to Pb stress.
In summary, this is the first report on systematically investigating and characterizing the protein abundance changes under Pb stress in radish root with the iTRAQ technique. A total of 3898 protein species were successfully detected and 2141 were quantified. A subset of 721 protein species were differentially accumulated upon at least one Pb exposure, and 135 ones showed common abundance alterations under two Pb-stressed conditions. GO and pathway enrichment analysis revealed that these 135 common DAPS were strongly enriched in the categories of cell structure, carbohydrate and energy metabolism-related pathways and antioxidative defense. Furthermore, the integrative analysis of transcriptomic, miRNA, degradome and metabolomic with proteomic data provided a strengthened understanding of radish root response to Pb stress, and a schematic genetic regulatory network was put forward. The genes associated with signal sensing, protein synthesis and degradation, carbohydrate and energy metabolism and ASA-GSH cycle for ROS scavenging as well as several key miRNAs and metabolites were crucially responsible for radish in response to Pb stress. Overall, the results of this study would be beneficial for further dissecting molecular mechanisms underlying other plant responses to HM stresses, and facilitate genetically effective management of HM contamination in root vegetable crops.

AUTHOR CONTRIBUTIONS
YW designed the experiments and drafted the manuscript. LX, MT and WZ participated in the design of the study and performed the data analysis. MT, RW and HJ planted radish seedlings and collected samples. LX and WC reviewed and edited the manuscript. LL conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.  Table S1 | The DAPS involved in radish root response under Pb500 exposure.