iTRAQ-based quantitative proteomic analysis of tea plant (Camellia sinensis (L.) O. Kuntze) during cold acclimation and de-acclimation procedures

During the life cycles of plants, high salinity, desiccation and low temperatures are the major environmental stresses. These stresses impair the development and productivity of crops, horticultural plants and wild species, and even threaten their survival. For the perennial plants grown in frigid and temperate zones, low temperature in winter is a serious stress to overcome. Based on long term evolution, plants have developed complicated mechanisms to perceive and respond to low temperature stresses. In higher plants, cold acclimation (CA) is an important mechanism to defend against cold temperatures in winter^[1]. Many physiological, biochemical, and structural changes happen during plant CA process^{[2, 3]}. The increasing number of transcriptome analyses of gene expression at the RNA transcript level have helped us to better understand the molecular basis of CA^{[4, 5]}. However, the transcript level of a given RNA does not always strictly correlate with the corresponding protein level in plant cells. And little has been done to elucidate the protein abundance changes during CA and de-acclimation (DA) processes^{[6, 7]}. As a valuable approach to study the stress responses at the level of protein abundance, proteome analysis has been implemented in many studies. In Arabidopsis, putative plasma membrane proteins associated with CA were identified using a mass spectrometric approach^[8]. Based on the two-dimensional difference gel electrophoresis (2-D DIGE) analysis, 26 differently expressed proteins were identified in rice, and the cellular phospholipase Dα1 protein was proven as a key candidate involved in the CA signaling pathway^[9]. Balbuena et al.^[10] found that the tolerant lines of sunflower showed a higher number of differentially expressed proteins in leaves, compared with freezing susceptible lines. As an important semi-permeable cellular membrane, plasma membrane plays vital roles in response to abiotic stress such as low temperature stress in plant. And plasma membrane proteins may change during CA^[11]. Although many cold-stress-related proteins have been identified and analyzed, and much knowledge about CA has been added recently, the understanding of the CA mechanism is still limited, especially in woody and evergreen plants.

Tea plant [Camellia sinensis (L.) O. Kuntze] is a woody, perennial, evergreen plant and is widely planted in developing counties of the tropics and sub-tropics as an important cash crop^[12]. Low temperature in winter is a key environmental factor restricting the growth of tea plants, which could lead to damage to tea plantations, decline of production, and even plant death. Therefore, understanding the tea plant CA mechanism and functional genes, and application in tea plant breeding is a crucial way to improve tea plant cold tolerance. Similar to many other plants, huge changes happen at cellular, physiological and metabolic levels during the tea plant CA process, such as the relative electrical conductivity, concentration of malondialdehyde and relative water content decrease^{[13, 14]}. Oppositely, the palisade tissue thickness increased and the plasma membrane stability enhanced through the increasing of total proteins and unsaturated fatty acids^[15]. The content of soluble sugars also increased in winter^[14]. Additionally, cold induced or related genes were identified and functions were validated by different technologies in tea plant^[16−20]. To further highlight the mechanisms of CA, we performed a transcriptome analysis based on RNA-seq. Many differentially expressed genes were identified and confirmed using quantitative RT-PCR analysis. These genes were grouped into signal transduction genes, cold-responsive transcription factor genes, plasma membrane stabilization related genes, osmosensing-responsive genes and detoxification genes, etc^[4]. The transcriptome analysis provided a valuable chance to look into global gene expression changes at RNA level during the CA process in tea plant. Transcriptomic data are not often consistent with protein or metabolism data due to post-transcriptional modifications. So to provide new ideas towards the CA mechanism in tea plant, we examined the proteome during CA and DA processes in tea plant. Tea plant, being a broad leaved woody evergreen, may provide novel information on cold resistance mechanisms of other broad leaved evergreen plants in winter.

In the present study, leaf samples at the non-acclimated (NA), fully acclimated (CA) and de-acclimated (DA) stages were collected according to our previous work^[4] and analyzed using isobaric tags for relative and absolute quantitation (iTRAQ) quantitative proteomic approach following the workflow shown in Supplemental Fig. S1. Finally, over 1,300 differentially expressed and functioning in varied biological processes proteins were identified.

The tea plants of C. sinensis (L.) O. Kuntze 'Longjing 43' grown in the field of the Tea Research Institute, Chinese Academy of Agricultural Sciences (TRI, CAAS) (N 30°10', E 120°5'), were used in this study. Intact mature leaves at NA, CA and DA stages were collected for further study according to Wang et al.^[4]. Each biological replicate contained ten intact leaves collected from ten individual plant and three biological replicates of each stage were collected. Those collected samples were frozen in liquid N₂, and stored at −80 °C for protein extraction and iTRAQ assay and qRT-PCR analysis.

Protein extraction and gel eletrophoresis

Protein extraction was performed according to the method of Yang et al.^[21] with minor revision. Frozen leaves were ground to a fine power and weighted 1.0 g, and then 5 ml lysis buffer (7 M Urea, 2 M Thiourea, 4% CHAPS, 40 mM Tris-HCl, pH 8.5) was added. After that, 1 mM and 2 mM PMSF and EDTA were added respectively. Five minutes later, DTT was added with a final concentration of 10 mM. Then the above suspension was sonicated (200 Watts) for 15 min and centrifuged (30,000× g) at 4 °C for 15 min. Then those new supernatant was transferred into another tube and five-fold 10% chilled TCA acetone was added and incubated overnight at −20 °C. The precipitate was washed three times with chilled acetone for at least 30 min, and harvested by a centrifugation at 4 °C, 30,000× g for 15 min after each washing. The pellet was air dried before being dissolved in 500 μl 0.5 M TEAB, followed by sonication at 200 Watts for 15 min. Sonicated solution was centrifuged at 4 °C, 30,000× g for 15 min. Finally, supernatant was transferred to a new tube and quantified using BSA as standard protein with GE Healthcare's 2-D Quant Kit (Code No. 80-6483-56) following the instructions. After the quantification, SDS-PAGE and standard colloidal staining were performed to measure the quality of the protein sample.

In gel trypsin digestion, iTRAQ labeling and LC-MS/MS analyses

One hundred micrograms of protein were digested using Trypsin Gold (Promega, Madison, WI, USA) with a ration of protein : trypsin = 30 : 1, and incubated at 37 °C for 16 h. Then, the digested proteins were dried using vacuum centrifugation and reconstituted in 0.5 M TEAB (Triethylammonium bicarbonate buffer). The reconstituted NA, CA and DA samples were labeled with different isobaric tags according to the munufacturer's protocol for 8-plex iTRAQ (Applied Biosystems).

The labeled peptides were separated using Strong Cation Exchange Choematography (SCX), and analyzed using liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS)^[10]. Shimadzu LC-20AB HPLC Pump system was used for SCX chromatography as follows: the digested peptides were dissolved with 4 mL buffer A (25 mM NaH₂PO₄ in 25% ACN, pH 2.7) and then were loaded onto a Ultremex SCX column (4.6 mm × 250 mm) containing 5 μm particles (Phenomenex). Parameters of peptide separating procedures were setted as follows: first buffer A for 10 min, then 5%−35% buffer B (25 mM NaH₂PO₄, 1 M KCl in 25% ACN, PH 2.7) for 40 min, and thirdly 35%−80% buffer B for 1 min, with a flow rate 1 mL/min. The absorbance at 214 nm was monitored and separated samples were collected every 1 min. Finally, each sample were separated into 20 fractions, then each fraction was desalted using Strata XC18 column (Phenomenex) and dried by vacuum.

The desalted and dried fractions were re-suspended using buffer A (2% ACN, 0.1% FA) and centrifuged (20,000× g) for 10 min. Then, 8 μl of 0.5 μg·μl⁻¹ supernatant was injected into a Shimadzu LC-20AD nano HPLC system with an analytical C18 column (inner diameter 75 μm). The flow rate was as follows: 8 μL·min⁻¹ for 4 min, 300 μL·min⁻¹ for 40 min with 2% to 35% B (95% ACN, 0.1% FA), followed by 5 min linear gradient to 80%, and mainted for 4 min with 80% B, finally returned to 5% within 1 min.

After that, the samples were nanoelectrospray ionized (1.6 kV) and put into an Q EXACTIVE (Thermo Fisher Scientific, San Jose, CA, USA) coupled online with HPLC system. And the high-energy collision dissociation (HCD) operating mode with a normalized collision energy setting of 27.0 was applied for peptides selected for MS/MS. An Orbitrap with a resolution of 70,000 automatic gain control (AGC) was applied for detection and spectra optimization. For AGC target, the parameters of MS and MS 2 were setted as 3e⁶ and 1e⁵ respectively. While for MS and MS 2 scans, the m/z parameters were setted as 350 to 2,000 Da and 100−1,800 Da respectively. The 15 most abundant precursor ions which have a threshold ion count above 20,000 and 15 s dynamic exclusion duration time were identified based on a data-dependent procedure.

Protein identification, bioinfomatics and differential abundance protein analysis

Proteome Discoverer 1.2 (PD 1.2, Thermo Fisher Scientific) was used for the conversion of raw data files acquired from the Orbitrap and the converted MGF file was used for protein identification using Mascot 2.3.02 (www.matrixscience.com) software (Matrix Science, London, UK). Tea plant unigene (from NCBI) translation database, theaceae_txid27065 database in NCBInr, and sequences from the tea plant transcriptome studies implemented by Wang et al.^[4] were chosen as databases. The mass tolerance was set as '± 0.05 Da' and '± 0.1 Da' for intact peptide and fragmented ion identification respectively and one missed cleavage was allowed in trypsin digests. The charge states of peptides were set to +2 and +3. Deamidated (NQ), Gln->pyro-Glu (N-term Q) and Oxidation (M) were deemed as potential variable modifications. Moreover, Carbamidomethyl (C), iTRAQ8plex (K) and iTRAQ8plex (N-term) were defined as fixed modifications. A random sequence of database and the real database were used for raw spectra test in the decoy checkbox. An automatic decoy database search was conducted based on the decoy checkbox in Mascot. To improve the accuracy of peptide identification, only those peptides with significance scores (≥ 20) at the 95% confidence interval and detected as greater than 'identity' by Mascot probability analysis could be counted as identified. At least one unique peptide should be identified for each confident protein^[22].

Protein amount was estimated by spectral counts according to Balbuena et al.^[10] A pairwise comparison was performed for NA, CA and DA and 1.5 fold cutoff wit p-value < 0.05 was defined as up-accumulated or down-accumulated proteins.

Functional classification, protein abundance profiling and association analysis

All identified proteins were classified using standard Gene Ontology (GO) online tool (www.geneontology.org) analysis. The enriched GO terms among the comparisons were identified using the statistical method described by Zheng & Wang^[23]. The KEGG pathways analysis was carried out by sequence alignment against the Kyoto Encyclopedia of Genes and Genomes database^[24] using BLASTP algorithm (E-value threshold 10⁻⁵). Differentially accumulated protein species among the different samples were grouped by Cluster 3.0 software and the output files were read by javaTreeview.

To detect the correlation between protein level and corresponding transcript level at different CA stages, all the identified protein species were matched to the corresponding transcripts in the transcriptome database^[4], and then correlation analyses were implemented according to Zheng et al.^[25].

Quantitative real-time PCR (qRT-PCR) analysis

RNAprep pure Plant Kit (Tiangen, Beijing, China) was used for total RNA extraction from the samples that were used for protein extraction according to the manufacturer's instructions. RNA concentration and integrity estimation, cDNA synthesis and real-time PCR were performed according to previous descriptions^[4], and the polypyrimidine tract-binding protein gene of the tea plant (CsPTB) was used as an internal reference^[26]. The gene-specific primers for qRT-PCR were designed according to the corresponding coding sequences in the genome database^[4], and all the primer sequences were provided in Supplemental Table S1.

During the tea plant CA procedure, its relative electrical conductivity following subsequent freezing decreased, then its cold tolerance could be improved. When tea plant was de-acclimated, its relative electrical conductivity increased and its cold tolerance became weaker^[4]. We collected tea plant leaves from three different stages, NA, CA and DA for iTRAQ assaies.

In total, 2,573 unique peptides out of 2,751 peptides were harvested. Based on Mascot searching in tea plant unigene translation database, NCBInr, theaceae_txid27065 database, and 51,940 sequences database 1,331 proteins were identified. Peptide length was generally 7 to 16 amino acids, and protein mass distribution was concentrated mainly at 10 to 40 kDa. Approximately 25% proteins had 5%−10% coverage by peptides, and 712 of 1,331 identified proteins were only represented by a single peptide (Supplemental Table S2). Above analyses suggested that high quality protein abundance libraries with low redundancy were constructed successfully.

Protein quantification revealed 407 differential accumulated proteins between CA compared with NA, and among those proteins 202 were up-accumulated, while 205 were down-accumulated. Compared with CA, 115 up-accumulated and 136 down-accumulated proteins were detected in DA. In addition, compared with NA, 477 differentially accumulated proteins, including 253 up-accumulated proteins and 224 down-accumulated proteins, were identified in DA (Fig. 1). The distribution analysis of differential abundance proteins indicated that, when comparing CA with NA samples, generally down-accumulated proteins had greater abundance differences. When comparing DA with NA samples, both up-accumulated proteins and down-accumulated proteins had large abundance difference (Supplemental Fig. S2).

Figure 1. 

Number of differentially accumulated proteins among different samples.

GO enrichment analysis of differential abundance protein species

Based on the GO enrichment analysis, those differently enriched GO terms (P-value < 0.05) among the comparisons between NA, CA and DA were listed in Table 1. Most of the significant differently enriched GO terms were grouped in biological process in the comparisons. Moreover, in cellular component category, more enriched GO terms were detected in DA vs CA and DA vs NA than in CA vs NA. Similar numbers of enriched GO term were grouped in molecular function category among the three different comparisons. Briefly, in the comparison of CA vs NA, extracellular region, plastid stroma, ammonia ligase activity, acid-ammonia (or amide) ligase activity, ribosome biogenesis and glutamine metabolic process were the most highly enriched; in DA vs CA, cell wall, plastid part, adenylyltransferase activity, racemase and epimerase activity, electron transport chain, and oxidation reduction were primarily enriched; in DA vs NA, extracellular region, cell wall, binding, protein kinase activity, electron transport chain and oxidation reduction were mainly enriched. Furthermore, comparing all the enriched GO terms in CA vs NA with DA vs CA, only one common enriched GO term, plastid stroma, was found. Six terms were found when comparing CA vs NA and DA vs NA, namely acid-ammonia (or amide) ligase activity, extracellular region, ammonia ligase activity, ribosome biogenesis, binding, glutamine metabolic process. Fourteen terms were found in the comparison of DA vs CA and DA vs NA, including cell wall, oxidation reduction, nucleoside phosphate metabolic process, electron transport chain, photosynthetic electron transport chain, nucleobase, nucleoside and nucleotide metabolic process, nucleotide metabolic process, photosynthesis, light reaction, purine nucleotide metabolic process, oxidoreduction coenzyme metabolic process, pyridine nucleotide metabolic process, nicotinamide nucleotide metabolic process, photosynthesis, and generation of precursor metabolites and energy.

KEGG pathway mapping

Pathway analysis is an important approach to expose the crucial biochemical metabolism and signal transduction pathways including given proteins^[27]. The identified protein species in this study were annotated based on KEGG database. Generally, more differential abundance protein species were annotated and assigned to larger number of pathways in CA vs NA and DA vs CA, compared with DA vs CA. Moreover, the differential abundance protein species in three different comparisons were mainly mapped onto carbon fixation in photosynthetic organisms, metabolic pathway, ribosome, starch and sucrose metabolism, biosynthesis of secondary metabolites and microbial metabolism in diverse environments, protein processing in endoplasmic reticulum, photosynthesis, plant-pathogen interaction and oxidative phosphorylation. Interestingly, the pathways of glycolysis/gluconeogenesis, starch and sucrose metabolism and pyruvate metabolism related to glycometabolism were the largest proportion of differential abundance protein species in CA vs NA and DA vs NA. In addition, lysosome, glutathione metabolism, peroxisome and ascorbate and aldarate metabolism and phagosome pathways had dramatic difference in the number of differential abundance protein species in the three comparisons. As listed in Table 2 and Supplemental Tables S3 & S4, many pathways only had detectable changes in one or two specific samples, for example, proteasome, fatty acid metabolism, streptomycin biosynthesis, biosynthesis of unsaturated fatty acids and so on.

Expression profile of differential abundance protein species

The differential abundance protein species in the three different comparisons were sorted on the basis of the abundance patterns during different CA stages (Fig. 2). The patterns of differential abundance protein species varied widely and were clustered roughly into 11 groups. Among these clusters, cluster C only contained three protein species, namely pertin acetylesterase family protein, macrophage migration inhibitory factor family protein and chloroplast nucleoid DNA binding protein. These protein species were down-accumulated in three comparisons, especially in DA vs NA and DA vs CA. Furthermore, the cluster G and H were formed only by FERONIA receptor-like kinase and ribosomal protein respectively. They were dramatically down-accumulated and up-accumulated respectively in DA vs NA and CA vs NA. However, their abundance had very small change in DA vs CA. Except above minor clusters, the other eight clusters constituted at least six protein species. The biological functions of the protein species involved in these large clusters were classified based on GO enrichment analysis (Supplemental Table S5). Interestingly, most of differential abundance protein species functioned in or composed cell part, cell, organelle, catalytic activity, metabolic process, binding, cellular process, organelle part and response to stimulus. Cluster A and cluster D were formed by a small number of protein species. The protein species in cluster A were up-accumulated in DA vs NA and DA vs CA, were down-accumulated in CA vs NA, while the protein species in cluster D had an opposite profile. Similarly, cluster F and I also had an opposite abundance profile. The protein species abundances in cluster F were increased in DA vs CA and decreased in both CA vs NA and DA vs NA. Cluster B and E were respectively consisted of 36 and 49 protein species that showed reduced abundance in three comparisons, especially for the protein species in cluster E among the comparisons of CA vs NA and DA vs NA. Cluster J and neighboring cluster K had similar abundance profiles because the protein species involved in these clusters were increased at different levels among the three comparisons.

Figure 2. 

Hierachical clustering of proteins showing different abundance profiles across different samples. The data were transferred using log2.

Among these differential abundance protein species, some stress-related or responding protein species were identified. These protein species included ribosomal proteins (RPs), photosynthesis related proteins, energy metabolism proteins, osmosensing-responsiveness proteins, antioxidation-related proteins, and some signal transduction, transporter and post-translationally modified proteins, etc. These protein species were grouped into different clusters. Most of these protein species were up-accumulated in CA and/or DA stage compared to NA, indicated that complex proteomics changes were happened during CA procedure in tea plant.

Analysis of the association between proteome and transcriptome

To illuminate the potential correlation between proteome and the corresponding transcriptome, all the identified proteins were correlated to corresponding transcriptome first, and thereafter association analyses were carried out between identified proteins and corresponding differentially expressed genes among the comparisons between different CA stages. Results showed that 1310 proteins out of total 1331 identified proteins were successfully associated with transcripts (Supplemental Fig. S3). Unexpectedly, the association analysis results showed that the identified proteins had low correlation coefficients (r) with the cognate genes level in the comparisons of CA vs NA, DA vs CA and DA vs NA in transcriptome analyses, with r values of 0.0023, 0.1519 and 0.1120 respectively (Fig. 3). Because the value of correlation coefficient was close to zero, the protein level was poorly correlated with transcript levels. To further show the details about the expression patterns of identified proteins and its corresponding associated gene, clustering analyses of expression patterns was implemented (Supplemental Fig. S4). The clustering results showed that the identified protein species and differentially expressed genes were mainly grouped into three kinds of cluster. First is positive correlation, such as the cluster F in CA vs NA, cluster C in CA vs DA and cluster D in DA vs NA. Second is confused correlation, which includes both week negative and positive correlation, such as cluster G in CA vs NA, cluster D in DA vs CA and cluster E in DA vs NA. These correlation relationships took a large share. Third is clear negative correlation, such as the cluster B, C, D, E and I in CA vs NA, cluster E in DA vs CA, cluster C, F, G in DA vs NA. Main proteins included in these groups have functions that include beta-primeverosidase, glycyl-tRNA synthetase, alpha-glucan water dikinase, AT-HF, phosphoenolpyruvate carboxylase, carbonic anhydrase, hydrolase family protein and plastocyanin-like domain-containing protein.

Figure 3. 

Correlation analysis of transcript (log2 FPKM value) and protein (log2 iTRAQ value) among different samples. (a) CA vs NA; (b) CA vs DA; (c) DA vs NA.

Transcriptional expression analysis by qRT-PCR

In order to validate the correlationship between proteome and the corresponding transcriptome, 20 protein species were chosen for qRT-PCR analysis, including 17 differential abundance protein species and three unchanged abundance protein species. The expression pattern analysis showed that only seven genes had similar patterns with iTRAQ results (Fig. 4). The results were consistent with the association analysis results between proteome and transcriptome. These results may be due to various post translational modifications and other complex regulatory networks in tea plant response to cold stress.

Figure 4. 

Analysis of transcript levels of the selected proteins among NA, CA and DA stages by qRT-PCR. Those genes which showed similar patterns are marked using black boxes. All data are the mean ± SD (n = 3).

As a perennial, evergreen and originating from tropical regions, low temperature is widely accepted as the most critical factor limiting tea plant growth and geographical distribution. An understanding of the adaptive mechanisms under low temperature stress of tea plant is necessary to enhance its cold tolerance. Although some studies have focused on the cellular, physiological, metabolic and transcrioptomics changes in tea plant during CA procedures^{[4, 13, 14, 16, 28−30]}, but the molecular mechanism at the proteome level remains unclear. In the present study, we conducted comparative proteomic analysis and iTRAQ labeling method to examine the whole protein profile changes at the different CA states in tea plant leaves. A large number of differential abundance proteins, putatively related to cold stress, were identified. Comparisons of the differentially accumulated proteins revealed that more differentially accumulated proteins were detected among the comparisons of CA vs NA and DA vs NA. The result indicates that many changes which resulted from CA remained unchanged in DA. This may indicate that de-acclimated plants might have greater cold resistance than NA plants, and might be more amenable to re-acclimation than un-acclimated plants. A similar response was observed in Arabidopsis with regards to transcript changes of those genes involved in photosynthesis, calcium signaling and general stress responses maintaining acclimated expression patterns^[31]. However, more work needed to be done to test this hypothesis. GO enrichment analysis is a commonly used tool to determine the potential function of differentially accumulated proteins in different data set comparisons^[23]. In this study, the significant different GO terms among different comparisons between NA, CA and DA were obtained using GO enrichment analysis. As has been observed previously in other systems^{[21, 32−35]}, our results showed that the differentially accumulated protein species during CA and DA mainly involve in cell wall, photosynthesis, energy, protein synthesis, metabolism, antioxidation, carbohydrate metabolic process, and binding. Recent studies showed that these biological processes or cellular components were common to CA. Degand et al. reported that the differentially accumulated proteins identified from chicory root after CA were mainly classified into the functional categories of protein synthesis, metabolism, energy or cell structure^[36]. Kosmala et al. found that proteins related to photosynthetic machinery, cell energy and cell metabolic pathways played important roles in the CA procedure of Festuca pratensis^[32]. The proteomic analysis of Thellungiella rosette leaves under cold stress revealed that most identified proteins were involved in photosynthesis, defense response, cell wall and cytoskeleton, RNA metabolism, energy pathway, protein metabolism and signal transduction pathways^[37]. The proteomic studies of the CA mechanism in sunflower found that those cold-responsive proteins from three different cold tolerant lines were mostly involved in metabolism, protein synthesis, energy, and defense processes^[10]. Moreover, the proteomic results in plantain also indicated that the majority of differentially accumulated proteins were involved in oxidation-reduction, photosynthesis, and several primary metabolic processes^[21]. Our results also showed that 20 pathways were significantly enriched in tea plant during CA process at transcriptome level, and the metabolism was the largest category, which included 'carbohydrate metabolism pathway', 'energy metabolism pathway', 'xenobiotics biodegradation and metabolism' and 'lipid metabolism', etc^[4]. Therefore, it can be assumed that the CA of tea plant is characteristic of many previously identified integral metabolic changes and suggest that many of the previous regulatory mechanisms controlling CA and DA can be used to direct future research into improving cold tolerance of tea plant.

Interestingly, according to the significantly different GO terms listed in Table 1, cell wall is a significantly over-represented ontology that is specifically associated with DA. Cell wall is the first physical barrier and it plays a vital role in plant responses to abiotic stress^[38]. Previous research suggested that many changes happen in the cell wall when plants were placed under cold stress, such as the increase in weight of cell walls, cell wall composition changes, and expression of cell wall-related gene changes and those changes showed close relationships with plant cold resistance^[39−41]. Our results were consistant with previous studies and provided some novel findings for tea plant CA mechanism research.

Tea plant is an evergreen woody plant, and thus adjusting photosynthetic processes to deal with alterations in membrane fluidity and structure are important^{[37, 42]}. We observed the terms of plastid part, thylakoid light-harvesting complex, chloroplast thylakoid membrane, light-harvesting complex, plastid thylakoid membrane, NADP metabolic process and DANPH regeneration being over-represented among differentially accumulating proteins in DA vs CA, plastid envelope in DA vs NA, plastid stroma in both CA vs NA and DA vs CA, and photosynthetic electron transport chain and photosynthesis in both DA vs CA and DA vs NA in concurrence with the hypothesis that modifications in photosynthesis are required for CA and DA in tea plant. These observations are consistent with the fact that photosynthetic processes can lead to greater production of damaging oxidative radicles when thylakoid membranes undergo changes of state during chilling^[43]. Furthermore, reactive oxygen species scavenging is an important mechanism needed for coping with oxidative stress under cold stress in plants^[44]. Consequently, over-representation of the terms of oxidation reduction and oxidoreduction coenzyme metabolic process suggest involvement of antioxidative mechanism. The terms of electron transport chain was significantly over-represented in both DA vs CA and DA vs NA, while in DA vs CA was energy reserve metabolic process. Moreover, the alterations in protein synthesis were indicated by over-representation of the terms of ammonia ligase activity, acid-ammonia (or amide) ligase activity and ribosome biogenesis in both DA vs CA and DA vs NA comparisons. Energy production related proteins were significantly up-regulated in rice leaf blades and were down-regulated in poplar leaves under a chilling environment^{[45, 46]}. The terms of metabolic process and other terms involved in metabolism were over-represented among protein annotations in DA vs NA. The terms of carbohydrate metabolic process in CA vs NA and racemase and epimerase activity, acting on carbohydrates and derivatives in DA vs CA are involved in carbohydrate metabolism. The term of binding has significant difference in both CA vs NA and DA vs NA. Significant different terms of membrane part and fatty acid process in DA vs NA were also detected. Membrane modifications – particularly regarding alterations in fatty acid composition, have long been associated with CA processes^[47]. Likewise, it is well known that carbohydrate accumulation can be protective against cold stress^[48]. And these results were consistent with our previous study at the transcriptome level^[4].

CA is an important mechanism for perennial plants to obtain or enhance freezing tolerance^[1]. Extended freezing temperatures in winter pose a great challenge for the survival of evergreen perennials such as tea plant, and, along with dormancy formation, such plants develop physiological and molecular changes to successfully archieve overwintering^[49]. Because freezing of extra-cellular water pose significant challenges in regards to dehydration, cellular changes in water content, water status and osmotic potentials are essential events^[50]. Consequently, cells decrease water content during fall or early winter and accumulate osmoprotectants such as specific storage proteins, sugars, starch and alter membrane chemistries to better deal with dehydration^[51]. These adaptive mechanisms rely in part on gene induction and regulation, resulting in related protein enrichment in the fully acclimated stage. Therefore, the functional category analysis revealed that cell wall, photosynthesis, energy, protein synthesis, metabolism, antioxidation, carbohydrate metabolic process, and binding take crucial roles in tea plant CA.

Pathway enrichment analysis was conducted to determine the main signal transduction and biochemical metabolic pathways involved by those differentially accumulated proteins in a previous study. The results showed that microbial metabolism in diverse environments, metabolic pathways and biosynthesis of secondary metabolites were the top three enriched pathways. In concordance with GO analyses above, carbon fixation in photosynthetic organisms, ribosome, carbon fixation in photosynthetic organisms, starch and sucrose metabolism, protein processing in endoplasmic reticulum, plant-pathogen interaction, oxidative phosphorylation and photosynthesis also take a large proportion of the differentially accumulated proteins. The changes of general metabolism and photosynthesis are the main responses of tea plant to low temperature. Accumulation of secondary metabolic products is an important character for tea plant. Many of the differentially accumulating proteins are involved in secondary metabolism related pathways. The role of these secondary metabolites in CA and DA are difficult to envision. For example, we observed increases in gallic acid and its derivatives which are known to have anti-oxidant properties, and thus might have protective roles in freezing stress^[48]. However, increases in gallic acid production are not commonly observed during CA processes of other non-evergreen systems. It is well noted in the literature that protein content in plant cell can change within hours of low temperature exposure, and that protein metabolism play important roles in CA and freezing stress tolerance^[52]. The enrichment of ribosome and protein processing in endoplasmic reticulum in this study also indicates the close correlation between protein synthesis and CA. Interestingly, a large number of proteins were involved in the pathways of microbial metabolism in diverse environments and plant-pathogen interaction. This was consistent with the results found in the plasma membrane of oat and rye after CA^[53]. Recent studies report that some of the pathogenesis-related proteins are induced during winter months and have been shown to have antifreeze activity, cryoprotective activity, or antifungal activity^[49].

Furthermore, the signal transduction pathway plays a pivotal role in the response to the stress of low temperatures^[54]. Plasma membrane play an important role in response to low temperature stress^[47]. Plasma membrane can sense and transduce cold signals and then signal responses that alter its structure, chemical composition and function to improve cold resistance^{[11, 53, 55]}. Li et al.^[56] and Takahashi et al.^[53] had identified hundreds of differentially abundant plasma membrane proteins in Arabidopsis, oat and rye using shotgun proteomic technology suggesting that these proteins may have a role in CA and freezing tolerance development. These proteins included signal transduction, disease/defense-related, energy-related, transporter and post-translationally modified proteins, etc. Phosphatidic acid is one of the major membranous second-messenger molecules and is produced by phospholipase D. During CA, plants generally increase their phospholipase D levels and phosphatidic acid content which is correlated with enhanced low temperature resistance^[8,57]. Heat shock proteins are a kind of membrane protein that act as molecular chaperones^[58], and usually play crucial roles in response to cold stress by re-establishing normal protein conformation and thus cellular homeostasis^[53,59]. In tea plant, we also found phospholipase D (comp1149_c0_seq1) and heat shock proteins (comp671_c0_seq1, comp1412_c0_seq1 and comp16542_c0_seq2) were up-accumulated in CA and DA compared with NA. The plasma membrane contains many proteins, and some of those involved in cold-responsive process in tea plant were also reported in other plants^[53,56]. These include, but are not limited to, proteins such as early response to dehydration proteins (ERDs), fasciclin-like arabinogalactan proteins (FLAs), aquaporins, ATPases, clathrins weren't differentially accumulated in our results (Supplemental Table S1). Studies using isolated plasma membrane should be conducted to elucidate the accurate changes of plasma membrane-specific proteins during CA in the future.

In the present study, 1,331 proteins were identified from NA, CA and DA tea leaves using iTRAQ analysis. 407 and 477 proteins were differently accumulated in comparison NA vs CA and DA Vs CA respectively. Function and KEGG pathway analysis revealed that those differently accumulated proteins were mainly mapped onto the metabolic, biosynthesis of secondary metabolites, microbial metabolism in diverse environment, ribosome, sugar metabolism, protein processing, photosynthesis and plant-pathogen interaction pathways. Further GO enrichment analysis indicated that those proteins were mainly involved in protein synthesis, photosynthesis, energy, sugar metabolism, antioxidation and stress defense. Correlation analysis showed that the proteome changes were not well-correlated with corresponding gene transcription changes. Overall, our study revealed general information about the proteome changes in tea plant leaf during NA, CA and DA procedures and provided some new insights on cold tolerance mechanism in tea plant.