Comparative proteomic analysis provides insight into a complex regulatory network of taproot formation in radish (Raphanus sativus L.)

The fleshy taproot of radish is an important storage organ determining its yield and quality. Taproot thickening is a complex developmental process in radish. However, the molecular mechanisms governing this process remain unclear at the proteome level. In this study, a comparative proteomic analysis was performed to analyze the proteome changes at three developmental stages of taproot thickening using iTRAQ approach. In total, 1862 differentially expressed proteins (DEPs) were identified from 6342 high-confidence proteins, among which 256 up-regulated proteins displayed overlapped accumulation in S1 (pre-cortex splitting stage) vs. S2 (cortex splitting stage) and S1 vs. S3 (expanding stage) pairs, whereas 122 up-regulated proteins displayed overlapped accumulation in S1 vs. S3 and S2 vs. S3 pairs. Gene Ontology (GO) and pathway enrichment analysis showed that these DEPs were mainly involved in several processes such as “starch and sucrose metabolism”, “plant hormone signal transduction”, and “biosynthesis of secondary metabolites”. A high concordance existed between iTRAQ and RT-qPCR at the mRNA expression levels. Furthermore, association analysis showed that 187, 181, and 96 DEPs were matched with their corresponding differentially expressed genes (DEGs) in S1 vs. S2, S1 vs. S3, and S2 vs. S3 comparison, respectively. Notably, several functional proteins including cell division cycle 5-like protein (CDC5), expansin B1 (EXPB1), and xyloglucan endotransglucosylase/hydrolase protein 24 (XTH24) were responsible for cell division and expansion during radish taproot thickening process. These results could facilitate a better understanding of the molecular mechanism underlying taproot thickening, and provide valuable information for the identification of critical genes/proteins responsible for taproot thickening in root vegetable crops.


Introduction
Radish (Raphanus sativus L., 2n = 2× = 18), one of the most important worldwide root vegetable crops, is an annual or biennial herb belonging to Brassicaceae family. The fleshy taproot is the significant edible part in radish plants, which contains abundant nutrient substances such as carbohydrates, crude fiber, vitamin C, protein, and other secondary metabolites including glucosinolate, and it determines the final yield and quality 1 . In addition, it is of significant value on diet and medicine [2][3][4][5] . Therefore, it is urgent to clarify the molecular mechanism underlying taproot thickening in radish plants.
In the past years, the morphological, physiological, and anatomical characterization of taproot thickening have been extensively studied on radish 6,7 . The morphogenetic process displays dynamic changes in the period of taproot thickening, which is determined by the interactions of genetic, environmental, and physiological factors 8,9 . Essentially, fleshy taproot development is the result of related genes programmed expression 8,9 . In recent years, with the rapidly developed "omics" technology, the draft genome sequences and transcriptome studies of R. sativus have been reported, which provide a valuable database for identification of the critical genes and genetic manipulation in radish [10][11][12][13] . Using RNA-Seq technique, characterization of transcriptome and miRNA may dissect the molecular mechanism underlying taproot thickening, with several miRNAs and differentially expressed genes (DEGs) identified during three different stages of taproot thickening (pre-cortex splitting stage, cortex splitting stage, and expanding stage) 8,14 . However, the molecular mechanism underlying taproot thickening in radish has not been comprehensively uncovered at the proteome level.
Proteomics studies provide a powerful tool for exploring related proteins in specific tissues with different stages in the post-genomic era, therefore comparative proteomic analysis for identifying proteins involved in radish taproot thickening would be an indispensable way for complementing genomics analysis to further explore molecular mechanism governing radish taproot thickening. Currently, the two-dimensional gel electrophoresis  in combination with mass spectrometry (MS) has been employed for proteomic studies providing useful tools for protein separation and quantification [15][16][17][18] . However, it has some limitations in protein abundance identification. More recently, isobaric tags for relative and absolute quantitation (iTRAQ), a new and powerful technology, was used in quantitative proteomics by isotope tagging and high-performance liquid chromatography (HPLC) 19,20 . High sensitivity (less than 1 ppm), low detection limit (proteins of less than 10 kD or greater than 200 kD) in identification and quantitation of proteins aspects make iTRAQ a valuable approach for proteomics studies in some plant species including cassava 21 , potato 22,23 , and Brassica napus 24 .
Previous studies reported that cortex splitting stage is a sign of transition from primary growth to secondary growth, and it is of significance to investigate the expression changes of proteins during the pre-cortex splitting stage, cortex splitting stage, and expanding stage of taproot thickening in radish 9,11 . In this study, comparative proteomic analysis of monitoring differentially expressed proteins (DEPs) from three libraries of precortex splitting stage (S1, 10 days after sowing (DAS)), cortex splitting stage (S2, 20 DAS), and expanding stage (S3, 40 DAS) were conducted by iTRAQ-coupled LC-MS/ MS. Furthermore, association analysis between DEGs and DEPs in three development stages of taproot thickening was performed for further screening of critical genes involved in taproot formation in radish. The outcomes of this study could be beneficial for further dissection of the molecular mechanism governing radish taproot thickening, and provide a fundamental basis for genetic improvement of taproot formation in root vegetable crops.

Plant materials
Seeds of an advanced inbred line of radish (R. sativus L.) "NAU-YH" were germinated on moist filter paper in dark at 25°C for 3 days, and then cultured in plastic pots with 16 h light (25°C) and 8 h dark (18°C) for the rest of the experimental period. At pre-cortex splitting stage (S1, 10 DAS), cortex splitting stage (S2, 20 DAS), and expanding stage (S3, 40 DAS), taproots were sampled, respectively, and three biological replicates were collected for each stage. Equal amounts of taproot samples from three independent biological replicates of each stage were pooled and immediately frozen in liquid nitrogen, then stored at −80°C for proteomic analysis.

Protein preparation
The taproot samples were ground in liquid nitrogen to a fine powder and extracted in lysis buffer, then 1 mM phenylmethyl sulfonyl fluoride (PMSF), 2 mM ethylene diamine tetraacetic acid (EDTA), and 10 mM dithiothreitol (DTT) were added. The samples were sonicated for 15 min and centrifuged at 25,000×g for 20 min. The supernatant was mixed well with 5× volume of chilled acetone and incubated at −20°C for 2 h. After centrifugation at 16,000×g for 20 min, the supernatant was discarded and then repeated previous processes once. The protein concentrations were measured with the Bradford assay 25 , and 12% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) was used for separating proteins and checking the proteins.

iTRAQ labeling and SCX fractionation
The extracted total protein from each sample solution (100 μg) was digested with Trypsin Gold with a ratio of protein:trypsin = 20:1 (w/w) at 37°C for 4 h, and then Trypsin Gold (protein:trypsin = 20:1) added once more at 37°C for 8 h. After trypsin digestion, peptides were dried by vacuum centrifugation and reconstituted in 0.5 M tetraethyl-ammonium bromide (TEAB) with the next steps following the manufacturer's procedure of 8-plex iTRAQ reagent (Applied Biosystems). Briefly, 1 U of iTRAQ reagent was thawed and reconstituted in 70 μL isopropanol. Samples from pre-cortex splitting stage S1, cortex splitting stage S2, and expanding stage S3 were labeled with isobaric tags 116, 117, and 118, respectively, and incubated at room temperature for 2 h, and then pooled and dried by vacuum centrifugation [26][27][28] .
Then the labeled peptides were separated by SCX chromatography using separation column of 4.6 × 250 mm (Ultremex SCX column) with LC-20AB liquid system. In detail, the mixed labeled peptides were reconstituted in 4 mL buffer A (25 mM NaH 2 PO 4 in 25% ACN, pH 2.7), and gradient-eluted in buffer B at a flow rate of 1 mL/min in column using 5% buffer B (25 mM NaH 2 PO 4 , 1 M KCl in 25% ACN, pH 2.7), 7 min; 5-60% in 20 min; 60-100% and maintenance for 1 min; 5% maintenance for 10 min 29 . The eluted process was supervised under 214 nm absorption photometry. The columns were cleaned to eliminate salts using StrataX. The extracted liquid was lyophilized and stored at −80°C.

LC-ESI-MS/MS analysis
The peptide was dissolved in buffer C (5% ACN, 0.1% FA) and centrifuged at 20,000×g for 10 min. The supernatant (approximately 10 μL) was loaded on a trap column and then the peptides were eluted onto an analytical column by LC-20AD nanoHPLC (Shimadzu, Kyoto, Japan) according to the manufacturer's instruction. In detail, the peptides were loaded at 8 μL/min for 4 min, then the gradient program was performed as follows: it was started from 2 to 35% in buffer D (95% ACN, 0.1% FA) at 300 μL/min, followed by a 5 min linear gradient to 60%, then a 2 min linear gradient to 80% and maintenance for 4 min, and finally return to 5% in 1 min 29 . Data were acquired by a TripleTOF 5600 System. The MS was operated with more than 30,000 resolutions. A sweeping collision energy was set for 35 ± 5 eV with a dynamic exclusion setting of 1/2 of peak width (15 s), ensuring that the same precursor ion was fragmented not more than twice.

Protein identification
For peptide data analysis, raw data was employed with the Mascot search engine (Matrix Science, London, UK; version 2.3.02) against the Raphanus_sativus (224,406 sequences) protein sequence database (ftp://ftp. kazusa.or.jp/pub/radish/). For protein identification, a series of standard parameters were set as follows: 0.1 Da fragmented mass tolerance, 0.05 Da peptide mass tolerances and one max missed cleavage. To reduce the probability of false peptide identification, only peptides at the 95% probability were used for protein identification by a Mascot probability analysis, and each confident protein contains at least one unique peptide 28 . For protein quantitation, protein has at least two unique spectra, and the parameter "median" was considered as the standard of quantitative protein ratios 26,30 . The identification results of significant DEPs were filtered with P-value < 0.05 and fold change (FC) > 1.2 31 .

Bioinformatics analysis
A hierarchical cluster analysis of the DEPs was carried out using Cluster 3.0 software. The similarity of proteins was calculated with Euclidean distance, and the averagelinkage method was selected for clustering 32 . Functional annotations of the DEPs were employed by KOBAS 2.0 program search against the non-redundant (nr) protein database deposited in NCBI. The database of Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), and Cluster of Orthologous Groups (COG) were performed to categorize and group the DEPs 31,33 .

RT-qPCR analysis
Quantitative real-time PCR (RT-qPCR) was employed to validate the quality of iTRAQ results. Total RNAs from taproot samples including S1, S2, and S3 stages was extracted and reverse transcribed to cDNA following the manufacturer's instructions (Tiangen Biotech Co., Ltd., China). Each reaction was carried out using 10 μL 2× SYBR green reaction mix, 2.0 μL diluted cDNA, and 0.2 μM of each primer in a total volume of 20 μL system. RT-qPCR amplification reactions were conducted on a Real-Time PCR Detection System (Bio-Rad iQ5, USA) following the reported protocol 34 . The specific primers used for RT-qPCR were designed with Beacon Designer 7.0 software (Premier Bio-soft International, USA), which are listed in Supplementary Table S1. Three replicates for each gene assay were performed and mRNA expression levels were normalized by Rs-Actin gene. The relative gene expression value was calculated with the 2 −ΔΔCT method.

Protein identification
Based on LC-ESI-MS/MS analysis, a total of 411,811 spectra were generated in this study, which included 95,124 matched spectra and 46,912 unique spectra. Totally 27,080 peptides containing 16,722 unique peptides as well as 6342 proteins were identified. Moreover, protein mass was dominantly enriched in 20-60 kDa, and the peptide number and peptide length were mainly distributed at the range of 1-10 and 6-29 (in amino acids), respectively (Supplemental Figure S1).

Functional annotation of proteins
GO function classification showed that the proteins were mainly annotated to the terms of cellular process (14.77%, in biological process ontology), cell part (22.82%, in cell component ontology), and binding (43.40%, in molecular function ontology) (Fig. 1a-c).
KEGG pathway analysis revealed two dominant pathways: "Metabolic pathways" (ko01100) and "Biosynthesis of secondary metabolites" (ko01110). In addition, "Starch and sucrose metabolism" (ko00500) pathway was the seventh top of all, which was one of the most important pathways involved in taproot thickening (Fig. 1d) 13 .
The COG database was employed for protein orthologous classification, and a total of 5260 proteins were aligned to 24 COG terms. These results showed that the category of "general functions prediction only" (17.60%) was the largest, followed by "Posttranslational modification, protein turnover, chaperones" (12.36%), "Translation, ribosomal structure and biogenesis" (9.13%), "Carbohydrate transport and metabolism" (7.89%), and "Energy production and conversion" (7.13%). In contrast, the categories of "Cell motility" and "Nuclear structure" were less than five proteins with high homology (Supplemental Figure S2).
To further identify the DEPs involved in taproot thickening in radish, two comparison pairs were overlapped. For instance, we overlapped between S1 vs. S2 and S1 vs. S3 pairs for identifying proteins related to taproot thickening initiation, and overlapped between S1 vs. S2 and S2 vs. S3 pairs for identifying proteins related to taproot thickening process. More interestingly, totally 256 up-regulated proteins including 77 unknown proteins displayed overlapped accumulation in S1 vs. S2 and S1 vs. S3 pairs, whereas 122 up-regulated proteins including 43 unknown proteins displayed overlapped accumulation in S1 vs. S3 and S2 vs. S3 pairs (Supplementary Table S3). Our findings suggest that these proteins may play critical roles in taproot thickening initiation and taproot formation in radish.

GO and pathway enrichment of DEPs
We used GO and pathway enrichment analysis of DEPs to determine significantly enriched DEPs compared with the background of all DEPs at P-value < 0.05 as the threshold value, and thereby explored the main biological functions, biochemistry metabolism, and signal transduction pathways.
In the S1 vs. S2 comparison, GO enrichment analysis showed that 1154, 1005, and 1090 DEPs were assigned to 271, 582, and 1708 GO terms of cell component, molecular function, and biology process, respectively. Among these, 83, 110, and 377 GO terms were significantly enriched in cell component, molecular function, and biology process, respectively (Supplementary Table S4). In the S1 vs. S3 comparison, 993, 864, and 944 DEPs were annotated to 254, 549, and 1689 GO terms of cell component, molecular function, and biology process, respectively. Of these, 76, 125, and 408 GO terms were significantly enriched in cell component, molecular function, and biology process, respectively (Supplementary Table S4). In the S2 vs. S3 comparison, 882, 742, and 833 DEPs were assigned to 250, 513, and 1601 GO terms of cell component, molecular function, and biology process, respectively. Of all, 65, 76, and 219 GO terms were significantly enriched in cell component, molecular function, and biology process, respectively (Supplementary Table S4). In addition, the shared significantly enriched GO terms of three libraries (Cluster frequency ≥50%) were calculated, and cytoplasm, cytoplasmic part, catalytic activity, metabolic process, cellular metabolic process, response to stimulus, and single-organism metabolic process were obtained.

Associated analysis of mRNA and proteins during taproot thickening in radish
Proteome and transcriptome analysis showed two difference levels, which reflect genes expression. In order to mutual corroborate the data reliability, association analysis of DEPs data with previous transcriptome data was performed in this study (Supplementary Table S7) 14 . The  (Fig. 3g-i). The result may suggest that gene regulation mechanisms were different in mRNA and protein level, and the associated proteins/genes were critical for exploring the molecular mechanism of taproot thickening in radish. Furthermore, the cluster analysis was performed to identify the characteristic of correlation between transcriptome and proteome (Fig. 4).

RT-qPCR validation
To evaluate the validity of iTRAQ results, 17, 11, and 9 proteins from three pairwise comparisons (S1 vs. S2, S1 vs. S3, and S2 vs.  Table S1). As shown in Fig. 5, the transcript levels of BG1, PPC, INV, CYP79F1, LOX1, GLS1, GLN, PGDH, and EXPB1 were higher in the pre-cortex splitting stage than in the cortex splitting stage and expanding stage, whereas SNX1 and PDC2 were higher in the cortex splitting stage and expanding stage than in the pre-cortex splitting stage at mRNA levels. Interestingly, the transcript levels of FP6, SEX4, HSC70, and GLN were higher in the expanding stage than in the pre-cortex splitting stage and cortex splitting stage, while the transcript of CDC5 was higher in cortex splitting stage than in the expanding stage. Furthermore, the mRNA levels of signaling transduction proteins including MAPK3, MAPK5, and AUR were higher in the pre-cortex splitting stage than in the cortex splitting stage. Overall, it was found that there was a good concordance between iTRAQ and mRNA expression levels (Fig. 5), indicating the reliability of iTRAQ-based quantitative proteomic data analysis.

Discussion
As an important storage organ, radish taproot thickening determines its final yield and quality of the product. miRNA and mRNA of radish taproot thickening at precortex splitting stage (S1), cortex splitting stage (S2), and expanding stage (S3) have been studied, which provided novel insights into the genetic regulatory network during taproot thickening 9,14 . However, the molecular mechanism underlying taproot thickening in radish is far from being fully clarified. In this study, iTRAQ-based proteomic approach was employed to monitor DEPs from three different stages of radish taproot thickening (S1, S2, and S3). Furthermore, based on the integrative analysis of transcriptomic and proteomic data, an overview of proteome changes involved in radish taproot thickening at three developmental stages was put forward (Fig. 6). To the best of our knowledge, this is the first comprehensive investigation to characterize the potential functional proteins involved in taproot thickening in radish.

Proteomic studies could be a supplement to transcriptome analysis
In this study, iTRAQ-based proteomic approach was employed to identify DEPs from three different stages of radish taproot thickening (S1, S2, and S3), as a follow-up to previous study 9,14 , which could be beneficial to further Fig. 4 Cluster diagram for the concordance between proteome and transcriptome based on pairwise comparison among the three libraries. a S1 vs. S2 pair. b S1 vs. S3 pair. c S2 vs. S3 pair dissecting the genetic molecular mechanism underlying radish taproot thickening. Furthermore, for further screening critical genes involved in taproot formation in radish, association analysis between DEGs and DEPs in three development stages of taproot thickening was also performed in this study. The relatively weak correlation of genes at mRNA and protein level indicates that an alteration in a proportion of transcripts may not be translated into changes in protein abundance. Low congruency of genes in mRNA and protein level was similar with other previous studies 20, 35 .
As shown in Fig. 6, the blue font indicated the proteins with the same expression patterns in mRNA and protein level, whereas the normal font represented unique proteins. For example, proteins involved in cell wall loosening and reconstruction including XTH24 and XTH32 and proteins involved in signal transduction such as PP2C58, BSK8, and BSL1 were unique proteins that were only identified in this study but not gathered in transcriptome analysis. The results showed that the transcriptome and proteome levels were various under taproot thickening in radish, and it could be presumed that proteomic studies was used as an effective supplement to transcriptome analysis.

Phytohormone is a basic regulatory factor for taproot thickening
Phytohormone is critical for the growth modulation of the root system 36 . Five main hormones including auxin, cytokinin, ethylene, abscisic acid (ABA), and gibberellic acid (GA) were extensively investigated in plant root system. Auxin-to-cytokinins ratio regulates root and shoot meristem. In general, the high ratio is favorable for the root meristem and control of root architecture, while the lower ratio is conducive to shoot meristem. The auxin influx carriers including AUX1, PGP, LAX, and auxin efflux transporters such as ABC transporters affect the distribution of auxin and regulate root morphogenesis 37 . Previous studies showed that changing the level of auxin influences cytokinin concentration, which functioned by the mechanism of relevant proteins controlling hormone accumulation 38 . Several other hormones including ethylene, ABA, and GA are also responsible for root development. For example, ethylene inhibits adventitious root formation in Arabidopsis 39 , ABA could modulate main root elongation in response to drought, and GA could enhance cell elongation of root. In radish plants, cytokinin plays critical role in controlling developmental process of taproot thickening 7 . In this study, "plant hormone signal transduction" is one of the most enriched pathways, and several DEPs including Gns7 (glucan endo-1,3-beta-glucosidase 7, Rsa1.0_00048.1_g00023.1), AHP2 (histidine-containing phosphotransfer protein 2, Rsa1.0_00568.1_g00001.1), PYR1 (abscisic acid receptor, CL14734.Contig1_CKA, CL3244.Contig2_NAU-LB), and PP2C58 (protein phosphatase 2C 58, CL6992.Contig1_NAU-LB) were upregulated in this process, suggesting that these proteins would be likely involved in the regulation of radish taproot thickening (Supplementary Table S6). However, a system network of hormone interactions in radish taproot thickening needs to be further investigated.
Interestingly, the basic formation mechanism of the root and tuber crops (RTCs) could be interlinked, but the genes and proteins involved in these processes may be diverse. For example, sporamin including sporamin A and sporamin B were specifically identified to be up-regulated from storage roots in sweet potato 17 , and patatin proteins (patatin protein 13 and patatin protein 15) were specifically identified from potato 16 . The sporamin and patatin protein were not identified from taproot thickening in radish, whereas several proteins including SUS1, EXPL2, XTH32, and CDC5 were specifically identified to be significantly differentially expressed in radish, suggesting that there is certain specificity in different types of RTCs ( Fig. 6; Supplementary  Table S2). However, carbohydrate metabolism was one of the most important pathways in the development of RTCs, and the proteins involved in this process including fructokinase, phosphoglycerate kinase, NADH dehydrogenase, glyceraldehyde-3-phosphate dehydrogenase, fructosebisphosphate aldolase, ADP-glucose pyrophosphorylase, glucose-1-phosphate adenylyltransferase, and alcohol dehydrogenase were identified to be differentially expressed both in potato 16 and radish (Supplementary Table S2). Meanwhile, several proteins including protein disulfide isomerase and anionic peroxidase were identified to be significantly expressed both in storage root of sweet potato 17 and radish taproot thickening (Supplementary Table S2), indicating that the protein expression level in the energy supply was relatively conservative. Therefore, we speculated that the possible reason is the difference in Fig. 6 An overview of proteome changes involved in radish taproot thickening at three developmental stages. The blue font represents overlap proteins between DEPs and DEGs, and the normal font represents unique proteins in this study. MAPK4 mitogen-activated protein kinase 4, CAM putative calmodulin, SRP54 signal recognition particle subunit, PP2C58 probable protein phosphatase 2C 58, SPP1 sucrose-phosphatase 1, PFK3 6-phosphofructokinase 3, PFK5 6-phosphofructokinase 5, bglL beta-glucosidase L, UGD1 UDP-glucose 1-dehydrogenase, SUS1 sucrose synthase 1, PME51 phosphoglucomutase 1, bZIP59 bZIP transcription factor 59, EF1BB elongation factor 1B alpha-subunit 2, PP2C59 protein phosphatase 2C 59, Gns7 glucan endo-1,3-beta-glucosidase 7, MAPK6 mitogen-activated protein kinase 6, BSK8 brassionsteroid-signaling kinase 8, BSL1 serine/threonineprotein phosphatase, EXPL2 expansin-like A2, EXPB1 expansion B1, XTH24/XTH32 xyloglucan endotransglucosylase 24/32, FRUCT5 betafructofuranosidase, bglH beta-glucosidase H, PGM1 phosphoglucomutase 1, APL3 glucose-1-phosphate adenylyltransferase, bglI beta-glucosidase I, CDC5 cell division cycle 5-like protein proteins to be involved resulting in diverse accumulation of substances in potato, sweet potato, and radish.

Cell division and cell expansion determining taproot thickening
Taproot thickening is a process characterized by changes in substances and energy. During this period, the basic metabolism of the taproot was reduced, and it would be act as a storage bank. A series of proteins are involved in the biosynthesis and metabolism of sugar, starch, and protein, which play critical roles in the initial step of taproot thickening in radish 13,40,41 . In the present study, cell division cycle 5-like protein (CDC5), a DNA binding protein belonging to a member of MYB3R-and R2R3type, was up-regulated in S1 vs. S2 pair, consistent with its action in Arabidopsis 42-45 . Interestingly, CDC5 gene was also identified to be up-regulated in S1 vs. S2 pair 14 . Meanwhile, cell division and cell enlargement tend to be along with the cell wall loosening and reconstruction, and expansin (EXP) and xyloglucan endotransglucosylase/ hydrolase (XTH) proteins played critical roles in these processes. Previous studies showed that down-regulation of EXPB1 gene enhanced storage root development in sweet potato 46 , and XTH24 gene involved in cell wall elongation in Arabidopsis 47 . In this study, the expansin proteins including EXPB1 (expansin B1) and EXLA2 (expansin-like A2) were identified to be differentially expressed, among them EXLA2 was all up-regulated in S1 vs. S2 pair and S1 vs. S3 pair, whereas EXPB1 was downregulated in S1 vs. S2, all in line with gene expression pattern at mRNA level ( Fig. 6; Table 2) 14 . Moreover, XTH32 was up-regulated in S2 vs. S3 pair, while XTH24 that involved in GA-mediated signaling pathway, cell wall loosening, organization or biogenesis, and xyloglucan metabolic process was down-regulated in S1 vs. S2 pair and S1 vs. S3 pair ( Fig. 6; Table 2). The results suggested that CDC5, EXPB1, and XTH24 might play critical roles in taproot thickening process in radish.
In addition, previous studies illustrated that transcription factors (TFs) including MADS-box, ABF/AREB, and homeobox were also responsible for the formation of modified stems or storage roots 38 . In the present study, MYB and bZIP TFs were both up-regulated in S1 vs. S2 pair, suggesting that MYB and bZIP might play critical roles in regulating the development of taproot thickening in radish, which were consistent with their action of rhizome formation in lotus root 14,42 .
In conclusion, the comparative proteome changes among three developmental stages of taproot thickening were systematically investigated in radish. A total of 1862 DEPs were identified during radish taproot thickening. GO and pathway enrichment analysis showed that several DEPs were mainly involved in "plant hormone signal transduction", "starch and sucrose metabolism" and "biosynthesis of secondary metabolites". Furthermore, the integrative analysis of DEGs and DEPs data enhanced our Table 2 continued understanding of taproot thickening molecular mechanisms in radish. Overall, it was concluded that taproot thickening initiation was triggered by phytohormone, and several functional proteins including CDC5, EXPB1, and XTH24 could contribute to cell division and expansion during radish taproot thickening. Together, these findings would provide fundamental insights into comprehensive clarification of molecular regulatory network of taproot thickening in radish.