Transcriptome Analysis of Low- and High-Sucrose Pear Cultivars Identifies Key Regulators of Sucrose Biosynthesis in Fruits

Abstract

Sucrose accumulation is one of the important factors that determine fruit enlargement and quality. Evaluation of the sugar profile of 105 pear cultivars revealed low-sucrose and high-sucrose (HS) types of pear fruits. To better understand the molecular mechanisms governing the sucrose content of pear fruits, this study performed transcriptome analysis during fruit development using low-sucrose ‘Korla’ fragrant pear and HS ‘Hosui’ pear, and a coexpression module uniquely associated with the control of high-sucrose accumulation was identified by weighted gene coexpression network analysis. These results suggested that there are seven candidate genes encoding key enzymes (fructokinase, glucose-6-phosphate isomerase, sucrose phosphate synthase and sucrose synthase) involved in sucrose biosynthesis and several transcription factors (TFs) whose expression patterns correlate with those of genes associated with sucrose biosynthesis. This correlation was confirmed by linear regression analysis between predicted gene expression and sucrose content in different pear cultivars during fruit development. This study provides insight into the molecular mechanism underlying differences in sucrose content across pear cultivars and presents candidate structural genes and TFs that could play important roles in regulating carbohydrate partitioning and sucrose accumulation.

Introduction

Pear (Pyrus spp.) is one of the vital Rosaceae trees cultivated worldwide, with an annual harvest of approximately 24 million tonnes (Reuscher et al. 2016). Poor-flavor quality of pears has often been recognized as a major limiting factor for acceptance by consumers. Sugar content is one of the critical flavor quality traits perceived by consumers, and the development of novel cultivars with sugar-enhanced content is a primary objective of breeding programs (Kitch et al. 1998, Roth et al. 2007, Cirilli et al. 2016). The level of sugar is mainly regulated by carbohydrate supply and metabolic transformation in sink tissues, such as fruit (Teo et al. 2006). As suggested in other Rosaceae families, sorbitol and sucrose are the major products of photosynthesis distributed throughout pear trees, with sorbitol accounting for approximately 70% of the transported carbohydrates (Zhang et al. 2014a). The transported sugars are stored or metabolized in cells, where they can affect fruit quality.

The comparison of individual sugar among pear cultivars suggests that sucrose content considerably varies among different varieties (Moriguchi et al. 1992, Zhang et al. 2014b), with Chinese cultivars exhibiting low-sucrose content. By contrast, Japanese pears (Pyrus pyrifolia Nakai) have been characterized as sucrose-accumulating fruit species because of their high-sucrose content (Moriguchi et al. 1992). Understanding the relationship between sucrose accumulation and fruit quality is important not only for pears but potentially for other sorbitol-transporting fruit species that synthesize and accumulate sucrose (Cirilli et al. 2016, Vimolmangkang et al. 2016). To date, a number of physiological and molecular mechanism studies have been conducted to investigate these differences, and limited information is available on the enzymatic control of sugar content in fruits (Moriguchi et al. 1992, Suzue et al. 2006, Yamada et al. 2007, Zhang et al. 2014b, Shen et al. 2017). The findings reported by Reuscher et al. (2016) provided novel insights into sugar accumulation as well as candidate genes for key reactions and transport steps in the development of pear fruit (Pyrus communis) based on quantitative proteomics. It is reported that histone acetylation at the promoter for the transcription factor (TF) PuWRKY31 is associated with sucrose accumulation in the Ussurian pear (Pyrus ussuriensis) fruit by activating the expression of a sucrose transport gene (Li et al. 2020). However, sugar accumulation in fruits is a complex process controlled by a regulatory network of multiple genes, affected by environmental conditions, and dependent on many interconnected physiological and metabolic processes (Cirilli et al. 2016). Although sugar accumulation has been extensively studied, it remains unclear which molecular mechanisms are responsible for phenotypic and genetic variability in sugar content of pear cultivars and what plays a key role in carbon metabolism and quality in fruits.

The accessibility of transcriptome and draft genome sequences of pears (Wu et al. 2013) along with next-generation RNA sequencing (RNA-seq) data provides an opportunity to reveal the molecular mechanism underlying sugar accumulation in pear fruits. As of now, comparative dynamic transcriptome analysis of genotypes/cultivars with different sugar accumulation patterns has not yet been performed in pears. Here, we tried to dissect the molecular mechanism of sugar accumulation in pear fruit using the RNA-seq-based transcriptome analysis and weighted gene coexpression network analysis (WGCNA) mainly by comparing two cultivars (low-sucrose ‘Korla’ fragrant pear and high-sucrose ‘Hosui’ pear) in three fruit-ripening stages. This work provides new insights into assimilating supply in fruit carbon metabolism and its effects on pear quality, opening up the possibility of altering metabolic flux toward sugars by the manipulation of specific genes in biosynthetic pathways.

Results

Analysis of sugars in various pear cultivars and preparation of transcriptome samples

Pear fruits accumulate different types of soluble sugars and sugar alcohols, mainly fructose, glucose, sucrose and sorbitol (Fig. 1A). A variable ratio indicates that sucrose content has higher variability among pear cultivars, ranging from 0.98% to 62.02% of the total sugar content, with the lowest content being 1.14 mg·g⁻¹ FW and the highest content being 72.92 mg·g⁻¹ FW (Fig. 1B;Supplementary Table S1).

Fig. 1

Samples and analysis strategy for RNA-seq. (A) Range and distribution of soluble sugar in mature pear fruit. These data are summarized alteration in four kinds of sugar contents among 105 cultivars. The horizontal lines in the interior of the box are the mean values. The box indicates the distribution for 50% of the data. Approximately 99% of the data fall inside the whiskers. The data outside these whiskers are indicated by circles. (B) The relative content of four soluble sugars in total sugar from 105 pear cultivars. (C) Samples for RNA-seq. The panel shows the changes in fruit fresh weight (mean of 10–20 fruits) and sucrose content, the phenotype of fruits at three development stages and three biological replicates. Data are the mean of three biological replicates ± standard deviation.

To understand the molecular basis of sucrose polymorphism in pear fruits, ‘Korla’ fragrant pear with low-sucrose (LS) content (8.90 ± 0.18 mg·g⁻¹ FW) and ‘Hosui’ pear with high-sucrose (HS) content (51.89 ± 4.39 mg·g⁻¹ FW) were selected for comparative dynamic transcriptome analysis. As shown in Fig. 1C, the two cultivars showed similar sigmoidal growth curves in terms of fruit weight but differed significantly in terms of changes in sucrose content with fruit development. The highest difference was observed at the mature stage when the sucrose content of ‘Hosui’ pear was at least five times higher than that of ‘Korla’ fragrant pear, which might also explain why the former is considered sweeter than the latter in sensory evaluation. Next, LS and HS samples were collected at three critical stages of fruit development for RNA-seq (Fig. 1C). The samples were designated LS1 [30 days after bloom (DAB)] and HS1 (30 DAB) representing the early stage, LS2 (90 DAB) and HS2 (90 DAB) representing the middle stage and LS3 (150 DAB) and HS3 (142 DAB) representing the mature stage. At each stage, the fruit phenotype was determined and three independent biological replicates were used.

RNA-seq of developing fruits of two pear cultivars with contrasting sucrose content

After RNA-seq, the quality of data collected from 18 cDNA libraries was assessed (Table 1). The number of raw reads for each library ranged from 43 to 64 million, with an average Q20 (sequencing base calls with an error rate of <1%) of >95% and a mapping rate of 72.55–77.41%. The Pearson correlation coefficient for the biological replicates of different samples varied from 0.981 to 0.99 (Supplementary Fig. S1). These data indicated that throughput and sequencing quality were high enough to warrant further analysis.

Table 1

Quality of reads for 18 RNA-seq samples

Sample name	Raw reads	Clean reads	Total mapped	Uniquely mapped	Q20 (%)	GC content (%)
HS1_1	43988012	43028848	31,662,696 (73.58%)	28,832,770 (67.01%)	95.56	46.55
HS1_2	46883708	45776008	33,574,272 (73.34%)	30,576,784 (66.8%)	95.37	46.74
HS1_3	43707514	42623776	31,276,002 (73.38%)	28,480,859 (66.82%)	95.32	46.67
HS2_1	45725016	44890508	33,553,184 (74.74%)	30,346,765 (67.6%)	98.21	46.93
HS2_2	57227574	56111096	42,658,785 (76.03%)	38,590,497 (68.78%)	98.07	46.99
HS2_3	45229524	44325688	33,636,390 (75.88%)	30,459,157 (68.72%)	98.07	46.87
HS3_1	56301764	55806690	42,249,278 (75.71%)	38,275,138 (68.59%)	97.71	47.46
HS3_2	50364886	49920208	38,097,615 (76.32%)	34,527,661 (69.17%)	97.73	47.40
HS3_3	61162956	60558572	45,687,616 (75.44%)	41,415,431 (68.39%)	97.62	47.42
LS1_1	57868850	57542236	42,442,531 (73.76%)	38,496,721 (66.9%)	97.70	46.78
LS1_2	59753810	59318476	43,050,761 (72.58%)	39,035,210 (65.81%)	97.33	46.78
LS1_3	59462246	59142288	43,413,777 (73.41%)	39,455,560 (66.71%)	97.33	46.69
LS2_1	60696720	60326718	44,730,164 (74.15%)	40,657,944 (67.4%)	96.90	46.29
LS2_2	60373854	60083126	43,529,620 (72.45%)	39,519,784 (65.78%)	96.99	46.04
LS2_3	55384878	55126430	39,339,305 (71.36%)	35,707,363 (64.77%)	97.40	46.28
LS3_1	64811566	64497148	45,768,179 (70.96%)	41,709,391 (64.67%)	97.10	46.74
LS3_2	58818804	58437120	42,942,025 (73.48%)	39,186,478 (67.06%)	97.77	46.94
LS3_3	54282554	53958976	38,640,760 (71.61%)	35,251,510 (65.33%)	97.43	46.82

‘Raw reads’ means the number of paired-end reads, ‘spliced reads’ means one read spliced-mapped to two exon reads and the rate (%) means the percentage of the spliced reads/clean reads.

Q20, percentage of bases with a Phred value of >20; HS1, high-sucrose sample at 30 DAB; HS2, high-sucrose sample at 90 DAB; HS3, high-sucrose sample at 142 DAB; LS1, low-sucrose sample at 30 DAB; LS2, low-sucrose sample at 90 DAB; LS3, low-sucrose sample at 150 DAB.

Density distribution profiles of fragments per kilobase of transcript per million mapped reads (FPKM) were constructed to reflect the gene expression pattern of each sample (Fig. 2A). Quartile–quartile (Q–Q) plotting of FPKM values indicated that log₁₀-transformed FPKM values approximately followed a normal distribution; therefore, we used log₁₀-transformed data for further analysis. To investigate the global difference in genes of LS and HS samples during fruit development, we performed principal component analysis (PCA) based on Pearson correlation coefficient analysis of FPKM values (Fig. 2B). The result revealed two distinct groups within these samples, although HS1 functioned differently, suggesting that significant changes in transcriptional programs occur across the cultivars. The sample similarity based on PCA was confirmed by hierarchical clustering (Supplementary Fig. S2). In further accordance with the principle of correlation between genes and phenotypes, we believe that differentially expressed genes (DEGs) related to sucrose accumulation mainly exist between HS3 and LS3. Venn analysis was performed to identify genes showing a difference (Fig. 2C). Comparative analysis revealed that 792 and 659 DEGs were uniquely upregulated and downregulated in the HS3 vs. LS3 comparison group, respectively. In addition, three comparison groups (HS1 vs. LS1, HS2 vs. LS2 and HS3 vs. LS3) shared 2,508 DEGs including 1,066 upregulated DEGs and 1,272 downregulated DEGs.

Fig. 2

Comparison of transcriptomes of different stages of fruit development in the two pear cultivars. (A) Boxplot showing the distribution of FPKM values after the imputation of missing values, fraction assignment, log₁₀ transformation and linear regression normalization. Samples are shown that average three replications for each developmental stage in HS and LS samples. The thick black line represents median, box represents 25th–75th percentile, full line represents 1.5 times the interquartile range and dark spots represent outliers. (B) Principal component analysis plot showing clustering of transcriptomes at different development stages of HS (HS1–HS3) and LS (LS1–LS3) samples. Shown are the positions of each sample along with the first two principal components including their contribution to the total variation in percentage. Colors indicate the developmental stage of each sample. (C) Number of DEGs identified by pairwise comparison between HS and LS samples at each time point (adjusted P < 0.001 and absolute log₂ fold change ≥2). Up, upregulated; Down, downregulated.

We randomly selected 12 DEGs to verify the RNA-seq expression data from three stages of fruit development for both cultivars by real-time quantitative PCR (qRT-PCR). Although the exact fold change of DEGs at individual data points varied between RNA-seq and qRT-PCR, the expression profiles of the tested genes were largely consistent in the two approaches (Supplementary Fig. S3), indicating the accuracy of RNA-seq data in reflecting the abundance of transcript levels.

Comparison of transcriptional profiles of genes involved in sugar metabolism

Previous reports have described the complex biological processes involved in sugar metabolism including key enzymes involved in biosynthesis, degradation and transport (Supplementary Fig. S4). Therefore, candidate genes were searched from the transcriptome data based on standard gene names and synonyms in the combined functional annotations of the pathway. In total, 68 DEGs were associated with the sucrose metabolism pathway, which displayed a different expression pattern (Fig. 3). Not surprisingly, a higher transcriptional activity of six known genes associated with sucrose biosynthesis was found in HS3 but not in the other samples (LS1, LS2, LS3, HS1 and HS2), which was consistent with the tendency of sucrose accumulation in the studied cultivars (Figs. 1C, 3). These include two sucrose phosphate synthase (SPS)-encoding genes (Pbr009578.1 and Pbr042506.1), two sucrose synthase (SUS)-encoding genes (Pbr003394.1 and Pbr035996.1), a fructokinase (FRK)-encoding gene (Pbr012599.1) and a glucose-6-phosphate isomerase (PGI)-encoding gene (Pbr006178.1). However, other genes (Pbr032770.1, Pbr032771.1 and Pbr024722.1) encoding SDH and S6PDH enzymes, which are known to be involved in sorbitol metabolism, displayed opposite expression patterns between HS and LS samples (Fig. 3). This suggests that these genes play an important role in regulating sucrose accumulation.

Fig. 3

Expression profiles of sugar metabolism-related genes in different pear samples. Hierarchical clustering was performed to group DEGs with similar expression levels and patterns based on the normalized log₂-transformed FPKM values of DEGs. The blue line represents the DEGs that are consistent with the tendency of sucrose accumulation in two cultivars. CINV, cytoplasmic invertase or neutral invertase; HXK, hexokinase; PFK, phosphofructokinase; PLT, polyol transporter; SDH, sorbitol dehydrogenase; S6PDH, sorbitol-6-phosphate dehydrogenase; STP, sugar transporter protein.

In addition, we detected four sugar transporter subfamilies in all DEGs: polyol/monosaccharide transporter subfamily, sugar transporter protein subfamily, sugars will eventually be exported transporters (SWEET) and sucrose transporter subfamily. It was found that most transporters were highly expressed at the early stage, but the expression level of Pbr004497.1 encoding a putative SWEET protein was significantly higher in HS3 (Fig. 3) and might be essential as an efflux for sucrose accumulation.

Identification of gene coexpression modules related to high-sucrose accumulation

We first analyzed network topology of various soft-thresholding powers to have relatively balanced scale independence and mean connectivity of WGCNA. As shown in Fig. 4A, power 13 (the lowest power for which the scale-free topology fit index reaches 0.90) was selected to construct a hierarchical cluster tree (dendrogram) of 11,868 genes. This analysis identified 13 distinct modules (labeled with different colors), as shown in the dendrogram in Fig. 4B, in which major tree branches define the modules. The largest module contained 4,250 genes, whereas the smallest module contained only five genes. According to a correlation study using a network heatmap (Supplementary Fig. S5), all modules had almost no correlation with each other but had high correlation within the respective module. Furthermore, we associated each coexpression module with stages of fruit development via Pearson correlation coefficient analysis (Fig. 4C), in which nine coexpression modules showed a relatively higher correlation (r ≥ 0.60) with all stages. However, only the ‘green’ module was highly correlated with sucrose content (r = 0.95, P = 1 × 10⁻⁹) in the six samples (Fig. 4C).

Fig. 4

WGCNA of DEGs identified from HS and LS samples. (A) Network topology of different soft-thresholding powers. The left panel displays the influence of soft-thresholding power (x-axis) on the scale-free fit index (y-axis). The right panel shows the influence of soft-thresholding power (x-axis) on mean connectivity (degree, y-axis). (B) Hierarchical cluster tree showing coexpression modules identified by WGCNA. Dissimilarity was based on the topological overlap. Each leaf in the tree represents one gene. Major tree branches constitute 13 modules labeled by different colors. The y-axis is the distance determined by the extent of topological overlap. (C) Module–tissue association. WGCNA grouped genes into modules based on the patterns of their coexpression. Each of the modules was labeled with a unique color as an identifier, and each row corresponds to a module. Module–trait correlations and corresponding P-values (in parentheses). The left panel shows the 13 modules, and the right panel shows the number of member genes. The color scale on the middle shows module–trait correlations from −1 (blue) to 1 (red). At the bottom is the sample name. Each row corresponds to a module. The number of genes in each module is indicated on the left. Each column corresponds to a specific tissue. The color of each cell at the row–column intersection indicates the correlation coefficient between the module and tissue type.

Accordingly, we focused on the green module, which the gene heatmap showed to be overrepresented in HS3 (Fig. 5A;Supplementary Table S3). The result was consistent with the sucrose accumulation phenotype, demonstrating the accuracy of the screening module. Furthermore, gene ontology (GO) and pathway enrichment analyses were conducted (Supplementary Tables S4, S5) and the top 20 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways are presented in Fig. 5B, sorted by the number of genes. Remarkably, starch and sucrose metabolism (ko00500) were highly enriched in the KO pathway for 16 genes and six candidate genes (Pbr003394.1, Pbr009578.1, Pbr012599.1, Pbr035996.1, Pbr042506.1 and Pbr006178.1) from Fig. 3 were also present in this gene set of sucrose metabolism pathways (Fig. 5B), showing that these genes are most likely responsible for high-sucrose accumulation. Meanwhile, 42 TF-encoding genes from 22 families were identified in the green module (Fig. 5C). We constructed the transcriptional network linking TFs with their coexpressed genes related to sucrose biosynthesis. An active coexpression network of 24 TFs with four structural genes (FRK_Pbr012599.1, PGI_Pbr006178.1, SPS_Pbr027868.1 and SUS_Pbr035996.1) was found (Fig. 5D). Next, we analyzed the potential binding motifs of bZIP, C2H2, GRAS, MADS and WRKY for promoter sequences of their coexpressed genes (TF-binding motif-target genes) and the motifs of the TFs were enriched (Supplementary Table S6).

Fig. 5

Expression profile and transcriptional regulatory network of the green module associated with high-sucrose accumulation. (A) Eigengene bar plot of the green module. Heatmaps show the expression profile of all coexpressed genes in the module. The color scale represents Z-score. Bar graphs (below heatmaps) show the consensus expression pattern of coexpressed genes in the green module. The y-axis indicates the value of the module eigengene, and the x-axis indicates the sample type. (B) KEGG enrichment of DEGs in the green module. The top 20 enriched pathways are sorted by the number of genes. Starch and sucrose metabolisms were highly enriched with 16 genes, and seven candidate genes related to sucrose biosynthesis are marked with red. (C) Distribution of the number of TFs from the green module. (D) A consolidated transcriptional regulatory network of TFs with DEGs related to sucrose biosynthesis displayed by Cytoscape software. Cytoscape representation of coexpressed genes with edge weights of ≥0.60 in the green module. The rhombus stands for the sucrose metabolism gene, and the circle stands for TFs. The size of the node in the interactive network is proportional to the degree of the node, i.e. the more edges are connected to the node, the larger the degree and the larger the node. The color gradient of the node from blue to red corresponds to the value of the clustering coefficient from low to high. The clustering coefficient represents the connectivity between adjacent points of the node. Interaction between structural genes and TFs is indicated by dotted lines.

Evidence of the close connection of predicted genes and sucrose content

To verify the close connection between genes and sucrose content, sugar profiles of different pear cultivars were determined (Fig. 6A). Meanwhile, we obtained the transcript profiles of seven predicted genes at different fruit development and ripening stages of these cultivars (Fig. 6B). Linear regression analysis between gene expression and sucrose content, shown in Fig. 6C, revealed that almost all gene transcript levels were correlated with sucrose content across cultivars. The transcript levels of the genes encoding FRK and SPS have a strong positive correlation with the content of sucrose in all studied cultivars, while SUS2 and three TFs (bZIP, C2H2 and ERF) showed a negative correlation with sucrose content in individual cultivars. This analysis of different cultivars provides correlative evidence that they are directly or indirectly involved in the regulation of sucrose biosynthesis (Fig. 6).

Fig. 6

Connection of predicted genes and sucrose content. (A) Content of soluble sugar analyzed and include the sorbitol as well during fruit development in different pear cultivars. (B) Relative expression level of seven selected genes during fruit development in different pear cultivars. (C) Linear regression coefficient between gene expression and sucrose content. Results are presented as means of three biological replicates ± standard deviation, and significant differences are marked with *P < 0.05 or **P < 0.01.

Discussion

Evaluation of the sugar profile (relative content of each individual sugar) of 105 pear cultivars at maturity revealed the LS and HS types of pear fruits (Fig. 1). We propose two reasons for this phenomenon.

Functional structural genes related to sucrose content

SPS is a pivotal enzyme that catalyzes sucrose formation, which converts UDP: uridine diphosphate glucose and fructose-6-phosphate into sucrose-6-phosphate and further converts sucrose-6-phosphate into sucrose by sucrose phosphate phosphatase (Wind et al. 2010). Previous reports have shown that SPS activity is highly correlated with the sucrose content of fruits, such as melon (Hubbard et al. 1989, Burger and Schaffer 2007), watermelon (Yativ et al. 2010, Fan et al. 2014), tomato (Dali et al. 1992, Islam et al. 1996) and peach (Moriguchi et al. 1990). SPS gene overexpression in Arabidopsis leads to elevated sucrose levels in sink tissues and increased total dry weight (Park et al. 2008), and overexpressed maize SPS gene in tobacco leads to increases in the sucrose–starch ratio and overall shoot biomass (Baxter et al. 2003, Seger et al. 2015). The effect of SUS on sucrose content has been reported in Japanese pear fruit (Tanase and Yamaki 2000), in which there are two isoforms of SUS (SS I and SS II). SS I plays a role in the degradation of sucrose, whereas SS II plays a role in sucrose biosynthesis in mature fruits because SS II has a higher affinity for UDP glucose than for UDP. An increase in SUS and SPS was found to influence sucrose accumulation in Japanese pear fruits, but SUS activity contributed more than SPS activity toward this (Moriguchi et al. 1992). A study found that FaSPS3 (gene31122) and FaSUS1 (gene12940) contribute to ripening-related sucrose accumulation (Vallarino et al. 2015).

In the present study, the expressions of SPS genes were particularly high in the mature fruit of HS samples and presented a noticeable positive correlation with sucrose content (Figs. 3, 5, 6) and the expression levels of two SUS-homologous genes (Pbr003394.1 and Pbr035996.1) were found to be upregulated during ripening of HS samples that coincided with relatively higher sucrose accumulation (Figs. 3, 5). The results were similar to those of our previous study showing that the activities of SUS and SPS were highly correlated with sucrose levels in pears (Zhang et al. 2014b).

Another key candidate gene was Pbr012599.1 (Figs. 3, 5, 6), which is the coding enzyme of FRK. FRK is an enzyme that mainly phosphorylates fructose and plays a major role in sucrose metabolism and the distribution in sink tissues (Granot 2007). In the pericarp of tomato, FRK can effectively control the carbon flux distribution of sucrose to starch synthesis (Petreikov et al. 2001). Therefore, we speculated that high Pbr012599.1 (FRK) expression induced higher sucrose biosynthesis by converting fructose into fructose-6-phosphate.

PGI directly or indirectly participates in sucrose biosynthesis and catalyzes reversible reactions (Ferreira and Sonnewald 2012). We detected one PGI-encoding DEG (Pbr006178.1) that exhibited higher expression in HS3, which was consistent with the tendency of sucrose accumulation in HS and LS samples (Figs. 3, 5).

We also found Pbr004497.1, a gene that belongs to clade III of the AtSWEET family (Li et al. 2017), with a noticeable positive correlation with the sucrose content in HS samples but downregulated during fruit ripening in LS samples (Fig. 3). SWEET proteins have been widely identified as sugar transporters in plants, especially for sucrose transport (Fenske et al. 2000, Chen et al. 2012). Coincidentally, PuSWEET15 was observed to transport sucrose in pear fruit, overexpression of PuSWEET15 increased sucrose content, while silencing of PuSWEET15 decreased sucrose content (Li et al. 2020). This suggests that this transporter regulates sucrose movement between compartments and exerts strict control of their fluxes.

Combining this information with data from GO and KEGG enrichment analyses of the green module highlighting the role of biological processes/pathways in sucrose accumulation (Figs. 3–5), it can be inferred that FRK (one transcript), SUS (two transcripts), SPS (three transcripts) and PGI (one transcript) are most likely responsible for the high-sucrose biosynthesis in sucrose-accumulating pears. However, we did not find that the invertase, (CWINV: cell wall-bound invertase, cytoplasmicinvertase/NINV and VINV: vacuolar invertase) is a major participant in sucrose accumulation in pear fruit. As suggested in other sorbitol-transporting species, such as apple and peach, the loss of soluble acid invertase is not absolutely required for sugar accumulation (Li et al. 2012, Cirilli et al. 2016). Thus, sucrose hydrolysis by invertase plays only a minor role in sugar assimilation in pear fruits, which is consistent with the results of European pears (Reuscher et al. 2016).

Factors implicated in the regulation of sucrose content

TFs play critical roles in regulating sugar metabolism in plants (Ruan 2014), and several important TFs involved in sucrose accumulation have been identified, such as MYB (Wei et al. 2018), AP2/ERF (Xiao et al. 2016), bZIP (Frank et al. 2018) and WRKY (Sun et al. 2003). The role of some members of TF families, such as ERF2, which affect sucrose accumulation by mediating the expression of key genes involved in sucrose metabolism and hormone signaling pathways, is well known in Oryza sativa (Xiao et al. 2016). A new study shows that the higher expression level of PuHLS1 might cause a higher histone acetylation level of PuWRKY31, resulting in higher sucrose accumulation in the Ussurian pear fruit (Li et al. 2020). To better understand sucrose accumulation in pear fruits, we identified unique gene sets/modules associated with high-sucrose accumulation in both cultivars and performed gene coexpression network analysis. In our dataset, 42 TF-encoding genes that were highly correlated with sucrose content (r = 0.95, P = 1 × 10⁻⁹) (Fig. 5C) were found in the green module. It is proposed that a high-sucrose content is involved in the upregulation of C2H2, bZIP, GRAS, MADS and WRKY, which interact with specific binding sites on the promoters of target genes related to sucrose biosynthesis (Fig. 5D;Supplementary Table S6).

It has been reported that sucrose content is affected by multiple signals including light, hormones and stresses (Vassey 1989, Guy et al. 1992, Perez et al. 1997, Jia et al. 2011). In addition, ubiquitination and phosphorylation of sucrose transporter SUC2 might involve the regulation of sucrose metabolism (Xu et al. 2020). In this study, the promoters of different genes contain multiple cis-acting elements, such as light-responsive elements, low temperature-responsive elements and hormone response elements (Supplementary Table S7), and 20 genes involved in plant hormone signal transduction were clearly identified in the green module related to high-sucrose accumulation (Fig. 5B). We think that the sucrose accumulation in pear fruits is also regulated by plant hormones and other signals.

In conclusion, the comprehensive information presented here will serve as a robust resource for understanding sugar metabolism, particularly the synthesis, accumulation, and regulation of sucrose in pear fruits. A model was proposed to provide a mechanistic framework to understand the difference in sucrose content in pears and explain high-sucrose accumulation in HS3 (Fig. 7), in which sucrose deposited in the cell is resynthesized by SUS and SPS and high productive forces (signals → TFs → target genes) may promote flux through mainly sorbitol-derived fructose to enhance sucrose metabolism in pear fruits. It seems to be a very attractive means to change sugar composition by editing the target gene and/or regulating external factors in the pathway to achieve taste control. However, further functional validation of candidate genes and regulatory elements is required to elucidate a complete picture.

Fig. 7

A proposed model of sucrose content difference in pear fruits. Several types of TFs are involved in the upregulated transcription of FRK, PGI, SPS and SUS genes in the nucleus, and high sucrose then accumulates in vacuoles. Sucrose accumulation is affected by multiple signals, such as plant hormones. CINV, cytoplasmic invertase or neutral invertase; Glc, glucose; SDH, sorbitol dehydrogenase. Fig. 7

Materials and Methods

Plant materials

For the determination of sugar composition, 105 mature pear cultivars from five cultivated species, namely, P. pyrifolia, Pyrus bretschneideri, P. ussuriensis, P. communis and Pyrus sinkiangensis (Supplementary Table S1), were collected at the Jiangpu pear germplasm resource nursery of Nanjing Agricultural University, Nanjing, China. Two pear cultivars with contrasting sucrose contents, namely, ‘Hosui’ (HS-accumulating type) and ‘Korla’ fragrant (LS-accumulating type), were used in this study for transcriptome analysis and RNA-seq. Fruits of uniform size and free from visible defects were collected in three biological replicates at 30, 60, 90, 120 and 142 DAB for ‘Hosui’ pear and at 30, 60, 90, 120 and 150 DAB for ‘Korla’ fragrant pear. Approximately 10–20 pear fruits were collected for each biological replicate depending on the fruit size at different development stages. Finally, the fruits collected at the same sampling time for each cultivar were mixed after washing, and the peel and core were removed, the flesh were cut into pieces, immediately frozen in liquid nitrogen and stored at −80°C until further analysis.

Determination of sugar content

The content of sucrose, fructose, glucose and sorbitol was determined using HPLC as previously described (Cheng et al. 2018). Sugar was extracted from pear flesh by grounding, and pear flesh was then dissolved and filtered through a SEP-C18 cartridge (WAT021515; Waters, Shanghai, China) and Sep-Pak filter. A Waters 1525 HPLC system (Waters) was used to determine sugar content. The column had an inner diameter of 6.5 mm × 300 mm and a particle size of 10 μm (Waters), with a guard column from Sugar-pak 1 Guard-Pak Holder and Insert (Waters). The column temperature was maintained at 85°C. The injection volume and flow rate of the mobile phase were 10 μl and 0.6 ml min⁻¹, respectively. All sugar contents were determined according to an external standard solution. The concentration of each sample was calculated by the comparison of peak areas and retention times with those of calibrated sugar solutions of known concentrations.

RNA extraction and RNA-seq library construction and sequencing

Total RNA was extracted using a Plant RNA Isolation Kit (Auto Lab), followed by RNA purification using the RNeasy MiniElute Cleanup Kit (TIANGEN) according to the manufacturer’s instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer^® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using the Qubit^® RNA Assay Kit with the Qubit^® 2.0 Flurometer (Life Technologies, CA, USA). The RNA integrity number (RIN) was determined using the RNA 6000 Nano Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction and RNA-seq were performed by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). Briefly, high-quality total RNA (RIN ≥ 8) of three biological replicates of six samples (18 in total) was processed using a TrueSeq RNA Sample Prep kit (Illumina Technologies) for library construction. mRNA was enriched and purified from total RNA using oligo(dT)-rich magnetic beads and then cleaved into short fragments. Using these cleaved mRNA fragments as templates, first- and second-strand cDNA samples were synthesized. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation at 3′ ends of DNA fragments, the NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments of preferentially 150–200 bp in length, library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Libraries were sequenced using Illumina HiSeq™ 2000.

Sequence data processing and mapping reads to the pear genome

Raw data (raw reads) in the FASTQ format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing an adapter, reads containing ploy-N and reads of low quality from the raw data. At the same time, Q20, Q30 and GC contents in the clean data were calculated. All downstream analyses were based on the high-quality clean data and mapped to the assembled pear genome data (Wu et al. 2013). An index of the reference genome was built using Bowtie v2.2.3, and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12 (Trapnell et al. 2009). Transcripts that did not exist in the CDS database of the pear genome were extracted to predict new genes using the EMBOSS package (http://emboss.open-bio.org/). HTSeq v0.6.1 was used to count the read numbers mapped to each gene; then, the FPKM of each gene was calculated based on the length of a gene and the read counts mapped to it.

DEGs, functional annotation and pathway enrichment analysis

Differential gene expression analysis of samples was performed using the DESeq R package (1.18.0) based on clean-read counts. The resulting P-values were adjusted using Benjamini and Hochberg’s (2000) approach for controlling the false discovery rate. Hierarchical clustering was performed to identify DEGs. DEG screening was conducted by setting the probability of an adjusted P-value of <0.001, absolute log₂ fold change of ≥2 or ≤−2 and gene with the FPKM of >1 in at least one sample.

To investigate the function of DEGs, GO enrichment, a standard international gene function classification system based on molecular functions, cellular components and biological processes was used for further analysis (Conesa et al. 2005). GO terms with the corrected P-values of ≤0.05 were considered significantly enriched. Annotation was analyzed using KOBAS (v2.0) (http://kobas.cbi.pku.edu.cn/), which detects significantly enriched pathways using a hypergeometric test and has been effectively used for the differential pathway analysis of living entities, such as plants, animals and bacteria (Wu et al. 2006). All DEGs mapped in KEGG pathways with the P-values of ≤0.05 were considered significantly enriched.

qRT-PCR validation of DEGs

In total, 12 DEGs were randomly selected for verification by qRT-PCR. Their coding sequences were acquired from the pear genome (http://peargenome.njau.edu.cn). Primers specific to the 12 DEGs were designed (Supplementary Table S2) using Primer Premier 5.0. The LightCycler 480 SYBR Green Master (Roche, USA) system was used, with tubulin (Pbr042345.1) as the housekeeping gene. Three biological replicates were processed for qRT-PCR, which was performed using a 20-μl reaction mixture containing 80–100 ng of cDNAs, 200 nM of each primer and 10 μl of LightCycler 480 SYBR Green I Master Mix (Roche, Basel, Switzerland). All reactions were run as duplicates on 96-well plates. PCRs were performed under the following conditions: preincubation at 95°C for 5 min; 45 cycles of 95°C for 3 s, 60°C for 10 s and 72°C for 30 s; and a final extension at 72°C for 3 min. Expression levels were calculated using the 2^−ΔΔCt method (Livak and Schmittgen 2001) for each sample.

Gene coexpression network

A WGCNA package was used to find clusters (modules) of highly correlated genes and to construct a coexpression network (Langfelder and Horvath 2008). Log base 2 (1 + FPKM) values with a matrix of pairwise Pearson correlations between all pairs of genes were generated and further transformed into a weighted matrix with a soft-threshold power of 13, which gave greater weight to the strongest correlations while maintaining gene–gene connectivity. To decrease data noise, adjacency matrices were transformed into topological overlap matrices (Yip and Horvath 2007). Next, we used hierarchical clustering with the dynamic tree cut algorithm to generate a tree to identify similar modules. A minimum module size of 30 and a height cutoff of 0.25, corresponding to a correlation of 0.8, were used to merge similar transcripts. To determine the association of a module with stage-specific expression for each cultivar, we determined the correlation between each ME and the binary indicator (tissue/stage = 1 and all other samples = 0) as described previously (Downs et al. 2013). A positive correlation indicated that the genes in a module had higher/preferential expression at a particular stage relative to all other samples. Furthermore, we followed a cross-tabulation approach to create a contingency table that reported the number of genes that fell into modules. Finally, networks were visualized using Cytoscape v.3.0.0.

TF analysis

Promoters (2,000 bp upstream) were obtained from the pear genome (http://peargenome.njau.edu.cn/), and TF-binding motifs were aligned to Arabidopsis using PlantPAN 2.0 (http://plantpan2.itps.ncku.edu.tw/) and PlantRegMAP (http://plantregmap.cbi.pku.edu.cn/binding_site_prediction.php). We analyzed the DEGs determined by the annotated databases to mine for TFs using a Linux server; then, data were further aligned with the protein sequences of the plant TF database PlantTFDB (version 3.0). Plant TFs were obtained from the TF database website (http://plntfdb.bio.uni-potsdam.de/v3.0/), and we then analyzed the annotated TFs with an e-value of ≤10⁻¹⁰ using BLASTX. Finally, annotated TFs from pear fruits were analyzed and classified into different groups based on the standards of PlantTFDB (version 3.0) (Pérez-Rodríguez et al. 2010). We then used PlantCARE, a plant cis-acting regulatory element database, to analyze promoters (Rombauts et al. 1999).

Data analysis

Data processing was performed using R package and Microsoft Excel. For linear regression normalization, PCA, Q–Q plot analysis and hierarchical clustering were performed. A Venn diagram was constructed using the Venn diagram and R packages. Mean values and standard deviations were obtained from three biological replicates.

Supplementary Data

Supplementary data are available at PCP online.

Funding

Fundamental Research Funds for the Central Universities in China (KYZ201510, KYZ201834).

Disclosures

The authors have no conflicts of interest to declare.

Author notes

Jiahong Lü and Xin Tao contributed equally.