Proteogenomic analysis of pitaya reveals cold stress-related molecular signature

Plant material and RNA extraction

The pitaya cultivar “Tiegusu” was used in this study. “Tiegusu” has light-green flowers and is one of most widely grown commercial cultivars in China. The pitaya plants were grown in a greenhouse at the Guizhou Academy of Agricultural Sciences (Guizhou, China) under a temperature range of 23–28 °C and natural light. To identify as many transcripts involved in cold stress as possible, RNA was extracted from 12 samples of six different tissues (young roots, tender shoots, stems, flower buds, new stems and fruits) from plants that had been grown at 0 °C (cold treatment) for 3 days or grown at a normal temperature (25 °C treatment). A TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate total RNA from each sample and RNase-free DNase I (TaKaRa) was used to treat the samples by following the manufacturer’s protocol. Equal amounts of total RNA from each treatment were pooled together for cDNA library construction and Illumina deep sequencing.

cDNA library preparation and Illumina sequencing for transcriptome analysis

RNA-Seq was performed at Shenting Genomics Institute (Hangzhou, China) using an Illumina HiSeq™ 2000 (Illumina Inc., San Diego, CA, USA). Briefly, poly (A)⁺ mRNA was isolated from the pooled total RNA sample using Oligo (dT) magnetic beads. The mRNA was fragmented into short fragments using a fragmentation buffer to avoid priming bias. A SuperScript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad, CA, USA) with a random hexamer-primer (Illumina) was used to synthesize double-stranded cDNA. The synthesized cDNA was then subjected to end-repair and phosphorylation, and the repaired cDNA fragments were 3′-adenylated with Klenow exo- (3′–5′ exo minus, Illumina). The ends of the 3′-adenylated cDNA fragments were ligated by Illumina paired-end adapters. The ligation products were purified on 2% agarose gel to select the appropriate templates for downstream enrichment. The cDNA fragments (approximately 200 bp) were excised from the gel. After end reparation and ligation of the adaptors, the products were amplified by PCR using PCR primers PE 1.0 and 2.0 (Illumina) with fusion DNA polymerase and purified using the QIAquick PCR Purification Kit (Qiagen, Valencia, CA, USA). Finally, the cDNA library was constructed using 200 bp insertion fragments. The library was validated using an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), then sequenced using Illumina HiSeq™ 2000 (Illumina Inc., San Diego, CA, USA). The raw data was deposit at www.iprox.org with accession project ID: IPX0001296002.

Data filtering, de novo assembly and annotation

The raw reads were cleaned by removing adaptor sequences, empty reads and low-quality reads using Trimmomatic (v0.32). The cleaned reads (a minimal sequence length of 200 bp with identity value above 95%) were assembled by the de novo software Trinity (v2.0.6) with default parameter. Then the TGICL and Phrap Clustering tools were used to obtain sequences that could not be extended at either end. The obtained sequences were defined as unique transcripts (or unigenes) (Grabherr et al., 2011; Pertea et al., 2003; Vogel et al., 2006). To obtain annotation of transcripts, all unigene sequences were searched against NCBI nr (non-redundant protein) database using the BLASTX algorithm (http://www.ncbi.nlm.nih.gov/) with cutoff of E-value 10⁻⁵. Functional annotation was performed using Gene Ontology (GO) and analyzed using the Blast2go (http://www.blast2go.org) software (Conesa et al., 2005). The COG and KEGG Orthology (KO) annotations were performed using the A. thaliana and Oryza sativa genome sequences as reference data in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2012; Wixon & Kell, 2000; Ogata et al., 1999). HMMER (http://hmmer.org) was used to obtain a domain-based annotation with Pfam (http://Pfam.sanger.ac.uk) database (Finn et al., 2016). All unigene sequences were further scanned by TransDecoder (v2.0.1) to identify Coding sequence (CDS). In totally, 44,337 CDS were predicted.

Protein extraction and profiling

Each replicate of pitaya samples was pulverized with mortar and pestle to a fine powder in liquid nitrogen, respectively. About 1 g of sample was used for protein extraction using a filter-assisted sample preparation method. Briefly, the sample was suspended in five mL protein extraction buffer (0.5 M Tris–HCl (pH 7.5), 0.7 M sucrose, 0.1 M KCl, 50 mM EDTA, and 40 mM dithiothreitol (DTT)) for 10 min at room temperature. After that, equal volume (five mL) of Tris-phenol was added. After 30 min of shaking, the suspension was centrifuged at 8,000×g and 4 °C for 5 min. After centrifuging, the upper phenolic phase was collected for further extraction by adding an equal volume of extraction buffer to the supernatant. Then, four volumes of 0.1 M ammonium acetate in methanol were added and kept the mixture overnight at −20 °C for protein precipitation, then centrifuging at 4 °C, 8,000×g for 10 min and discarded the supernatant. The pellet was washed thrice at 4 °C with ice-cold acetone. Finally, the pellet was dried for 2 h in a vacuum drier. After drying, the rehydration solution (100 μL; 8 M (w/v) urea, 0.1 M (w/v) Tris, and 10 mM DTT) was used to solubilize the pellet. The concentration of protein was determined by Bradford method.

In the centrifuge tube, the deposit was buffer exchanged with 8 M urea containing 10 mM DTT and 100 mM Tris–HCl. Then, the deposit was alkylated with 55 mM iodoacetamide. The urea concentration in the extract was then diluted to 1 M using the Tris–HCl buffer (pH 7.6). Protein samples were digested by trypsin (enzyme to substrate ratio = 1:20) in a thermomixer (1,000 rpm) overnight at 37 °C. Nano LC-1DTM plus system (Eksigent, Dublin, CA, USA) combined with AB Triple TOF 5600 MS (Foster City, CA, USA) were used to analysis digested peptides. Firstly, eight μL crude polypeptide was injected using a full sample loopin. Crude polypeptide was then desalted on a ChromXP trap column (NanoLC TRAP Column, 3 μm C18-CL, 120 A, 350 μm × 0.5 mm; CA, USA) and then eluted into a second analytical column (Nano LC C18 reversed-phase column (3C18-CL, 75 μm × 15 cm; CA, USA)) using a linear gradient formed by mobile phases A (5% ACN and 0.1% FA) and mobile phases B (95% ACN and 0.1% FA) for 120 min gradient at a flow rate of 300 nL/min. AB Triple TOF 5,600 system was operated in data-dependent acquisition mode to automatically switch between TOF–MS and product ion acquisition using Analyst Software (TF1.6). β-Galactosidase digestion was used to calibrate every two samples by 10 min for elution and 30 min for identification.

Proteome data processing

Raw MS files from AB Triple TOF 5600 were processed by MaxQuant version 1.5.2 (http://www.maxquant.org). MS/MS spectra were searched against the transcriptome derived database above mentioned. Precursor mass and fragment mass were identified with an initial mass tolerance of 6 ppm and 20 ppm, respectively. The search included variable modifications of methionine oxidation and N-terminal acetylation, and fixed modification of carbamidomethyl cysteine. Minimal peptide length was set to seven amino acids and a maximum of two mis-cleavages was allowed. MS runs were analyzed with the “match between runs” option. For matching, a retention time window of 20 s was selected. Proteins matching to the reverse database were filtered out. The false discovery rate was set to 0.01 for peptide and protein identifications. The raw data of proteome were all deposit at http://www.iprox.org with accession project ID: IPX0001296001.

Illumina paired-end sequencing and de novo assembly of reference pitaya transcripts

In this study, from one plate (8 lanes) in a single sequencing run, a total of 25.3 million reads were obtained, generating approximately 4.2 giga base pairs (Gbp) of raw data (Table 1). After the removal of adaptor sequences, low-quality reads (Q-value < 25) and ambiguous reads, 21.3 million high-quality clean reads (3.1 Gbp, 84.4% of the raw data) remained. The quality of the clean reads data was assessed based on the base-calling quality scores. The scores obtained using Illumina’s base-caller Bustard. Phred-like quality scores at the Q30 level were obtained for 97.23% of the clean reads data. All the high-quality reads were assembled into 110,330 isotigs (81,252 unigenes) with a maximum size of 200 bp (Fig. 1A). The isotigs were more abundant than those previously reported for a pitaya transcriptome dataset by Hua et al. (2016). The greater number of isotigs obtained in our study may be the result of trying to acquire the most comprehensive coverage possible by sequencing RNA extracted from six different types of pitaya plant tissues. The assembled isotigs have an average contig length of 934 bp (868 bp, unigene) and an N50 of 1,445 bp (1,373, unigene) (i.e., 50% of the assembled bases were incorporated into contigs of 1,445 bp or longer). Although a large proportion of the contigs (48.77%) were between 200 to 500 bp, we obtained 62,413 contigs (51.23%) that were more than 500 bp in length (Fig. 1A). Most of the highly expressed unigenes were approximately 2,500 bp in length (Fig. 1B).

Functional annotation and classification of the pitaya transcripts

In order to acquire complete functional information, sequence-based alignments of the unigenes were performed against different public databases, including the NCBI non-redundant protein database (NCBInr), Swiss-Prot/UniProt, KEGG pathway, GO and KOG cluster using the BLASTX algorithm with a significant E-value threshold of < 10⁻⁵. The distribution of the log10 (e-value) of the NCBInr BLAST and Swiss-Prot searches were normally, and the center of distribution were all at approximately 10⁻⁵⁰ (Fig. 2A), which indicated that the annotation had a high degree of similarity with known sequences. The Hidden Markov Model method was used to search both the InterPro and Pfam databases for Domain/family searches, and BLASTX was used to search alignments against the Clusters of Orthologous Groups (COGs) database at NCBI, with the E-value thresholds also set at ≤ 1 e⁻⁵ (Table S1). Out of 56,388 hits in the Nr database, 43,127 unigenes also had hits in the Swiss-Prot database (1,487 unigenes only had hits in Swiss-Prot) (Fig. 2B). The number of sequences that were annotated by searching the GO, KEGG and Pfam databases are shown in Fig. 2C. In total, 53,462 sequences were assigned to 23 of the EuKaryotic Orthologous Group (KOG) categories (Fig. 3). The KOG tool is a eukaryote-specific version of the COG tool for identifying orthologous and paralogous proteins (26). The KOG database also provides information about the classification of gene products, including their evolutionary relationships (Tatusov et al., 2003; Koonin et al., 2004). Based on the assumption that every protein evolved from an ancestor protein, pitaya isotigs were compared with known sequences in the KOG database to predict and classify their possible functions. Among the 23 KOG categories, “S: Function unknown” was the largest group (12,569 isotigs; 23.51% of all isotigs), followed by “R: General function prediction only” (6,608; 12.36%) and “O: Posttranslational modification, protein turnover, chaperones” (4,395; 8.2%). “N: Cell motility” was the smallest group (25; 0.04%), followed by “D: Cell cycle control, cell division, chromosome partitioning (65; 0.12%), “B: Chromatin structure and dynamics” (225; 0.42%) and “V: Defence mechanisms” (234; 0.43%) (Fig. 3). To retrieve function of pitaya, we mapped the annotated sequences to GO and canonical KEGG pathways (Kanehisa et al., 2008). The results showed that the most highly represented GO term was the “GO:0006950, response to stress,” with 451 members, and “path:ko04626, plant–pathogen interaction” was enriched at KEGG pathways, with 114 members (Figs. 4D and 5).

Proteomics characterization of pitaya under cold stress

In order to analyze the mechanism of pitaya response to cold stress, proteomics approach was used to identify and determine the abundance of proteins in all of tissues. The proteins of those samples (cold treatment vs. control) were extracted and each group treatment conducted three biological replicates. iBAQ values derived from MaxQuant software were used to represent the abundance of identified protein. Hierarchy clustering analysis, a quality-control measure based on expression profiles among replicates, indicates that the results are highly reproducible (Fig. 6). A total of 1,712 non-redundant proteins were identified (including those samples with 24 h cold treatment and control; Table S2).

Comparison of the abundance profiles of proteins between control and cold treatment of pitaya

In this study, the label-free quantitative proteomic analysis characterized the differences in protein synthesis between cold treatment and control. We used volcano plot to show differential expressed proteins. Compared with the control group, the up-regulated and down-regulated proteins (fold change ≥ 2, p ≤ 0.05) in cold treatment group were 18 and 98, respectively (Fig. 7; Table S3). Total 116 DAPs were subjected into Blast2GO and sorted by major enrichment of biological processes (Fig. 8), including proton transport, ATP hydrolysis coupled proton transport, glycolytic process, carbon fixation, ATP synthesis coupled proton transport, ATP metabolic process, proteolysis, metabolic process, ion transport, and intracellular protein transport. All DAPs were successfully annotated with KEGG pathways, sorted by enrichment score (Fig. 8), including AMPK signaling pathway (path:ko04152), fructose and mannose metabolism (path:ko00051), pentose phosphate pathway (path:ko00030), glycolysis/gluconeogenesis (path:ko00010), starch and sucrose metabolism (path:ko00500), carbon fixation in photosynthetic organisms (path:ko00710), ascorbate and aldarate metabolism (path:ko00053), pentose and glucuronate interconversions (path:ko00040), amino sugar and nucleotide sugar metabolism (path:ko00520), carbon metabolism (path:ko01200). It was found that proteins under ATP related GO catalogues were decreased, which indicated that the cold stress reduced the activity of energy generate, storage or release, which was one of the key aspects of life activity.

Overview and analysis of differentially abundant proteins related to cold stress

Compared with the control group, there are 18 proteins were increased expression abundance after 24 h cold treatment. The expression level of VPS15 Serine threonine-kinase, TRPA2 alpha chloroplastic tryptophan synthase and PPDK1 chloroplastic phosphate dikinase were significantly higher than other proteins. However, majority of the DAPs were down regulated and the 20 most decreased DAPs were listed in the Table S3. Among them, it is interesting that a number of them are chloroplastic or mitochondrial related proteins. Of which, 6 DAPs, including Ketol-acid reductoisomerase, Glutamate decarboxylase, Malate dehydrogenase, Pentatricopeptide repeat-containing protein, Phosphoglycolate phosphatase and Monodehydroascorbate reductase, are related to chloroplastic proteins, and 5 DAPs, including Aldehyde dehydrogenase family 2, Monodehydroascorbate reductase, Glycine dehydrogenase, and Probable mitochondrial-processing peptidase are related to mitochondrial proteins. Some studies provide evidence that chloroplastic and mitochondrial proteins play an important role in response to cold stress.