Total flavonoid content
The total flavonoid content in apple fruit was determined using the aluminum chloride colorimetric method, and the results showed that the total flavonoid content in BG and BGBB was significantly higher than that in JKG, with ratios of 4.28 (BG/JKG) and 4.68 (BGBB/JKG), respectively, while the total flavonoid content in JKBF was lower than that in JKG, with a ratio of 0.57 (JKBF/JKG) (Fig. 1). At the same time, we found the total flavonoid content in BGBB was a little higher than that in BG, with ratio of 1.09 (BGBB/BG), and the total flavonoid content in JKBF was significantly lower than that in BG and BGBB, with ratios of 7.55 (BG/JKBF) and 8.26 (BGBB/JKBF), respectively.
RNA-Seq and quality filtering
Healthy and bitter-pit apple fruit of uniform ripeness and size were chosen to study the influence of bitter pit on different expression of flavonoid biosynthesis-related genes. Twelve cDNA libraries of the whole transcriptome were generated for RNA-SEq. We obtained 963.948 M raw reads containing 144.592 G nucleotides by RNA-SEq. After clean-up and quality filtering, 862.772 M clean reads containing 124.435 G nucleotides were obtained (Table 1). The clean Q30 percentages were 94.1%, indicating that the RNA-seq results were of high quality and suitable for use in further analyses. By miRNA sequencing, 4331 reads associated with 3753 target genes were obtained, among which 2571 genes were upregulated and 1182 genes were downregulated. Of the 41894 total RNA sequences subjected to sequence analysis and annotation, 15852 were known, and 26042 were novel (Table 2).
Table 1
The quality control of whole-transcriptome RNA-Seq stats
Quantity (Percentage) | Raw Reads (M) | Raw Bases (G) | Raw Q20 (G) | Raw Q30 (G) | Clean Reads (M) | Clean Bases (G) | Clean Q20 (G) | Clean Q30 (G) | Average Length (bp) |
JKG_1 | 93.383 | 14.007 | 13.308 (95.0%) | 12.528 (89.4%) | 82.488 (88.3%) | 11.988 (85.6%) | 11.787 (98.3%) | 11.388 (95.0%) | 145.3 |
JKG_2 | 75.388 | 11.308 | 10.721 (94.8%) | 10.095 (89.3%) | 66.332 (88.0%) | 9.559 (84.5%) | 9.412 (98.5%) | 9.108 (95.3%) | 144.1 |
JKG_3 | 88.231 | 13.235 | 12.573 (95.0%) | 11.856 (89.6%) | 77.882 (88.3%) | 11.260 (85.1%) | 11.084 (98.4%) | 10.730 (95.3%) | 144.6 |
BG_1 | 81.606 | 12.241 | 11.703 (95.6%) | 11.063 (90.4%) | 73.073 (89.5%) | 10.652 (87.0%) | 10.481 (98.4%) | 10.140 (95.2%) | 145.8 |
BG_2 | 63.588 | 9.538 | 9.133 (95.8%) | 8.649 (90.7%) | 57.173 (89.9%) | 8.378 (87.8%) | 8.250 (98.5%) | 7.991 (95.4%) | 146.5 |
BG_3 | 80.519 | 12.078 | 11.575 (95.8%) | 10.966 (90.8%) | 72.705 (90.3%) | 10.596 (87.7%) | 10.433 (98.5%) | 10.103 (95.3%) | 145.7 |
BGBB_1 | 77.476 | 11.621 | 11.026 (94.9%) | 10.333 (88.9%) | 70.019 (90.4%) | 10.038 (86.4%) | 9.857 (98.2%) | 9.485 (94.5%) | 143.4 |
BGBB_2 | 78.607 | 11.791 | 11.187 (94.9%) | 10.464 (88.7%) | 69.857 (88.9%) | 10.070 (85.4%) | 9.871 (98.0%) | 9.485 (94.2%) | 144.1 |
BGBB_3 | 88.091 | 13.214 | 12.487 (94.5%) | 11.638 (88.1%) | 79.116 (89.8%) | 11.347 (85.9%) | 11.130 (98.1%) | 10.681 (94.1%) | 143.4 |
JKBF_1 | 76.206 | 11.431 | 10.844 (94.9%) | 10.159 (88.9%) | 68.899 (90.4%) | 9.927 (86.8%) | 9.747 (98.2%) | 9.378 (94.5%) | 144.1 |
JKBF_2 | 84.28 | 12.642 | 12.013 (95.0%) | 11.258 (89.0%) | 75.281 (89.3%) | 10.809 (85.5%) | 10.603 (98.1%) | 10.199 (94.4%) | 143.6 |
JKBF_3 | 76.573 | 11.486 | 10.984 (95.6%) | 10.359 (90.2%) | 69.947 (91.3%) | 9.811 (85.4%) | 9.638 (98.2%) | 9.297 (94.8%) | 140.3 |
Table 2
Transcript types of RNA-Seq stats
Class | mRNA | lncRNA | circRNA | miRNA | Others | Unknown | Total |
Known | 12645 | 0 | 52 | 3017 | 0 | 138 | 15852 |
Novel | 20937 | 1860 | 46 | 2877 | 322 | 0 | 26042 |
Total | 33582 | 1860 | 98 | 5894 | 322 | 138 | 41894 |
The RNA-Seq data were submitted to SRA database in NCBI with the accession number of PRJNA640254.
Analysis of sequence length distribution
The sequence length distribution of the transcripts is depicted by FPKM values in Table 3. FPKM was used to quantify gene expression, as this method eliminates the effect of different gene lengths and sequencing levels on the gene expression calculation. For the fruit samples of JKG, BG, BGBB and JKBF, there were 17086, 16396, 16736 and 17401 transcripts with an FPKM value < 5.0, accounting for 48.17%, 45.80%, 46.51% and 49.17% of the total, respectively. There were 3874, 4044, 4178 and 3688 transcripts with FPKM values > 50.0, accounting for 10.92%, 11.3%, 11.61% and 10.42% of the total, respectively.
Table 3
Sequence length distribution of transcripts and genes assembled from Illumina reads
FPKM Value (Percentage) | 0-0.5 | 0.5-1.0 | 1.0–5.0 | 5.0–10.0 | 10.0–50.0 | > 50.0 | Total |
Transcripts | JKG | 3402 (9.59%) | 2866 (8.08%) | 10818 (30.50%) | 5356 (15.10%) | 9156 (25.81%) | 3874 (10.92%) | 35472 |
BG | 2354 (6.58%) | 2801 (7.82%) | 11241 (31.40%) | 5664 (15.82%) | 9694 (27.08%) | 4044 (11.30%) | 35798 |
BGBB | 1363 (3.79%) | 3090 (8.59%) | 12283 (34.14%) | 5563 (15.46%) | 9506 (26.42%) | 4178 (11.61%) | 35983 |
JKBF | 3493 (9.87%) | 2868 (8.1%) | 11040 (31.2%) | 5304 (14.99%) | 8996 (25.42%) | 3688 (10.42%) | 35389 |
Genes | JKG | 1580 (6.00%) | 1348 (5.12%) | 6454 (24.52%) | 4282 (16.27%) | 8694 (33.03%) | 3963 (15.06%) | 26321 |
BG | 816 (3.07%) | 1208 (4.55%) | 6646 (25.01%) | 4554 (17.14%) | 9254 (34.82%) | 4099 (15.42%) | 26577 |
BGBB | 337 (1.27%) | 1214 (4.56%) | 7196 (27.03%) | 4519 (16.97%) | 9141 (34.33%) | 4219 (15.85%) | 26626 |
JKBF | 1699 (6.46%) | 1356 (5.16%) | 6536 (24.86%) | 4326 (16.46%) | 8604 (32.73%) | 3768 (14.33%) | 26289 |
Similar results were obtained for the sequence-length distribution of genes (Table 3). For the fruit samples of JKG, BG, BGBB and JKBF, there were 9382, 8670, 8747 and 9591 genes with FPKM values < 5.0, accounting for 35.64%, 32.62%, 32.85% and 36.48% of the total, respectively. There were 3963, 4099, 4219 and 3768 genes with FPKM values > 50.0, accounting for 15.06%, 15.42%, 15.85% and 14.33% of the total, respectively.
Analysis of DETs and DEGs
The expression patterns of transcripts in healthy and bitter-pit apple fruits were investigated using pairwise comparisons. By whole-transcriptome sequencing, 35920, 35990 and 35714 transcripts were obtained from JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF, respectively, among which 2577, 4081 and 776 transcripts, respectively, were differentially expressed. Compared with JKG, BG, BGBB and JKBF had 2085, 2774 and 325 upregulated differentially expressed transcripts (DETs), respectively, and 492, 1307 and 154 downregulated DETs, respectively (Table 4).
Table 4
The DETs and DEGs obtained from JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF
Treatment | DETs | DEGs | miRNA |
Total | Up | Down | Total | Up | Down | Total | Up | Down |
JKG vs. BG | 2577 | 2085 | 492 | 2498 | 2035 | 463 | 134 | 45 | 89 |
JKG vs. BGBB | 4081 | 2774 | 1307 | 3847 | 2680 | 1167 | 157 | 56 | 101 |
JKG vs. JKBF | 776 | 325 | 451 | 774 | 337 | 437 | 137 | 34 | 103 |
Correspondingly, 26600 shared genes were obtained from JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF, among which 2498, 3847 and 774 genes, respectively, were differentially expressed (Table 4). Compared with JKG, BG, BGBB and JKBF had 2035, 2680 and 337 upregulated DEGs, respectively, and 463, 1167 and 437 downregulated DEGs, respectively. For miRNA sequencing, 134, 157 and 137 miRNAs DEGs were obtained from JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF, respectively (Table 4). Compared with JKG, BG, BGBB and JKBF had 45, 56 and 34 upregulated DEGs, respectively, and 89, 101 and 103 downregulated DEGs, respectively.
In summary, a total of 2632, 4004 and 911 DEGs were obtained from JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF, respectively, of which 2080, 2736 and 371 DEGs were upregulated and 552, 1268 and 540 DEGs were downregulated.
GO enrichment analysis of DEGs
For the JKG vs. BG comparison, a total of 4643 genes (3922 upregulated and 721 downregulated) were highly enriched in 3 classes and 34 sub-classes, of which 1065, 2632 and 946 genes were annotated with GO terms related to “molecular function”, “biological process” and “cellular component”, respectively (Fig. 2a). For the JKG vs. BGBB comparison, a total of 4760 genes (3121 upregulated and 1639 downregulated) were highly enriched in 3 classes and 32 sub-classes, of which 1357, 2605 and 798 genes were annotated with GO terms related to “molecular function”, “biological process” and “cellular component”, respectively (Fig. 2b). For the JKG vs. BG comparison, a total of 2156 genes (501 upregulated and 1655 downregulated) were highly enriched in 3 classes and 34 sub-classes, of which 457, 1107 and 592 genes were annotated with GO terms related to “molecular function”, “biological process” and “cellular component”, respectively (Fig. 2c).
KEGG enrichment analysis of DEGs
In the KEGG pathway enrichment analysis, 60 pathways corresponding to 646 DEGs, 73 pathways corresponding to 788 DEGs and 31 pathways corresponding to 109 DEGs were significantly enriched in the JKG vs. BG, JKG vs. BGBB and JKG vs. JKBF comparisons, respectively, of which 26 DEGs (0 downregulated and 24 upregulated), 23 DEGs (0 downregulated and 23 upregulated) and 3 DEGs (3 downregulated and 0 upregulated) involved in flavonoid biosynthesis were enriched (Fig. 3).
The flavonoid biosynthesis pathway in KEGG database from JKG vs. BG (Fig. 4a), JKG vs. BGBB (Fig. 4b) and JKG vs. JKBF (Fig. 4c) indicated that the DEGs involved in flavonoid biosynthesis were upregulated in JKG vs. BG and JKG vs. BGBB while downregulated in JKG vs. JKBF, which was consistent with the KEGG pathway enrichment analysis.
qRT-PCR Validation and expression analysis of DEGs
To confirm the reproducibility and accuracy of the differential gene expression identified through Illumina analysis, 8 DEGs [CYP98A2 (1), CYP98A2 (2), BADH, DAT, MdHCT (1), MdHCT (2), CHI (1) and CHI (2)] related to flavonoid biosynthesis [ko00941] with fold-changes ≥ 1 were selected for qRT-PCR validation (Table 5). The expression patterns of these genes according to RNA-Seq and qRT-PCR are shown in Table 6.
Table 5
The qRT-PCR primers for the 8 selected genes
Gene ID | Gene Name | Predicted Function | Primers (5′-3′) |
MD08G1242900 | CYP98A2 (1) | cytochrome P450 98A2-like (LOC103441965) | F: AACCACTGCACCAAACCTGA R: AGCGATCCGCCTATCTTCAC |
MD15G1436500 | CYP98A2 (2) | cytochrome P450 98A2-like (LOC103441965) | F: CGGTATCAACTTGGTGGCCT R: CGTCCCTGGATTTTCGGACA |
MD05G1219000 | BAHD | BAHD acyltransferase At5g47980-like | F: GGCGGATGTTTCCACACTCA R: CAGCCCCAAATGCTACGAGA |
MD16G1110600 | DAT | deacetylvindoline O-acetyltransferase-like | F: CAGTGAACCTGCGCACTAGA R: GCAACTCAATCACGCTGTCC |
MD17G1224900 | MdHCT (1) | Hydroxycinnamoyl-CoA Shikimate O-hydroxycinnamoyl transferase-like | F: AATAAGACCGCACAGACGGA R: GTCGAAGAAGGAGTTGCCCT |
MD17G1225100 | MdHCT (2) | Hydroxycinnamoyl-CoA shikimate O-hydroxycinnamoyl transferase-like | F: GGTGATTCCGAGCATCCACA R: GCCGCAGTCAATCTCGATAC |
MD01G1118000 | CHI (1) | chalcone–flavonone isomerase | F: CAGAGAAGCGGGTCTCTCCC R: GCGTTGGAGCCATTTGACAAT |
MD01G1118100 | CHI (2) | chalcone–flavonone isomerase | F: ATTTCCCACCCGGTGCTT R: CCGGCTTCAGCACCATTAGA |
| MdActin | endogenous control | F: TGACCGAATGAGCAAGGAAATTACT R: CTCAGCTTTGGCAATCCACATC |
Table 6
RNA-Seq FPKM values, qRT-PCR relative expression levels and ratios of 8 candidate DEGs
Genes Values | CYP98A2 (1) | CYP98A2 (2) | BADH | DAT | MdHCT (1) | MdHCT (2) | CHI (1) | CHI (2) |
RNA- seq | JKB | 61.80 | 2.32 | 0.37 | 5.67 | 60.45 | 4.04 | 10.29 | 0.71 |
BG | 88.61 | 25.19 | 4.25 | 9.66 | 405.70 | 9.13 | 69.70 | 20.83 |
BGBB | 112.91 | 42.38 | 6.81 | 18.12 | 780.75 | 16.20 | 104.78 | 36.87 |
JKBF | 27.55 | 1.62 | 0.30 | 1.94 | 31.10 | 1.82 | 7.95 | 0.46 |
BG/JKG | 1.43 | 10.86 | 11.48 | 1.70 | 6.71 | 2.26 | 6.77 | 29.42 |
BGBB/JKG | 1.83 | 18.27 | 18.42 | 3.20 | 12.92 | 4.01 | 10.18 | 52.08 |
JKBF/JKG | 0.45 | 0.70 | 0.82 | 0.34 | 0.51 | 0.45 | 0.77 | 0.65 |
qRT- PCR | JKB | 1.58 | 1.71 | 5.33 | 1.27 | 5.29 | 2.08 | 4.76 | 1.83 |
BG | 1.80 | 6.41 | 35.60 | 7.61 | 17.97 | 2.18 | 17.26 | 13.14 |
BGBB | 1.98 | 6.90 | 38.04 | 8.65 | 24.92 | 2.23 | 18.29 | 14.71 |
JKBF | 0.65 | 0.20 | 5.13 | 0.85 | 4.81 | 0.49 | 1.83 | 1.37 |
BG/JKG | 1.14 | 3.74 | 6.69 | 6.01 | 3.40 | 1.05 | 3.62 | 7.17 |
BGBB/JKG | 1.25 | 4.03 | 7.14 | 6.83 | 4.71 | 1.07 | 3.84 | 8.03 |
JKBF/JKG | 0.41 | 0.12 | 0.96 | 0.67 | 0.91 | 0.24 | 0.38 | 0.74 |
The expression patterns of the 8 candidate genes were determined by RNA-Seq and analyzed to determine the FPKM values. The results showed that these genes were upregulated in JKG vs. BG and JKG vs. BGBB, with FPKM ratios of 1.43–29.42 and 1.83–52.08, respectively, while they were downregulated in JKG vs. JKBF, with FPKM ratios of 0.34–0.82. The expression levels of the 8 candidate genes were validated by qRT-PCR analysis, and the results showed that the relative expression levels of these genes were increased in JKG vs. BG and JKG vs. BGBB comparisons, with ratios of 1.05–7.17 and 1.07–8.03, respectively, while the relative expression levels were decreased in JKG vs. JKBF comparison, with expression ratios of 0.12–0.96.