iTRAQ-Based Quantitative Proteome Analysis Insights into Cold Stress of Winter Turnip Rapa (Brassica rapa L.) Grown in the Field

: Winter Turnip rapa ( Brassica rapa L.) is a major oilseed crop in Northern China, where its production was severely affected by chilling and freezing stress. Previous studies have demonstrated that differentially accumulated proteins (DAPs) were expressed in roots and leaves under control chilling stress. In this study, the isobaric tag for relative and absolute quantification (iTRAQ) technology was performed to identify DAPs under freezing stress. Two winter rapaseed varieties, Longyou 7 (cold-tolerant) and Lenox (cold-sensitive), were used to investigated morphological, physiological, cell and protein levels in the shoot apical meristem (SAM) of field-grown Brassica rapa to reveal the molecular mechanisms of cold stress tolerance. Compared to Lenox, Longyou 7 had a lower SAM of height, higher collar diameter. The level of malondialdehyde (MDA), SAM of height and IAA content were repressed. At the same time, compared to Lenox, the soluble sugars (SS) content, superoxide dismutase (SOD) activity, peroxidase(POD)activity, soluble protein (SP) content and collar diameter increased in Longyou 7. In total, we identified 6330 proteins, among this, 98 DAPs were expressed in L7 CK/Le CK, 107 DAPs were expressed in L7 d /Le d 183 DAPs were expressed in Le d /Le CK, 111 DAPs were expressed in L7 d /L7 CK. Quantitative real-time PCR (RT-qPCR) analysis of the coding genes for seventeen randomly selected DAPs were performed for validation. These DAPs were identified from the two winter rapa seed cultivars involved in the biological process, cellular component and molecular function analysis, which revealed glutathione transferase activity, carbohydrate-binding, and glutathione binding, glutathione metabolic process and response IAA were closely associated with the cold stress response. Some cold-induced proteins, such as glutathione S-transferase phi 2(GSTF2), might play essential roles cold acclimation in Brassica rapa of SAM. Our work will help to provide valuable information for responding to the cold stress in Brassica rapa L. iTRAQ-Based quantitative proteome analysis. Previous RNA-seq and iTRAQ-based quantitative proteome analysis applied on leaves and roots tissues were limited to obtaining relative differences of the genes or proteins expression under controlled environment. Our study had demonstrated that DAPs used to evaluate relative differences of protein expression as cold acclimation progressed under field conditions from early autumn to winter. The field trials were designed to capture the big picture of global changes of the proteomics during cold acclimation in Longyou 7 and Lenox. iTRAQ-based quantitative proteome analysis was performed in Longyou 7 and Lenox. The analysis of protein was performed using 8-plex iTRAQ-based comparative proteome analysis. Three independent biological replicates were simultaneously performed. After using Mascot to search against the Brassica-rapa 1.0 database. In total, 6 328 protein were identified (Supplementary Table S1). increase of osmotic adjustment substance content could reduce the freezing point to increase plant cell membrane stabilization and protect membrane integrity during dehydration caused by freezing treatment in B. rapa . The aboveground tissues of winter rapa seed were damaged by low temperatures. Longyou 7 has strong protective enzyme activity and high content of soluble regulators that can remove ROS. These data showed that osmometric and ROS were adjusted in B. rapa to adapt to the cold stress.


Introduction
Cold stress, including chilling (0-15℃) and freezing (<0 ℃) temperatures, dramatically affects plant growth, yield, and quality of crop species. It is one of the major threats to crop production in northwestern China [1,2]. Plants adopted several strategies to develop complicated and effective coldresponsive mechanisms to survive cold stress, such as the level of chaperones and antioxidants variations in their leaf tissue structure, maintaining osmotic balance by altering membrane structure, and the activation of cold-related genes [3,4]. Some plants abandon their cold sensitive structure to protect themselves in winter, such as above-ground or avoid freezing injury by shrinking vegetative organs into underground organs. Although the way in which plants perceive and respond to temperature cues is not well understood, wheat winter survival involves two important evolutionary adaptive mechanisms: cold acclimation and vegetative/ reproductive transition [5,6]. The most important growth site in the nutritional compared with reproductive transformation is the SAM, it play an important role in plant growth. Meristem function is essential for local auxin synthesis, and auxin biosynthesis genes also show specific as well [7]. So, cold stress of molecular mechanisms will shed light on the regulatory mechanisms for cold acclimation in Brassica rapa and provide an effective approach to select targeted candidate genes for manipulation and and/or cross-breeding of agronomic [8]. Our study will shed light on regulatory mechanisms of cold acclimation and provide an effective approach to select targeted candidate genes for manipulation and cross-breeding of agronomic in B rapa.
iTRAQ is a gel free mass spectrometry technique for the corresponding protein. It applies isobaric amine specific tags to compare peptide intensities between samples, followed by inferring quantitative values for the corresponding proteins [9]. iTRAQ labelling combined with two-dimensional liquid chromatography-tandem mass spectrometry(2D-LC-MS compared with MS) is a high through put quantitative proteomics technique widely used in identifying and quantifying cellular metabolic changes in proteome-related analysis [10]. This technique has been successfully applied to identify plant protein responses to stress in various plants such as Vitis amurensis [11], Triticum aestivum [12], Gossypium hirsutum L [13], Lycopersicon esculentum Mill [14], Brassica rapa L [15], Spica Prunellae [16]. Differential proteomics analysis focused on screening and identifying differences and changes in different species or growth stages, revealing and validating proteomics changes.
Brassica winter rapa (Brassica rapa L.) is an oil crop that can survive the winter in the cold and arid region of northern China and is one of the effective surface cover crops in winter and spring. It has significant ecological and economic benefits in agricultural production [17]. However, the extremely low temperature in northern China negatively influences the regional distribution and safe production of winter rapa in winter. Previous studies have shown that strong cold-resistant plants show that their leaves crawl and grow close to the ground in morphology. Before winter, organic matter is preferentially distributed to the root, storing enough organic matter and establishing an extensive root system, guaranteeing its safe overwintering and vegetativeness construction after winter. Previously, our research team has used high-throughput sequencing technology and successfully identified cold-stressresponsive microRNAs(miRNAs) [18], differentially expressed genes(DEGs) [19] and differentially accumulated proteins(DAPs) [20] in roots and leaves of B. rapa. The SAM depression and bulge in the seedling stage might be associated with their strong cold resistance. The functions of SAM target protein regulations have not yet been understood. The objectives were to explore the relationship between SAM-target proteins regulations and their relationship to cold stress tolerance. The proteins identified can be used to inform breeding programs and improve our understanding of the molecular mechanisms underlying cold resistance.

Results
Physiological development and cold tolerance. Longyou 7 and Lenox was used to study the morphologyical and physiological differences at the seedling stage under cold stress. Lenox had an upright petiole with several leaves growing in the leaf positions. However, Longyou 7 was always opposite (Fig 1A). When growing in an open field, the height of SAM levels was measured at the different growth stage. Lenox of average height of SAM since the first sampling date and displayed a high apical meristem at the next stage. Compared with the first sampling date (CK), Lenox of average height of SAM was significantly increased in the the next stage, Longyou 7 of average height of SAM was increased in the next stage. Compared with Longyou 7 of average height of SAM, Lenox of average height of SAM was significantly higher. Compared with the first sampling date (CK), Longyou 7 of collar diameter was significantly increased in the the next stage (d), Longyou 7 of collar diameter was increased in the next stage. Compared with Longyou 7 of collar diameter, Lenox of collar diameter was significantly higher (Fig 1B). Compared with the L7 d and Le d, the level of SOD, POD activity and soluble protein content significantly decreased in the L7 CK and Le CK samples. The MDA content significantly increased in L7 CK and Le CK. Under freezing treatment (0℃ and -11℃), the MDA level of Le CK and Le d was significantly higher than that of L7 CK and L7 d. Still, other physiological indexes of Le CK and Le d were lower than those of L7 CK and L7 d (Fig 1C ). These results suggested that compared with Lenox, Longyou 7 had lower SAM but higher tolerance to freezing stress. Strong coldresistant varieties such as Longyou 7 were found have a hollow SAM in the samples collected from field experiments, which were 5-10 cm below the soil surface. Thus the SAM always located in the moist soil under a relatively stable temperature, which is beneficial for winter B. rapa to the overwinter. These findings showed the SAM seedling stage have been related to cold tolerance. iTRAQ-Based quantitative proteome analysis. Previous RNA-seq and iTRAQ-based quantitative proteome analysis applied on leaves and roots tissues were limited to obtaining relative differences of the genes or proteins expression under controlled environment. Our study had demonstrated that DAPs used to evaluate relative differences of protein expression as cold acclimation progressed under field conditions from early autumn to winter. The field trials were designed to capture the big picture of global changes of the proteomics during cold acclimation in Longyou 7 and Lenox. iTRAQ-based quantitative proteome analysis was performed in Longyou 7 and Lenox. The analysis of protein was performed using 8-plex iTRAQ-based comparative proteome analysis. Three independent biological replicates were simultaneously performed. After using Mascot to search against the Brassica-rapa1.0 database. In total, 6 328 protein were identified (Supplementary Table S1). iTRAQ quantitative identified DAPs. Based on fold change (FC) ＞1.5 (p＜0.05), DAPs were up accumulated, FC＜0.8 (p ＜ 0.05), DAPs were down accumulated. Here, up or down accumulated proteins were determined using the Longyou 7 and Lenox. The Venn diagram reflected the DAPs in two winter rapa seed varieties. 76 DAPs were found in between the Le d /Le CK and L7 d /L7 CK and 48 DAPs were found in both the L7 d /Le d and L7 CK /Le CK (Figure 2A and 2B). Of these, there were 35 DAPs unique to L7 d /L7 CK. This was consistent with the results of the principal component analysis (PCA) of the three biological replicates, which indicated that the three biological replicates of each sample had a good repeatability, L7 d /L7 CK as well as Le d /Le CK had a great difference ( Figure 2C ). As shown in Figure2D   To further characterize the functions of the DAPs, KEGG pathway mapping was also performed. Only significantly enriched categories (p < 0.05) were selected (Supplementary Table S7).

Discussion
Proteomic analysis has identified many differentially accumulated proteins mainly related to plant secondary metabolism in germinated seeds under cold stress. In this study, a large number of cold stress-responsive proteins were identified in winter turnip rapa. Previously several candidate genes has been identified by EST analysis which may be involved in the cold stress response in winter turnip rapa. However, the mechanisms underlying the effect of cold stress on SAM are largely unknown. Therefore, that physiological and proteomic analyses were performed in this study to examine the mechanism of the increased freezing tolerance in winter tunip rapa.

Freezing Stress Affects Morphology, Physiological and Biochemical Changes in Brassica rapa.
Freezing stress (<0°C) harmfully affects plant growth and development, limits their geographic distribution, and significantly reduces agronomic productivity [20]. In the natural cooling process, plants change in morphology, physiology and biochemistry due to low temperature stress, especially freezing stress [21]. The main target of freezing injury is cell membranes. Lipid, peroxide and MDA could reflect the damage of plants under low temperature [22]. The regulation of osmotic adjustment substances such as soluble protein, and free proline could keep the stability of cellular structure and osmotic balance [23]. Meanwhile, proline acts as a cell structure maintainer, a signaling molecule and a ROS scavenger in the response to cold stress, among these, SOD and POD are key factors of antioxidant enzymes under stress [24,25]. There was an increase in the content of MDA and SP in the CK compared to freezing stress. However, the content of Pro, SOD and POD activity decreased in in the CK compared to freezing stress ( Figure 1C). The increase of osmotic adjustment substance content could reduce the freezing point to increase plant cell membrane stabilization and protect membrane integrity during dehydration caused by freezing treatment in B. rapa. The aboveground tissues of winter rapa seed were damaged by low temperatures. Longyou 7 has strong protective enzyme activity and high content of soluble regulators that can remove ROS. These data showed that osmometric and ROS were adjusted in B. rapa to adapt to the cold stress.
DAPs involved in lipid metabolism. Several studies have shown that high concentration of unsaturated fatty acids are commonly presented in plants grown under cold stress, thus cold-tolerant plants often have higher unsaturated membrane lipids [26,27]. Linoleic an alpha-Linolenic acids play key roles in maintaining membrane integrity and fluidity, and have been demonstrated to contribute to the osmotic tolerance in rice [28,29]. Unsaturated sphingolipids content is reported to be higher in coldtolerant species [28,30]. Glycerolipid metabolism is one of the metabolic pathways related to fatty acid degradation. Glycerol plays vital roles in cell energy generation and lipid synthesis [31]. In previous study, lipoxygenase LOX could damage to cell membranes and other cell components [32]. LOX2 was downregulated under cold stress, Linolenate hydroperoxide lyase (CYP74B2) and LOX2 play a role in JA biosynthesis [33,34]. Cold hardiness increased in S. cerevisiae when acohol dehydrogenase class-3(ADH3) was overexpressed [35]. In this study, DAPs were found involve in Linoleic acid metabolism, alpha-Linolenic acid metabolism, glycolipid metabolism, sphingolipid metabolism, and fatty acid degradation. Of these, only fatty acid degradation was up regulated and rest were down regulated in L7 d /Le d. ADH3(Bra033706) was up regulated, LOX2(Bra004057, Bra003526, Bra003526) and CYP74B2(Bra012788) were down regulated in L7 d /Le d. It is known that ADH3 might play an important role in the cold stress of B. rapa.
DAPs involved in carbohydrate and energy metabolism. Carbohydrate and energy metabolism are the basis of plant growth, development and morphogenesis [36]. Carbohydrate metabolism occupies the core of energy metabolism since it provides the essential saccharides and energy that plants need [37]. It was reported that beta-amylase 5 (BAM5) might also play a important role in the petunia response to cold stress [38]. Down regulation of alpha-galactosidase(α-Gal) has been shown to increase plant freezing tolerance [39]. The pentose phosphate pathway(PPP) supplies NADPH, an electron donor [40]. Probable 6-phosphogluconolactonase 4(PGL4), an enzyme participating in PPP agreed with the increased needs for ROS scavenging [39]. 1,4-alpha-glucan-branching enzyme2-41(SBE2.1) increased in salt and drought stressed mulberry roots and leaves, are involved in starch biosynthetic metabolism [40]. In this study, α-Gal 1 (Bra028660) and PGL 4(Bra009759) were down regulated in L7 d /Le d BAM5(Bra038088) and SBE2.1 (Bra005269) were up regulated in L7 d /Le d. BAM5 and SBE2.1 were involved in the starch and sucrose metabolism. Previous study suggested that many starch grains accumulated in the leaf cells under stress. So, this might provide a new direction for future research in the B. rapa under cold stress.
Energy metabolism may play important roles in the development of chilling resistance in plant [41]. Nitrilase 2 (NIT 2) promotes nitrilase biosynthesis, and nitrilase leads to increased IAA content in Arabidopsis under cold stress [42]. Cytochrome oxidase subunit 6b-1(COX 6 B) was downregulated in rice under salt stress. Chlorophyll a-b binding protein increased in mulberry leaves under salt-drought stress conditions [43]. In our study, DAPs participated in energy metabolism, including nitrogen metabolism, photosynthesis-antenna proteins, oxidative phosphorylation were identified in L7 d/Le d. Our results showed that the expression levels of chlorophyll a-b binding protein CP 24 (Bra026745) and NIT 2 (Bra035006) increased, while and COX6B (Bra031391) decreased in L7 d / Le d. Our results also indicated that IAA content of the Longyou 7 higher compared to Lenox in the SAM of the B. rapa. Therefore the focus of our future research will be on NIT 2 and auxin in the B. rapa under cold stress.
DAPs involved in amino acid metabolism. It is reported that freezing stress response amino acids metabolic pathway in Eriobotrya japonica and Populus tomentosa under cold or freezing conditions [44,45]. Our study demonstrated that several amino acid metabolic pathways were identified after freezing stress. Seventeen DAPs were annotated inclduing Glutathione metabolism, Cysteine and methionine metabolism, Alanine, aspartate and glutamate metabolism, Arginine biosynthesis, Tryptophan metabolism, Glycine, serine and threonine metabolism and Tyrosine metabolism. Glutathione Stransferases (GSTs) are a major family of detoxification enzymes that are important in protecting plants against oxidative damages. Auxin induced the expression of several glutathione S-transferases (GSTs) such as AtGSTF 2, so GSTF 2 is up regulated might be linked to stress-mediated growth responses [46]. Our results showed that bigger of SAM indicated relatively high levels of IAA content, while a relatively smaller SAM indicates relatively low levels of IAA content. so auxin content and thus was related auxin content in relatively high of SAM. Our results showed that five GSTs (Bra000875, Bra021673, Bra000876, Bra036259, Bra012409) and glutathione S-transferase DHAR 1 (Bra025749) were down regulated, while GSTU (Bra012422, Bra000474) and GSTs (Bra032010) were up regulated in L7 d /Le d. In addition, antioxidant enzymes played important roles in the removal of excess ROS to reduce oxidative stress by GSTs. the activities of a number of antioxidant enzymes, such as SOD and POD can remove ROS. The increased accumulation of these proteins indicated that plant cells initiated their antioxidant mechanisms to maintain redox homeostasis and resist cold stresses. DAPs involved in translation and biosynthesis of secondary metabolite. Posttranslational modification (PTM) in Brassica juncea cold stress systems include S-nitrosylation [46]. Post-translational modification that confers the glutathione S-transferase (GST) under oxidative stress [47]. The biosynthesis of secondary metabolites are associated with their involvement in tolerance by cold stress in plant [48]. Secondary metabolites are involved in plant responses to abiotic stresses and provide a valuable contribution to the antioxidant activity of plant tissues [49]. Protein disulfide isomerase (PDI) gene was over-expressed that can enhance heat stress [50]. Tropinone reductase (TR1) significantly increased the low temperature in Arabidopsis [51]. ADH3(Bra033706) was up regulated and PDI(Bra020239) was down regulated in Longyou-7 d / Lenox d. Tropin one reductase homolog (Bra039974) was involved in biosynthesis of other secondary metabolites and folding, sorting and degradation. GST was involved in posttranslational modification in Longyou-7 d /Lenox d. Our study showed its great potentials in improving the genetics of low temperature tolerance in plants.

Conclusions
In this study, two winter cultivars (Longyou 7 and Lenox) with different phenotypes for cold stress were used to analyze molecular mechanisms for cold response. Base on proteomics data, major proteins involved in metabolic pathways through accumulating differentially accumulated proteins in both winter turnip rapa varieties under cold stress. DAPs were identified at treatment between Longyou 7 and Lenox.
In addition, based on the functional analysis, we concluded that DAPs involved in amino acid metabolism, carbohydrate and energy metabolism, lipid metabolism, translation and biosynthesis of secondary metabolites, which showed that freezing stress greatly changed the survival of winter turnip rapa. In summary, these findings serves as our understanding of the molecular mechanisms in freezing tolerance and increases protein resource for rapa seed for freezing tolerance breeding in winter turnip rapa.

Methods
Plant materials, field, sample collection. Two different winter turnip rapa cultivars were used in this study.
Longyou 7 and Lenox from Gansu Agricultural University (Gansu, China). Both were grown in normal agronomic field trials at Gansu Research Center of Rapaseed Engineering and Technology in Lanzhou, China. All method was performed in accordance with the relevant guidelines and regulations. Three replicates were prepared for each cultivar in the field. The field-collected tissue were used periodically for cold hardiness and morphological measurements form October 13 to December 16, when daily average temperature had reached 0℃. The daily average temperature is the mean value of the lowest temperature in 15 days. Each sampling point, plants were collected in the field between 10:30 a.m. and 11:30 a.m. to minimize circadian rhythm effect. Data was collected daily at average daily minimum temperature (Tmin) from October to December in 2018. The Tmin was 0℃ On October 13, December 16 was -11℃. The sample for analysis were named, Le CK (sample of "Lenox" at 0℃ October 13), Le d (SAM of "Lenox" at -11℃), L7 CK (SAM of "Longyou 7" at 0 ℃), and L7 d (SAM of "Longyou 7" at -11℃). To provide sufficient tissue for RNA extraction, SAM (less than a 1 cm long section at the base of the crown) was collected from four to six plants in each point and replicate. Samples were flash frozen in liquid nitrogen and stored in a -80℃ freezer.
Morphological and physiological analyses. The samples were performed three times with similar results. POD activity [51], CAT activity [52], SS content [53], IAA content [54], SP and MDA content [55] were analyzed as indicators of physiological response. The majuscules indicated a significant difference (p < 0.01) for the freezing stress-treated samples compared with chilling stressed samples. Values are means and SD of three biological replicates, each calculated from the mean of three technical replicates. Protein extraction. iTRAQ analysis were sent to Shanghai Luming biology (Shanghai, China). Three biological replicates were used for iTRAQ-based comparative proteomics analysis. Approximately 500 mg from each biological replicate was ground into powder using liquid nitrogen and dissolved (vortex blending) with 500 Μl extraction buffer. The samples were ground at 60 Hz for 2 min. Then supplemented with extraction buffer for 1 mL and mixed and added with Tris-phenol buffer and mixed for 30 min at 4 °C. The homogenates were centrifuged at 7100×g for 10 min at 4 °C. The supernatant was collected, added for five volumes of 0.1 M cold ammonium acetate methanol buffer and stored at −20 °C overnight. The homogenates were centrifuged at 12,000×g for 10 min to collect precipitations. The methanol was replaced with acetone to remove the methanol, and the steps used to obtain precipitate were repeated. The samples were centrifuged at 12000×g for 10 min at 4°C to collect supernatants. The supernatants were dried at room temperature for 5min and dissolved in sample lysate form 2 h. The supernatants were centrifuged again to remove precipitations completely. The protein concentration was quantified using BCA method and the protein purity was detected by SDS-PAGE.15 μg proteins of each sample were separated on 12% SDS-PAGE gel. Protein digestion, iTRAQ labeling and RP chromatography separation. Protein digestion was performed according to the FASP procedure [56]. IAA was added to the final concentration of 50 mM in the dark for 40 min. The protein solutions were centrifuged on the filters at 12,000×g for 20 min at 4 °C. Remove the supernatant and add TEAB to the solutions and centrifuged at 12,000×g for 20 min. The solutions were collected and lyophilized. The lyophilized samples were suspended in TEAB (100 μL, 50 mM) and 40 μL of each sample was transferred into new tubes for labeling. Each sample add iTRAQ label reagent (iTRAQ® Reagents-8plex kit, Sigma) following the manufacturer's protocol (Applied Biosystems, Foster City, CA, USA). All labeled peptides were pooled together. iTRAQ labeled peptides were fractionated by RP chromatography separation using the 1100 HPLC System(Agilent). RP separation was performed on the Agilent Zorbax Extend RP column (5 μm, 150 mm × 2.1 mm). Tryptic peptides were separated at an eluent flow rate of 300μL • min −1 and monitored at 210 and 280 nm. Dried samples were harvested from 8 min to 50 min and elution buffer were collected in every minute and numbered from 1 to 10 with pipeline. The separated peptides were lyophilized for MS detection [57]. Mass spectrometry analysis. All LC-MS compared with MS analyses were performed on a Q-Exactive mass spectrometer (Thermo, USA) equipped with a Nanospray Flex source (Thermo, USA). The peptides mixtures were loaded by a capillary C18 trap column (3 cm × 100μm, C18, 3μm, 150 A) and separated by a C18 column (15 cm × 75μm, C18, 3μm, 120 A) on an ChromXP Eksigent system (AB Sciex). Full MS scans were acquired in the mass range of 300-1600 m compared with a mass resolution of 70,000 and the AGC target value was set at 1,000,000. The 10 most in-tense peaks in MS were fragmented with higher-energy collisional dissociation (HCD) with collision energy of 30. MS compared with MS spectra were obtained with a resolution of 17,500 with an AGC target of 200,000 and a max injection time of 50ms. The Q-E dynamic exclusion was set for 15.0 s and run under positive mode [58]. Protein identification and function annotation. Raw data of iTRAQ-labeled proteins by was search against rapa genome protein database in National Center for Biotechnology Information (NCBI) using the Proteome Discoverer TM 2.2 (Thermo, USA). Database searches were performed with trypsin digestion specificity, and the cysteine alkylation was considered as parameters in the database searching. For protein quantification method, iTRAQ 8-plex was chosen. For protein identification, a decoy database search approach was used to determine the false discovery rate (FDR) with acceptance if their FDR < 1.0% while protein identification containing at least two peptides. The molecular functions of the identified proteins were classified according to their gene ontology annotations and their biological functions. Only the proteins identified with at least two different peptides and P value< 0.05, and quantified with a ratio of fold change > 1.5 or fold change < 0.8 were considered. The NCBI(https://www.ncbi.nlm.nih.gov/) and Uniprot databases (https://www.uniprot.org/) were chosen to the validation and annotation of the protein sequences. Gene Ontology (GO) annotation for the identified proteins was assigned according to Uniprot databas [59]. Statistical analysis. The control and treatment groups were analyzed for statistical significance of differences between multiple groups using one-way ANOVA followed by Duncan's multiple comparisons test. All calculations were performed using SPSS software (version 19.0, IBM, Armonk, NY, USA). All results are presented as mean ± SD from three independent biological replications. Treatment means were separated by the Duncan multiple range test at P ＜ 0.01. We use minmax normalization method through the R programming language to analysis transcriptional and proteomic represent expression values of heat map. RNA extraction and qPCR analysis of gene expression. There were three biological replicates, and three technical replicates were performed for each gene. The total RNA SAM were completed by using TRNzol Universal Reagent (Tiangen, China), according to the manufacturer's instructions. One microgram of total RNA was used for first-strand cDNA synthesis according to the protocol supplied with PrimeScript™ RT Master Mix (TaKaRa Biotechnology Dalian, China) and qRT-PCR amplification reactions were performed using a LightCycler®96 Real-Time PCR System (Roche, Basel, Switzerland). Primers used for qRT-PCR assay. The primer sequences for internal standard. The relative gene expression was calculated using the 2−∆∆Ct method [60]. The primer sequences were designed by