Identification of candidate genes influencing anthocyanin biosynthesis during the development and ripening of red and white strawberry fruits via comparative transcriptome analysis

Plant materials

Fruits from the three octoploid cultivated strawberry cultivars (‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’) were used in the present study (Figs. 1A–1I). Plantlets of the three cultivars were grown in a greenhouse at temperatures ranging from 8 °C (dark) to 28 °C (light) and a relative humidity ranging from 55–70%. All the experiments were performed in the Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences (Zhengzhou, Henan, China). The fruit development and ripening process was divided into three visual stages: S2, middle green fruit; S5, initial red fruit; S7, full red fruit (Figs. 1G–1I), as described in previous reports (Fait et al., 2008; Zhao et al., 2018). Three replicates experiments were performed, with each consisting of 30 fruits. Approximately 300 g of fruits was randomly sampled and pooled, constituting one replicate for subsequent analysis (RNA-seq, anthocyanin content, sugar content, and qRT-PCR). The sample from each stage was ground into a powder with liquid nitrogen and stored at −80 °C for later use in each experiment.

RNA extraction, library construction and RNA sequencing

According to a previous study (Estrada-Johnson et al., 2017; Hu et al., 2018), the development and ripening process of three strawberry cultivars was divided into three key important stages (S2, S5, and S7) (Figs. 1G–1I) and were chosen to examine the content of the major anthocyanins (Figs. 1J–1O) and sugars. For transcriptome analysis, three key stages (S2, S5, and S7) from the three strawberry cultivars were chosen to study fruit coloration in the present study. Total RNA was isolated from the fruit powders using an E.Z.N.A Plant RNA Kit (R6827-01, Omega, USA) following the manufacturer’s instructions. The concentration and quality of total RNA were analysed on a NanoDrop 2000 spectrophotometer (Thermo Fisher, USA), and RNA degradation and contamination were assessed using 1% agarose gels. Two biological replicates were performed per group. For each sample, 2 µg of total RNA was used for library construction and subjected to deep sequencing on an Illumina HiSeq X Ten platform (BerryGenomics Corporation, Beijing). All sequencing data have been submitted to the NCBI Sequence Read Archive (SRA accession number PRJNA552213).

Transcriptome and transcript analysis

Reads from each library were assembled separately. Raw sequence reads were filtered using the Illumina pipeline according to the methods of BerryGenomics, Beijing, China. Briefly, after removing the adaptor sequences and low-quality sequences (including the reads with percentages of N over 10%), the remaining high-quality clean reads were used for analysis (Tables S1–S3). The clean sequence data was compared to the Fragaria vesca diploid strawberry genome (National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/genome/3314), abbreviated as in the diploid genome) and the Fragaria x ananassa Camarosa genome (the genome database for rosaceae (GDR, https://www.rosaceae.org/species/fragaria_x_ananassa/genome_v1.0.a1), abbreviated as in the octoploid genome) to identify all DEGs in S2 vs. S5, S2 vs. S7, S5 vs. S7 stage comparisons for the ‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’ cultivars, respectively.

The programs Cuffquant and Cuffnorm used fragments per kilobase of transcript per million fragments mapped (FRKM) as a measure of transcript or gene expression levels (Mortazavi et al., 2008). DEseq2 was used to identify DEGs from two samples, and the results of all statistical tests were adjusted with a fold change ≥ 2 and false discovery rate (FDR) <0.01 (Love, Huber & Anders, 2014). Venn diagrams for the different DEGs between each combination compared to the diploid strawberry genome (Benihoppe vs. Xiaobai, Benihoppe vs. Snow White, and Xiaobai vs. Snow White) (Figs. 2A–2C). Three stages (S2, S5, and S7) of Benihoppe, Xiaobai, Snow White (Figs. 2D–2F) were evaluated. A number of common and specific DEGs were identified in the six assayed combinations compared to the diploid strawberry genome (Fig. 2G). Venn diagrams were generated for the different DEGs between each combination compared to the octoploid strawberry genome (Figs. 2H–2N). The significance of DEGs was determined by determining the FDR adjusted p-value. A log₂ value (one sample/one sample) of >1 or <-1 was the criterion used to select candidate genes (anthocyanin biosynthesis pathway (Fig. 3), anthocyanin, sugar and hormone -related pathways (Fig. 4) and TFs (Fig. 5) (Tables S4, S5). Genes were annotated according to BLAST search results that were compared to sequences in several databases, including the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), COG/KOG, Pfam, and NCBI non-redundant protein sequences (Nr) databases, as well as a manually annotated and reviewed section of the UniProt Knowledgebase database (Swiss-Prot). In addition, a heatmap was drawn using TBtools.

Validation of gene expression profiles by qRT-PCR

To test the reliability of the RNA-seq data, nine candidate genes were selected for qRT-PCR expression analysis. Total RNA was extracted from three strawberry cultivars at three developmental stages (S2, S5, and S7) using an E.Z.N.A. Plant RNA Kit following the manufacturer’s instructions. The concentration and quality of RNA were analysed on a NanoDrop 1000. Approximately 1.5 µg of total RNA was used for cDNA synthesis using a PrimerScriptTM RT Reagent Kit with gDNA Eraser (TaKaRa, China). All the primers used in the qRT-PCR analysis were designed using Vector NTI (Table S6) without any interference of the conserved region, with an amplified product length of 150–300 bp. The cDNA concentration was adjusted based on the strawberry housekeeping gene FvRib413 (Meng et al., 2018; Su, Wu & Cui, 2016). For qRT-PCR, each 20 µl reaction contained 10 µl of 2 × LightCycler 480 SYBR Green I Master mix (Cat# 4887352001, Roche), 2 µl of 50 × diluted cDNA, 0.4 µl of each primer, and 7.2 µl of ddH₂O, and the data were processed with LightCycler^® 480 software (Roche, China). The cycling program was as follows: 95 °C for 5 min followed by 45 cycles of 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 30 s, which was followed by a melting curve analysis at 60−95 °C and obtained one melting point then subsequently used for further study. The relative expression level of each gene was calculated using the comparative 2^−ΔΔCT method (Livak & Schmittgen, 2001). All analyses were repeated three times using biological replicates.

Measurement of anthocyanins and soluble sugars

For anthocyanin extraction, 5 g of frozen powders of fruits from three stages (S2, S5, and S7) for the three strawberry cultivars were used in the present study (Table S7). To extract the anthocyanin mother liquid, an ultrasonic extraction process was performed for 30 min, and the mixture was heated in a boiling-water bath for 1 h. Subsequently, the mixture was filtered through 0.75-µm filter and analysed by high-performance liquid chromatography. The absorbance was measured at 530 nm, and the solvent flow rate was 0.8 ml/min. Three biological replicates were performed for each analysis.

Subsequently, 5 g of frozen powders of fruits from three stages (S2, S5, and S7) for the three strawberry cultivars were used for soluble sugar analysis (Table S8). The powder was added to deionized water, after which the volume was brought up to 50 ml, and the solution was incubated at 75 °C for 30 min. Then, the mixture was centrifuged at 12,000 rpm for 10 min and then filtered through a Dionex OnGuard II column to acquire a liquid supernatant, which was subsequently analysed by high-performance liquid chromatography. The methods used to determine the levels of anthocyanins and soluble sugars were provided by the Institute of pomology of the Chinese Academy of Agricultural Sciences. All analyses were repeated three times using biological replicates.

Statistical analysis

Each experiment was independently repeated at least three times. The data were statistically analysed with a one-way analysis of variance of Duncan’s multiple range text using IBM SPSS Statistics 20. The results are presented as the mean value ± standard deviation of the mean (SD) (Figs. 1J–1O; Fig. 3; Figs. 4A–4C; Fig. 6), and significant differences relative to controls are shown in lowercase (<0.05) and capital (<0.01) letters (Figs. 1J–1O; Fig. 3; Figs. 4A–4C; Fig. 6).

Anthocyanin accumulation in fruits of three strawberry cultivars at different developmental stages

In the present study, the anthocyanin composition and contents of the red-fruited strawberry (‘Benihoppe’), the white-fleshed mutant of ‘Benihoppe’ (‘Xiaobai’), and the white-fruited strawberry (‘Snow White’) were analysed (Figs. 1A–1I). As expected, cyanidin (cyanidin-3-O-glucoside and cyanidin-3-malonylglucosyl-5-glucoside) and pelargonidin (pelargonidin-3-O-glucoside, pelargonidin-3-O-rutinoside, and pelargonidin-3-malonylglucosyl-5-glucoside) were the two major types of anthocyanins present in the strawberry fruits, which contributed to the red colour of the strawberry fruit (Figs. 1J–1O). Generally, the cyanidin and pelargonidin contents increased with the fruit development and ripening stages in the three cultivars. The contents of pelargonidin (167.5 mg/kg) and cyanidin (8.2 mg/kg) in ‘Benihoppe’ were significantly higher than those detected in ‘Xiaobai’ (13.3 mg/kg of pelargonidin and 2.1 mg/kg of cyanidin), while smaller amounts of cyanidin-3-O-glucoside (0.2 mg/kg) and pelargonidin-3-O-glucoside (0.7 mg/kg) were detected in ‘Snow White’ (Figs. 1J–1O, Table S7). Significant differences in the quantity of cyanidin and pelargonidin were observed between the red and white strawberry fruits, indicating the differential metabolism of anthocyanin in these cultivars.

RNA-seq and de novo transcriptome assembly

Eighteen cDNA libraries from 3 stages (S2, S5, and S7) of ‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’ with two repeats were used for sequence analysis with the Illumina HiSeq platform. These libraries were subjected to RNA-seq using an Illumina HiSeq X Ten, generating 27,673,862, 33,937,779, 27,331,497, 28,914,712, 32,588,549, 26,390,260, 25,061,869, 29,690,765, 25,703,115, 26,179,018, 26,487,231, 28,493,685, 25,410,599, 31,471,955, 30,928,839, 28,521,409, 28,248,711, and 25,042,090 150-bp paired-end raw reads (Table S1). All of the raw reads are available in the NCBI SRA database (accession number PRJNA552213). On average, 91.2% of the reads from the 18 libraries had a quality score over Q30. Moreover, over 99% and approximately 90% of the cleaned reads could be uniquely mapped to the diploid and octoploid strawberry genomes, respectively. Pearson’s correlation coefficients of the transcriptome profiles were 0.99 between each set of biological replicates, indicating that the sequencing quality was suitable for further analysis (Tables S2, S3).

The identification of differentially expressed genes (DEGs) between genotypes and developmental stages

Using the diploid strawberry genome database, 10,354 DEGs were identified by comparing our transcriptome data with the diploid strawberry genome. A Venn diagram showing the DEGs from the three strawberry cultivars (in the diploid and octoploid genomes) was drawn using TBtools (Fig. 2). In total, 1817, 4543, and 4767 DEGs were identified for the three stages (S2, S5, and S7) in comparisons between ‘Benihoppe’ and ‘Xiaobai’ (Fig. 2A), ‘Benihoppe’ and ‘Snow White’ (Fig. 2B), and ‘Snow White’ and ‘Xiaobai’ (Fig. 2C), respectively. In addition, 4836, 4773, and 6245 DEGs were identified among the three stages (S2 vs. S5, S2 vs. S7, and S5 vs. S7) for ‘Benihoppe’ (Fig. 2D), ‘Xiaobai’ (Fig. 2E), and ‘Snow White’ (Fig. 2F), respectively. Moreover, 5049 DEGs were identified between S5 and S7 stages when comparing the Benihoppe S5 vs. Xiaobai S5, Benihoppe S5 vs. Snow White S5, Benihoppe S5 vs. Xiaobai S5, Benihoppe S7 vs. Xiaobai S7, Benihoppe S7 vs. Snow White S7, and Benihoppe S7 vs. Xiaobai S7 samples (Fig. 2G). In addition, 53,200 DEGs were obtained by comparing our transcriptome data with the octoploid strawberry genome to identify candidate genes, with 33,397, 46672, and 46,610 DEGs identified for the three stages (S2, S5, and S7) in comparisons between ‘Benihoppe’ and ‘Xiaobai’ (Fig. 2H), ‘Benihoppe’ and ‘Snow White’ (Fig. 2I), and ‘Snow White’ and ‘Xiaobai’ (Fig. 2J), respectively. Furthermore, 37,242, 37,529, and 41,867 DEGs were identified when comparing three stages (S2 vs. S5, S2 vs. S7, and S5 vs. S7) for ‘Benihoppe’ (Fig. 2K), ‘Xiaobai’ (Fig. 2L), and ‘Snow White’ (Fig. 2M), respectively. Moreover, 53,200 DEGs were identified between the S5 and S7 stages for the Benihoppe S5 vs. Xiaobai S5, Benihoppe S5 vs. Snow White S5, Benihoppe S5 vs. Xiaobai S5, Benihoppe S7 vs. Xiaobai S7, Benihoppe S7 vs. Snow White S7, and Benihoppe S7 vs. Xiaobai S7 sample comparisons (Fig. 2N). The number of DEGs identified in the octoploid genome was 10-fold higher compared to that observed in the diploid genome between the S5 and S7 stages among the three cultivars (Fig. 2). In addition, the greatest changes occurred between the S2 vs. S5, and S2 vs. S7 stages, while the smallest changes occurred between the S5 vs. S7 stages in the three octoploid strawberry cultivars. The number of DEGs in the octoploid genome was approximately 10-fold higher than that observed in the diploid genome. These results demonstrated that the genetics of ‘Benihoppe’ and ‘Xiaobai’ are similar, and relatively few DEGs were involved in the fruit ripening process compared to the assayed other combinations.

Genes related to the anthocyanin biosynthesis pathway in the fruit development and ripening process

To investigate which anthocyanin synthesis step was blocked in ‘Xiaobai’ and ‘Snow White’, we compared the expression of anthocyanin-related genes in three developmental stages (S2, S5, and S7) of ‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’ to test whether the expression of some of anthocyanin-related genes affect fruit colour. Twelve DEGs were identified in the anthocyanin biosynthesis pathway among the three cultivars in the diploid genome (as mentioned above) (Fig. 3A, Tables S4, S5), of which eleven showed significantly higher expression levels in the fruits of ‘Benihoppe’ and ‘Xiaobai’ compared to those of ‘Snow White’. In addition, only one DEG involved in the anthocyanin biosynthesis pathway was upregulated in the fruits of ‘Snow White’ compared to those of ‘Benihoppe’ and ‘Xiaobai’ (Fig. 3A, Tables S4, S5). Forty DEGs were identified in the anthocyanin biosynthesis pathway among the three cultivars in the octoploid strawberry (as mentioned above) (Fig. 3B, Tables S4, S5), of which 24 involved in the anthocyanin biosynthesis pathway were downregulated in the fruits of ‘Snow White’ compared to those of ‘Benihoppe’ and ‘Xiaobai’. Moreover, the expression levels of four structural genes in ‘Xiaobai’ were lower than these in ‘Benihoppe’ in the diploid genome, including PAL (gene 22493), F3H (gene 1350), DFR (gene 21883, gene 236, gene23371, and gene 7069) and 3-GT (gene 25816 and gene 7173). One DEG (gene8111) involved in the anthocyanin biosynthesis pathway was upregulated in the fruits of ‘Snow White’ and ‘Xiaobai’ compared to those of ‘Benihoppe’, while the expression levels of these genes in ‘Snow White’ were similar to those observed in ‘Xiaobai’. In addition, the expression levels of eight structural genes in the octoploid strawberry, namely, PAL, CHS, CHI, DFR, F3H, ANS, UFGT, and 3-GT, were lower in ‘Snow White’ than in ‘Benihoppe’ and ‘Xiaobai’. These results suggest the involvement of these genes in the anthocyanin biosynthesis pathway, where enhancing or repressing their expression possibly blocks anthocyanin biosynthesis in ‘Xiaobai’ and ‘Snow White’.

We also analysed seven DEGs in the diploid genome and 30 DEGs in the octoploid genome involved in anthocyanin or flavonoid biosynthesis pathways during the fruit development and ripening process among the three cultivars (Figs. 3, 4, Tables S4, S5). Among the assayed cultivars, a DEG annotated as phenylalanine ammonia-lyase had the homologous gene of gene 22493 (diploid genome) and PAL (nine genes in the octoploid genome chromosome 6 (three) and chromosome 7 (six)). The expression profiles showed that PAL (one gene in the diploid genome and five genes in the octoploid genome) was highly expressed in ‘Benihoppe’ and ‘Xiaobai’ compared to ‘Snow White’ during the fruit development and ripening process. The expression of four genes encoding PAL in the octoploid genome showed no major differences among the three cultivars (Fig. 3). Three genes in the diploid genome and nine genes in the octoploid genome (in chromosome 2 (six genes) and chromosome 3 (three genes)) were annotated as DFR (Fig. 3). Among them, two genes in the diploid genome and two genes in the octoploid genome were all expressed at significantly higher levels in ‘Benihoppe’ and ‘Xiaobai’ than in ‘Snow White’ during the ripening stage. The expression of the DFR gene located at chromosome 3 was extremely low compared to that located at chromosome 2 among the three cultivars in the octoploid genome. Gene 20514 (homologous gene augustus_masked-Fvb6-4-processed-gene-148.10 in the octoploid genome), which was predicted to encode a polyphenol oxidase gene involved in the pigment biosynthetic process, was upregulated during the fruit development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to that observed in ‘Snow White’ (over 50-fold higher) (Fig. 4, Tables S4, S5). Polyphenol oxidase activity is responsible for red colour and stability due to the degradation of anthocyanins, which is consistent with the anthocyanin contents observed in the three assayed cultivars (Aguilar & Hernández-Brenes, 2015). Moreover, the expression profiles of gene 2729 (homologous gene augustus_masked-Fvb1-2-processed-gene-127.8 in the octoploid genome), which is predicted to encode a UDP-glycosyltransferase 91A1 (UGT 91A1) involved in the flavonoid biosynthetic process, was upregulated during the fruit development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to that observed ‘Snow White’ (over 2-fold higher) (Fig. 4, Tables S4, S5). These results demonstrated that these genes could be involved in the anthocyanin or flavonoid biosynthetic pathways to regulate anthocyanin biosynthesis in red and white strawberry fruits by altering their transcription levels.

DEGs involved in altering the sugar contents and sugar-related pathways in the three strawberry cultivars

To assess the association between sugar and fruit colour, we measured the contents of three sugars and analysed the expression profiles of DEGs involved in sugar-related genes according to a GO enrichment analysis (Fig. 4, Tables S4, S5). The concentration of sucrose was the highest (approximately 120 mg/g) in strawberry fruits at the S7 stage, followed by fructose (approximately 10 mg/g), and a relatively low concentration of glucose (approximately 5 mg/g) (Figs. 4A–4C, Table S8). The sucrose content was significantly different between ‘Benihoppe’ and ‘Snow White’ and was higher in these two cultivars than in ‘Xiaobai’, where the sucrose content increased with fruit development and ripening (Figs. 4A–4C, Table S8). In addition, significant differences were observed in the glucose contents of ‘Benihoppe’ compared to ‘Xiaobai’ and ‘Snow White’ at the S7 stage. It is interesting to note that the fructose content in ‘Xiaobai’ was significantly higher than that observed in ‘Benihoppe’ and ‘Snow White’ at the S2, S5, and S7 stages. The sucrose and glucose contents increased with fruit ripening stages, whereas that of fructose did not (Figs. 4A–4C). These results indicated that sugar is associated with fruit ripening.

Sugars are the primary components of fruit soluble solids that govern changes in fruit quality, which depends on starch and sucrose metabolism. Six DEGs in the diploid genome and 42 DEGs in the octoploid genome were identified as being involved in the sucrose, glucose, fructose-related pathways (Fig. 4, Tables S4, S5). Three DEGs (bidirectional sugar transporter, UDP-arabinose 4-epimerase 1, and polyol transporter 5-like (PLT 5-like)) have the same annotation and exhibited similar expression profiles between the octoploid and diploid genomes. It is interesting to note that the bidirectional sugar transporter is present in both the diploid (one gene, gene26963) and octoploid (three genes, maker-Fvb2-3-augustus-gene-86.43, maker-Fvb2-4-snap-gene-183.59, maker-Fvb3-2-snap-gene-98.73) genomes, and their expressions was upregulated in ‘Benihoppe’ and ‘Xiaobai’ (approximately 10-fold higher expression levels than in ‘Snow White’) compared to ‘Snow White’ (Fig. 4, Tables S4, S5). Moreover, UDP-arabinose 4-epimerase 1 expression was upregulated during the fruit development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to that observed in ‘Snow White’ (over 5-fold higher) (Fig. 4, Tables S4, S5). In contrast, PLT 5-like expression was downregulated during the fruit development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to that observed in ‘Snow White’ (over 3-fold higher) (Fig. 4, Tables S4, S5). The expression profiles of these DEGs were similar or opposite to the anthocyanin accumulation trend observed during the fruit development and ripening process among the three cultivars, indicating that these genes may be involved in anthocyanin biosynthesis.

Genes involved in hormone biosynthesis or signal transduction among the three strawberry cultivars

To investigate hormone biosynthesis or pathways related to the formation of fruit colour, we analysed the expression profiles of DEGs involved in hormone-related genes according to a GO enrichment analysis (Fig. 4, Tables S4, S5). Twenty-three DEGs in the diploid genome and 101 DEGs in the octoploid genome were identified as being involved in hormone-related pathways. Auxin-induced protein 15a-like, protein phosphatase 2c 37 (PP2C 37), and gibberellin 2-beta-dioxygenase (GA2Ox) have the same annotation and expression profiles in the octoploid and diploid genomes. Among these factors, auxin-induced protein 15a-like was encoded by one gene (gene 4393) in the diploid genome and seven genes in the octoploid genome, which were located at chromosomes 2 and 5, respectively (Fig. 4, Tables S4, S5). The expression of auxin-induced protein 15a-like was higher in ‘Snow White’ than in ‘Benihoppe’ and ‘Xiaobai’ (one gene in the diploid genome and five genes in the octoploid genome). In addition, gene 25059, annotated as PP2C 37, was highly expressed in three varieties (insignificant difference) in the diploid genome (Fig. 4, Table S4). However, PP2C had three geneIDs and was located at chromosome 7 in the octoploid genome (Fig. 4, Table S5). In the octoploid genome, two PP2C 37 genes were more highly expressed in ‘Benihoppe’ and ‘Xiaobai’ compared to ‘Snow White’ (approximately 20-fold higher). In contrast, only one gene (PP2C 37) was more highly expressed in ‘Snow White’ than in ‘Benihoppe’ and ‘Xiaobai’ (approximately 10-fold higher). With respect to GA2Ox, one gene (gene 7769 in the diploid genome) and four genes (in the octoploid genome) were upregulated during the fruit development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to that observed in ‘Snow White’ (two genes) (Fig. 4, Tables S4, S5). In the octoploid genome, the expression profiles of another two GA2Ox genes were downregulated during the development and ripening process in ‘Benihoppe’ and ‘Xiaobai’ compared to ‘Snow White’. The expression profiles of these DEGs indicated that these genes function in hormone biosynthesis or signal pathways by controlling hormone levels to regulate anthocyanin contents.

Relationship between transcription factors and anthocyanin

To elucidate the relationship between transcription factors and the anthocyanin contents of ‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’, we extracted 413 TFs from the DEGs identified at the three stages among the three cultivars and further divided them into 25 TF families according to the diploid genome (Fig. 5A). The majority of the TFs encoding DEGs were members of the MYB family, followed by the bHLH, ERF, WRKY, and NAC families (Fig. 5A, Tables S4, S5). In the present study, 17 TFs were differentially expressed in white and red fruit in the diploid genome, including the bHLH family (3 DEGs), B3 family (3 DEGs), and MADS family (2 DEGs) (Fig. 5B, Tables S4, S5). In contrast, 104 TFs were identified in the octoploid genome among the three cultivars and were further divided them into 13 TF families (contain others TFs). Most TFs belonged to bHLH family (23 DEGs), the ERF family (15 DEGs), the MYB family (11 DEGs), and the WRKY family (18 DEGs) (Fig. 5C, Tables S4, S5). Only three DEGs (bHLH130, MADS23, and WRKY22) were identified as having the same annotations in the diploid and octoploid genomes. The expression pattern of bHLH130 and MADS23 was higher in ‘Snow White’ compared to that observed in ‘Benihoppe’ and ‘Xiaobai’ and contrasted with the pattern of anthocyanin accumulation observed in the three cultivars (Tables S4, S5). However, WRKY22 was highly expressed in ‘Benihoppe’ compared to ‘Xiaobai’ and ‘Snow White’, consistent with the accumulation of anthocyanin in strawberry fruit observed at the three stages among the three cultivars (Tables S4, S5). It is interesting to note that MYB44 was located on three chromosomes (1, 2, and 5) and MYB44 on chromosome 1 were more highly expressed than MYB44 on chromosome 2 and 5. MYB44, which is located on chromosome 1, was upregulated in red fruit compared to white fruit. MYB44 located on chromosomes 2 and 5 was downregulated in red fruit compared to white fruit (Table S5). In our present study, we identified many types of TFs that were differentially expressed during fruit development and ripening among the three cultivars, which suggested a strong association between these TFs and anthocyanin in strawberry.

Validation of RNA-seq results using qRT-PCR

To validate the RNA-seq-based DEG data, we quantified the expression of nine DEGs related to anthocyanin biosynthesis using quantitative real-time PCR (qRT-PCR) for the ‘Benihoppe’, ‘Xiaobai’, and ‘Snow White’ cultivars (Figs. 6A–6I, Table S9). All selected DEGs included structural and regulatory genes involved in the anthocyanin biosynthesis pathway, such as FaPAL1, FaF3H, FaDFR, FabHLH, and FaTT12 (F. ananassa transparent testa 12). In addition, the expression patterns of FaPAL1, FaF3H, FaDFR, FabHLH, and FaTT12 were observed to be consistent with their expression profiles in diploid and octoploid genomes. The qRT-PCR results were similar to our RNA-seq results, supporting the reliability of our RNA-seq data.