Overview of the single-molecule real-time (SMRT) sequencing
To obtain the gene expression profiles under sheltered and unsheltered conditions, SMRT sequencing and Illumina RNA-Seq were conducted for leaves (Figure 1a) and fruits (Figure 1b) on the 35th, 45th, and 55th days after flowering (DAF35, DAF45, and DAF55). In total 1,048,866 post-filter polymerase reads (21.55G) were scanned. The subreads from the same polymerase read sequence generated a circular consistent sequence (CCS), which yielded 685,339 CCS sequences. Among them, 542,795 full-length non-chimera (FLNC) sequences with 5'-primer, 3'-primer, and poly-A were obtained, including 79.20% of all CCSs being FLNCs (300–22,293 bp). This proportion varied slightly between the two tissues, at 34.96% in the leaf data set and 44.25% in the fruit data set. The mean length of the FLNC reads is presented in Table 1, and the length distribution of the FLNC reads is shown in Figure S1. The FLNC length of leaf and fruit libraries of >2 kb accounted for 63.95% and 65.21% of the corresponding FLNC, respectively. Comparison between the protein-coding gene transcripts and those of FLNC revealed a strong concordance, which exhibited a better recovery of large transcripts than previous Illumina RNA-Seq data for gene model prediction, particularly in the size range of 2,500–4,000-bp size range [23].
Annotation and functional classification of unigenes
The functions of unigenes were annotated by BLAST comparison and were predicted by comparative analysis with five databases [NCBI Nonredundant Protein Sequences (NR), Clusters of Orthologous Groups of Proteins (COG), a manually annotated and reviewed protein sequence database (Swiss-Prot), Kyoto Encyclopedia of Genes and Genomics (KEGG), Gene Ontology (GO)]. Among the 45,825 unigenes distributed to each of the databases, 45,747 (97.35% of the total) for Nr, 28,481 (60.61%) for COG, 38,629 (82.20%) for SwissProt, 18,202 (38.73%) for KEGG, and 26,065 (55.46%) for GO were investigated (Figure 2a).
Analysis of the Nr database indicated that the highest homologies of cherry were with Prunus mume or with P. persica, with 22,919 (50.09%) and 18,545 (40.53%) unigenes annotated, respectively. In total, 28,481 annotated unigenes were classified into 37 functional groups of the three GO main categories: nine groups for molecular function (MF), 18 for biological process (BP), and 10 for cellular component (CC) (Figure 2b, Supplementary Table S2). The top three GO terms for the classified genes were “protein binding” (4,812), “ATP binding” (3,907), and “protein kinase activity” (2,272) for MF; “protein phosphorylation” (2,273), “oxidation-reduction process” (1,833), and “signal transduction” (1,255) for BP; and “membrane” (1,656), “integral component of membrane” (1,571), and “nucleus” (748) for CC.
In total, 18,202 unigenes were mapped into 265 KEGG database pathways. The top three KEGG pathways were “Metabolism, Carbohydrate Metabolism” (1,719); “Environmental Information Processing, Folding, Sorting, and Degradation” (1,563); and “Genetic Information Processing, Translation” (1,491) (Figure 2c, Supplementary Table S3). The maps representing the highest number of unigenes included those of carbon metabolism (map01200), RNA transport (map03013), biosynthesis of amino acids (map01230), protein processing in endoplasmic reticulum (map04141), Spliceosome (map03040), and starch and sucrose metabolism (map00500).
Functional classification of DEGs under sheltered covering
After filtering the low-quality reads, 17.25 billion clean reads were acquired by Illumina RNA-Seq (Table S4). According to the Illumina data, pair-wise comparisons (SL vs. UL, SF vs. UF) of gene expression among the three stages were performed. In response to microclimatic change with low-PAR conditions, in total, 38,621 (SL35 vs. UL35, SL45 vs. UL45, and SL55 vs. UL55) and 3,584 (SF35 vs. UF35, SF45 vs. UF45, and SF55 vs. UF55) DEGs were detected from the GO database. Additionally, 38,621 DEGs for SL45 vs. UL45 and 2,871 DEGs for SF45 vs. UF45 were acquired (Figure 3a). Throughout the three developmental stages, a total of 6,911 and 1,755 DEGs were detected from the KEGG database in leaves and fruits; 6,868 and 1,552 DEGs for SL45 vs. UL45 and SF45 vs. UF45 were obtained, respectively (Figure 3b). To assess the adaptation of leaves and fruits to the sheltered covering, the overall changes in biochemical pathways were analyzed using MapMan, including the DEGs of SL vs. UL and SF vs. UF at DAF45 (Figure 3c, d). The overall changes in the biochemical pathway of these DEGs were shown by MapMan annotation. The metabolic activities of the leaves under the sheltered covering were greatly increased at DAF45, and the DEGs were assigned to different functional categories, such as photosynthesis, secondary metabolism, coenzyme metabolism, RNA biosynthesis, etc (Figure 3c). Under shelter covering for 45 days, the DEGs of fruit characterized by the up-regulation of secondary metabolism, cell wall organization, amino acid metabolism, etc (Figure 3d).
Compared with UL and UF, the up-regulated genes from GO terms of SL were significantly enriched in “catalytic activity,” “biological process,” “oxidoreductase activity,” and “oxidation-reduction process,” among others (Figure 4a). The up-regulated genes of SF were primarily associated with “catalytic activity,” “biological process,” “oxidoreductase activity,” and “oxidation-reduction process,” among others, at DAF45 (Figure 4b). Up-regulated genes from KEGG terms of SL in contrast to UL included those associated with “circadian rhythm-plant,” “glyoxylate and dicarboxylate metabolism,” “porphyrin and chlorophyll metabolism,” and “carbon fixation in photosynthetic organisms,” among others (Figure 4c); moreover, the up-regulated genes of SF were particularly associated with “biosynthesis of amino acids,” “phenylpropanoid biosynthesis,” “phenylalanine metabolism,” and “flavonoid biosynthesis,” among others, at DAF45 (Figure 4d).
Common expression patterns were employed to further analyze the DEGs between SL vs. UL and SF vs. UF at three stages; overall, 7,244 (leaf) and 1,707 (fruit) DEGs were placed into four clusters (Figure 5). Most of the candidate DEGs was categorized into either leaf Cluster 1 (3796 genes) or fruit Cluster 1 (1271 genes). By KEGG pathway enrichment, the biological pathway of enrichment in each similar regulatory transcription cluster was described. The DEGs in leaf Cluster 1 exhibited peak expression at DAF45 of SL, including a broad range of genes responsible for porphyrin and cholorophyll metabolism, circadian rhythm-plant, photosynthesis-antenna proteins, carotenoid biosynthesis, peroxisome, and carbon fixation in photosynthetic organisms. The DEGs were mainly involved in photosynthesis, which reflected that leaves can maintain an effective adaptability from DAF35 to DAF45 after a long duration of rain-sheltered covering. For the top five accumulated pathways in leaf Cluster 2, the DEGs showed peak expression at DAF45 of UL, including peroxisome, phenylpropanoid biosynthesis and glutathione metabolism, and these DEGs were mainly involved in the antioxidant system (Figure 5a). Meanwhile, for the top four enriched pathways in fruit Cluster 1, the DEGs showed peak expression at DAF35 of SF, including a broad range of genes responsible for phenylpropanoid biosynthesis, starch and sucrose metabolism, ascorbate and aldarate metabolism, and phenylalanine metabolism. Finally, for the top five enriched pathways in fruit Cluster 2, the DEGs showed peak expression at DAF55 of SF, including flavonoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism, and phenylpropanoid biosynthesis (Figure 5b). Transcriptional results indicated that the adaptability of leaves to the microclimate was primarily attributed to the regulation of photosynthetic characteristics, assimilation, antioxidant status, and circadian rhythm. The genes implicated in the biosynthesis of sugars and anthocyanins in the sheltered fruits were up-regulated at DAF35 until DAF55, reflecting the prior accumulation of nutrition compared with that under shelter-free conditions.
To confirm their authenticity, 12 DEGs were randomly selected to analyze their expression profiles by qRT-PCR. The results of qRT-PCR analysis showed that the expression profiles of the 12 DEGs were similar to those obtained by high-throughput sequencing (Supplementary Table S2). These results confirmed the reliability of the genome-wide transcriptome profiling analysis.
To verify the authenticity of the RNA-Seq results, 12 DEGs were randomly selected and their expression profiles were analyzed by qRT-PCR. The results of qRT-PCR analysis showed that the expression profiles of these DEGs were similar to those assessed by RNA-Seq (Figure 6), confirming the reliability and accuracy of our RNA-Seq data.
Expression of genes involved in photosynthetic system in sheltered leaves
To clarify the molecular adaptability of cherry trees to the microclimatic change during fruit development upon exposure to the shelter covering, genes involved in environmental sensitivity were screened out from the filtered DEGs for further investigation (Figure 7, 8). Most of the DEGs encoding antenna proteins, or proteins involved in electron transport, reaction center in photosystem I (PSI) and photosystem II (PSII), and components of CO2 fixation, were highly expressed in SL compared with UL at DAF45 (Figure 7a).
More than 70 genes were annotated to three metabolic pathways: photosynthesis-antenna proteins (ko00196), photosynthesis (ko00195), and carbon fixation in photosynthetic organisms (ko00710). Therefore, we focused on the transcriptional levels of those genes closely related to the photosynthetic efficiency (Figure 7a). All annotated DEGs of light-harvesting Chl a/b binding protein complex I and II (LHCs), namely, 17 genes encoding the Chl A/B binding protein complex I (Lcha1, Lcha2, Lcha3, Lcha4, and Lcha5) and 20 genes encoding Chl A/B binding protein complex II (Lchb1, Lchb2, Lchb3, Lchb4, Lchb5, Lchb6, and Lchb7), were up-regulated in SL at DAF45 (Figure 7a, 8a). Moreover, 13 genes encoding proteins involved in reaction center and electron transport in photosynthesis, including the PSI reaction center subunit X (PsaK), reaction center subunit VI (PsaH), PSI reaction center subunit PsaN (PsaN), PSII oxygen-evolving enhancer protein 2 (PsbP), PSII 10 kDa protein (PsbR), PSII repair protein Psb27-H1 (Psb27), ferredoxin of photosynthetic electron transport (PetF), and H+/Na+-transporting ATPase subunit beta (AtpF), were up-regulated under sheltered covering (Figure 7a, 8b), whereas one gene of the cytochrome f complex (PetA) and one gene of H+/Na+-transporting ATPase subunit alpha (ATPF1AI) had lower transcription levels. In the carbon fixation pathway, D-ribulose 1,5-bisphosphate (RuBP) and CO2 produced 3-phosphate-glycerate (3-PGA) under the action of ribulose-bisphosphate carboxylase large chain (rcbL) and ribulose-bisphosphate carboxylase small chain (rcbS); then, 3-PGA was reduced to glyceraldehyde-3P (3-PGAld) by glyceraldehyde 3-phosphate dehydrogenase (GADPH) and glyceraldehyde-3-phosphate dehydrogenase (GAPA), which completed the energy storage process of photosynthesis and increased the production and accumulation of photosynthate. Among them, one gene of rbcL, nine genes of rbcS, one gene of GADPH, and seven genes of GAPA were significantly up-regulated under sheltered covering. Notably, the expression of rbcL and GAPA was up-regulated ≥10-fold (Figure 7a, 8c).
In combination with the PAR-Pn and CO2-Pn curves, the Pn of sheltered leaves was lower, but there were no significant differences between the first two stages of fruit development (Figure 9). The AQY and ACE of sheltered leaves visibly increased by 13.0% and 23.5%, respectively; meanwhile, the LCP and CCP parameters decreased to 13.87 and 75.62 µmol·m−2·s−1 at DAF45 (Table S5). In general, combined with photosynthetic characteristics and transcriptomic results of leaves, the findings illustrated that sheltered leaves had stronger abilities to capture and utilize the weak light, while maintaining stable CO2 utilization efficiency. This indicated the possibility of good adaptation to the weak light conditions under the sheltered covering in a short time.
Expression of genes encoding photosynthetic pigments in sheltered leaves
The adaptability of plants to weak light is inextricably related to the photosynthetic pigmentation synthesis pathway [24]. In this study, 34 genes involved in porphyrin and chlorophyll metabolism (ko00860) and one gene related to each of terpenoid backbone biosynthesis (ko00900) and carotenoid biosynthesis (ko00906) in sheltered leaves were up-regulated at DAF45 (Figure 7b). Regarding the DEGs annotated to Chl synthesis and carotenoid biosynthesis, one gene encoding geranylgeranyl diphosphate reductase (CHLP), one gene encoding oxygen-independent coproporphyrinogen III oxidase (hemN), two genes encoding glutamyl-tRNA reductase (hemA), 23 genes encoding magnesium chelatase subunit H (ChlH), and two genes encoding magnesium chelatase subunit I (ChlI), among others, were up-regulated under sheltered covering at DAF45 (Figure 8d). Additionally, four genes of 15-cis-phytoene synthase (crtB), two genes of 15-cis-phytoene desaturase (PDS), one gene of prolycopene isomerase (crtISO), and two genes of lycopene beta-cyclase (lcyB) were up-regulated under sheltered covering at DAF45 (Figure 8e). The contents of photosynthetic pigments Chl a, Chl b and Car tended to increase under sheltered covering. Moreover, the contents were higher than those under Cont; for example, Chl a increased by 14% to 16.7%, Chl b by 13.6% to 24%, and Car consistently increased by 20% to 22.6% under sheltered covering (Figure 10). The transcription levels of genes related to Chl and carotenoid synthesis in the sheltered leaves also fully confirmed that they would not be adversely affected by the sheltered covering (Figure 7b, 8e); conversely, sheltered covering would enhance the gene transcription levels, thereby increasing the pigment content.
Expression of genes involved in antioxidant systems in sheltered leaves
ROS scavenging enzymes such as SOD, catalase (CAT), peroxidase (POD), and ascorbate peroxidase (APX), as well as nonenzymatic antioxidants (glutathione, carotenoids, etc.) are essential for ROS detoxification [25]. Malondialdehyde (MDA) can be used as an indicator of lipid peroxidation under different stress conditions. Transcription levels of multiple genes relating to antioxidant capacity including peroxisome (ko04146), phenylpropanoid biosynthesis (ko00940), glutathione metabolism (ko00053), and carotenoid biosynthesis (ko00906), fluctuated slightly with the microclimatic conditions (Figure 7c). The expression of most antioxidant-related genes was clearly higher in SL than in UL at DAF45; for example, eight genes involved in peroxisome (encoding SOD, and CAT), 32 genes involved in phenylpropanoid biosynthesis [encoding trans-cinnamate 4-monooxygenase, CYP; phenylalanine ammonia-lyase, PAL; 4-coumarate-CoA ligase, 4CL; cinnamyl-alcohol dehydrogenase, CAD; caffeic acid 3-O-methyltransferase, COMT; 5-O-(4-coumaroyl)-D-quinate-3′-monooxygenase, C3’H; caffeoyl-CoA O-methyltransferase, CCoAOMT; and POD]; eight genes involved in glutathione metabolism (encoding glutathione synthetase, GSS; glutathione reductase, GSR; and glutathione S-transferase, GST); and 16 genes involved in carotenoid biosynthesis (encoding zeaxanthin epoxidase, ZEP; violaxanthin de-epoxidase, VED; beta-carotene 3-hydroxylase, CrtZ; beta-ring hydroxylase, LUT5; lycopene beta-cyclase, and lcyB; and abscisic acid 8′-hydroxylase, CYP707A) showed markedly higher expression in SL than in UL at DAF45 (Figure 7c).
The antioxidant enzyme activities and MDA content in leaves also changed with the change of microclimate associated with sheltered covering. At DAF35, there was no remarkable difference in the activities of SOD, POD, and CAT of SL and UL; however, the MDA content was higher under sheltered covering. Moreover, the activities of antioxidant enzymes were higher than under shelter-free conditions at DAF45, whereas the MDA content was lower than in UL, with the POD, SOD, and CAT activities increasing by 1.7-, 1.9-, and 1.3-fold, respectively (Figure 11). These results showed that the activities of antioxidant enzymes in leaves increased during RSC, which maintained the ability of plants to scavenge ROS, and prevented the damage of membrane lipid peroxidation to leaves, thus preserving photosynthetic efficiency. The results showed the same trend as the results of gene expression levels.
Expression of genes associated with anthocyanin synthesis in sheltered fruits
The anthocyanin content of mature cherry fruits under sheltered covering increased significantly compared with that of mature cherry fruits under unsheltered covering [22]. The up-regulated genes of fruits were mainly involved in phenylpropanoid biosynthesis (ko00940) and flavonoid biosynthesis (ko00941). The transcriptional levels of more than 25 genes in sheltered fruits began to increase from DAF35 until fruit ripening (Figure 7d); these included 12 genes involved in phenylpropanoid biosynthesis, with one gene of PAL, six genes of C4H, and five genes of 4CL. Among these, the expressions of C4H and 4CL were up-regulated more than 8.5- and 6.1-fold, respectively. Moreover, for 13 genes involved in flavonoid biosynthesis encoding chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3'-hydroxylase (F3'H), dihydro-flavonol 4-reductase (DFR), and anthocyanidin synthase (ANS), the expressions were remarkably higher in the sheltered fruits at DAF45 (Figure 12). The results showed that the synthesis of anthocyanin could be accelerated and the accumulation cycle could be prolonged in the sheltered fruits.