Gene expression profiling during peach ripening and novel transcripts annotation
To reveal the molecular responses to HT and CT in post-harvest peach fruits, samples with two biological replicates from 6 time points were processed to construct 24 RNA libraries. A total of 684 M reads were generated with an average length of 100 bp, ranging from 10.18 to 21.57 M raw reads per library (Supplementary Table S2 ). The reads data were pre-processed by trimming adapters and filtering out short reads (< 36 bp) and low-quality reads (Q < 30). After pre-processing, 562 Mb (82.22%) high-quality sequences were retained. The remaining paired reads were mapped to the Prunus persica reference genome, and the total mapped rate was approximately 92.26% in average (Supplementary Table S3 ).
Relative abundance of genes was calculated and normalized as FPKM, and the robustness and reproducibility of data were indicated by high Pearson correlation coefficients (> 0.97) among biological replicates for all comparison pairs (Supplementary Table S4 ). Besides, we identified 1,281 novel transcripts which were aligned to the reference genome but not included in the 26,873 protein-coding genes (Supplementary Data ). The lengths of those transcripts ranged from 46 to 215,643 bp, with an average length of 3,045 bp. The biological functions of novel transcripts were annotated with Gene Ontology (GO) terms by successively using BLASTx (against the NCBI database) and Blast2GO; they were also subjected to KAAS to obtain the KEGG Orthology (KO) identifiers. Finally, 353 transcripts were classified into 3 GO categories and 176 transcripts have KO annotations (Supplementary Figure S1 & Table S5 ).
Analysis of DEGs
To identify DEGs between HT and CT conditions, the samples of each time point (day) were compared and the DEGs were detected. Firstly, we counted the expressed genes in each sample and found that peach fruits at CT underwent a more substantial decrease in the number of expressed genes compared with HT (Figure 2a ). A total of 7,245 DEGs (P< 0.05 and | Log2 ratio | ≥ 1) were identified in all comparison pairs between HT and CT. The numbers of up- and down-regulated DEGs were shown inFigure 2b . Interestingly, there is a strong decline in the numbers of both up-regulated and down-regulated genes at ripening stage. The number of up-regulated DEGs reached summit on day 2 with 2,778 genes, showing the most significant response to HT in transcriptome.
Several groups of transcription factors (TFs) were identified involved in plant responses to diverse abiotic stresses including temperature and drought. Here, 301 TFs from all DEGs were identified belonging to 49 TF families. Among all the TF families, ethylene response factors (ERFs) showed the highest numbers (44) compared with others (Figure 2c ). These data suggest the important role of ERFs in high-temperature response.
GO analysis was performed to dissect the functions of DEGs. As a result, 2,665 up-regulated DEGs and 1,990 down-regulated DEGs were classified into 3 main GO categories: cellular component (CC), biological process (BP) and molecular function (MF). In up-regulated DEGs, the most abundant BP subcategories were ‘nucleobase-containing compound metabolic process’, ‘gene expression’ and ‘nucleic acid metabolic process’ (Figure 3a ). In down-regulated DEGs, the major BP subcategories were ‘organic substance biosynthetic process’, ‘cellular macromolecule biosynthetic process’ and ‘cellular nitrogen compound biosynthetic process’ (Figure 3b ). Besides, only down-regulated DEGs had enriched groups in CC and MF which are ‘cytoplasm’, ‘structural constituent of ribosome’ (Figure 3b ), etc.
Mitogen–activated protein kinase (MAPK) cascades response to HT stress
HT induces changes of cytoplasmic calcium signaling, which exerts regulatory effects on plant development and stress responses (Berridge et al., 2003; Perochon et al., 2011). Among the 7 genes of calmodulin (CALM), Prupe.3G266000 and Prupe.4G269100 showed higher transcript abundance in CT fruits, while Prupe.1G089100 was on the contrary (Figure 4a ).
The MAPK cascades are one of the major signaling pathways and highly conserved in evolution for the transformation of extracellular stimuli into intracellular responses in eukaryotic cells (Hirt, 1997; Tena et al., 2001; Zhang and Klessig, 2001; Nakagami et al., 2005). The MAPK cascades consist of three interconnected protein kinase modules: MAPK kinase kinase (MEKK), MAPK kinase (MKK), MAPK. In MAPK gene families, we found that MEKK1-MKK2-MPK4 and MPK6 had high expression levels responding to HT (Figure 4b ), which were also detected in A. thalianaresponding to salt, drought and cold (Teige et al., 2004). Therefore, MEKK1-MKK2-MPK4/6 module is the most likely signaling pathway activated during HT stress.