Expression Profiling Uncovers Novel Genes Associated with Plant Hormone Signal Transduction Involved in Regulation of Anthocyanin Biosynthesis in Red-Fleshed Apple Cultivars

1. Introduction

Apples are one of the most widely cultivated crops in the world temperate regions worldwide. Their cultivation is very popular, especially due to their high ecological adaptability and nutritional value of the fruits [1,2]. The global production of apples is about 78.4 million tons. Poland is the largest producer of apple fruits in the European Union (4.3 million tons), third in the world after China (35-40 million tons) and USA (over 4.5 million tons) [3]. The management of such massive fruit production raises challenges in terms of cultivation and storage. However, an opportunity to address this situation lies in increasing of consumption of fresh apples and their products. This can be achieved by introduction into apple industry red-fleshed varieties, rich in anthocyanins, having health-promoting properties [4]. In addition apples with red flesh are an interesting raw material for the production of cloudy apple juices [5].

Apple wildly spread breeding efforts allowed to confirm that Malus sieviersii is the most closely related to the commonly cultivated ‘Golden Delicious’ cultivar [6,7]. Moreover, the primary donors of genes for the red skin and red flesh of apple fruits are Malus sieversii f. niedzwetzkyana and Malus pumila var. niedzwetzkyana [7]. As red fruits are highly attractive to consumers and often associated with apple ripeness and good flavor, their pigmentation becomes an important trait in horticultural research [2,8]. However, there are only few breeding programs focused on developing new red-fleshed cultivars with high quality of the fruits [9]. An analysis of apple segregating F₁ populations combined with QTL mapping as well as transcriptomic studies, have highlighted the complexity of this trait [10,11,12,13].

The anthocyanins (glycosides of anthocyanidins), which are synthetized by flavonoid pathway, are the most important fruit component contributing in red-skinned and red-fleshed apple fruits, as well as red apple plant tissues (e.g. stems, leaves and seeds). They have been extensively studied especially in regard to human health, playing the huge role in reduction of risks diseases such as: cardiovascular, asthma, gastrointestinal, weigh regulation etc. [12,14,15]. The flavonoid pathway is well-recognized and consist the number of different enzymes catalyzes biosynthesis of anthocyanidins and anthocyanins [16]. In apple, anthocyanidins are divided into several subclasses: pelargonidin, cyanidin, delphinidin, peonidin, petunidin, malvidin, and almost 80% of them are cyanidins, thus the reconstruction of galactose into form cyanidin galactoside is mainly based on transformation of galactose to cyanidin-3-galactoside [17]. Determination of anthocyanidins coloration in specific tissues is revealed by different types and numbers of substituents in the ring of the anthocyanin [9,18,19,20,21,22].

Initial precursor in anthocyanins and flavonoids biosynthesis pathway is a phenylalanine. It is subsequently catalyzed by phenylalanine ammonialyase (PAL), 4-coumarate:coenzyme A ligase (4CL), chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), flavonoid 3’5′-hydroxylase (F3’5’H), and dihydroflavonol 4-reductase (DFR) leading to colorless metabolites production. Then, the colored anthocyanins are synthesized by anthocyanidin synthase (ANS). Finally, they are transformed into brick-red, magenta or even blue–violet glycosides through addition of UDP-glucose by flavonoid 3-O-glucosyltransferase (UFGT) [23,24,25,26,27].

Based on apple genome analysis, a number of genes encoding MYB transcription factors involved in anthocyanin biosynthesis regulation were identified [9,28,29,30]. Overexpression of MYB10 strongly correlates with apple type 1 flesh color regulation by anthocyanin increment in skin, shoots and fruit flesh (expressing the pigmentation of fruit flesh, cortex and white core). This can be a result of autoregulation of MYB10 leading to red flesh apple phenotype [19,31,32,33,34]. Meanwhile in the type 2, apple red-fleshed fruit phenotype did not co-segregate with MYB10 [9]. An analysis performed by Chagné et al. [11] and Mahmoudi et al. [26] confirmed the relationship between the MYB allele combinations and fruit flesh color. Characterization of 16 apple cultivars with red and white flesh as well as segregating population derived from the cross of ‘Geneva’ and ‘Braeburn’, confirmed that the seedlings possessing MYB1/MYB1 allele combination had red skin and white fruit flesh phenotype, while the presence of MYB10 alleles (MYB1/MYB10 or MYB10/MYB10) shaped the genotypes with red skin, red fruit flesh and red seeds [11,26].

Recent studies, established by the analysis of different set of apple genotypes as well as transgenic ‘Royal Gala’, allowed to investigate new genes such as MdJa2, MdNAC, WD40, bHLH or SEPALLATA and MADS-box transcription factors, involved in brassinosteroids accumulation and affecting anthocyanin and proantocyanidins (PA) biosynthesis. MYB-bHLH-WD40 (MBW; a composition of R2R3-MYB and WD transcription factors) was described as a part of the transcriptional complex of the anthocyanin structural genes regulation [2,13,19,35,36,37]. Simultaneously, current studies also confirmed that transcription factors such as MYB17, MYB111, MYBL2 and MYB16, interacting with HLH33 and HB1 may negatively regulate the anthocyanin biosynthesis [38,39].

Since the mechanism of regulation of red-fleshed apple fruits is complex and still unclear, the objective of this study was to investigate new, so far unrecognized, genes involved in apple fruit flesh coloration. They can become candidates to develop reliable molecular markers to accelerate an early selection of apple cultivars with red-fleshed fruits, applicable in apple breeding programs.

2. Results

2.1. Transcriptome Profiling

The major goal of this analysis was to determine the difference between gene transcription in fruits of red- and white-fleshed apple cultivars. RNA-seq experiment with fruits of Red Love^® ‘General’ and ’Early Fuji’, characterized by red and white fruit flesh respectively, has been performed. For each cultivar 3 independent fruit samples (3 replications) were used, six RNA-seq libraries were constructed and over 40G of clean reads were obtained. The effective data volume for each sample was in the range from 6.06 to 7.23G. The quality of reads was high, with Q30 ratio from 92.33 to 92.78% of inferred base call accuracy. The GC content was similar in all samples with an average of 47.46% (Supplementary Table S1). Clean reads were mapped on the apple reference genome at a similar ratio, on average 93.63%, with 89.7% of uniquely mapping reads and 85.91% of reads mapped in proper pairing.

The expression of protein-coding genes was calculated using FPKM method, which indicates the number of fragments per kilobase length of a protein-coding gene per million of sequenced fragments. For 35 890 genes annotated in apple genome, similar number in each analyzed sample, on average 25 666 genes, were found to be expressed (Supplementary Table S2). Sample-to-sample cluster analysis and principal component analysis (PCA) of expression profiles indicated that the results were consistent and investigated cultivars samples clearly are grouped together (Figure 1).

2.2. Differentially Expressed Genes

Based on gene expression profiling, differentially expressed genes (DEGs) in Red Love® ‘General’ as compared to ‘Early Fuji’ control were identified (log2 fold-change >1 false discovery rate FDR < 0.05) (Supplementary Table S3). In total 4 839 genes were found to be differentially expressed in Red Love® ‘General’ with 2 259 genes up-regulated and 2 580 down-regulated (Figure 2).

2.3. Functional Categories of DEGs

2.3.1. GO Enrichment Analysis of DEGs

After finding the differentially expressed genes in red-fleshed Red Love® ‘General’, we performed Gene Ontology (GO) enrichment analysis to describe their functions (Figure 3).

The GO enrichment analysis allowed classifying the general role of differentially expressed genes in: biological processes (BP), cellular component integration (CC) and molecular function (MF). The total GO enrichment categorization of 3 500 DEGs (Supplementary Table S4 and Supplementary Figure S1) led to assigned 1 635 genes as up-regulated and 1 865 as down-regulated in red-fleshed fruits of Red Love® ‘General’ as compared to white-fleshed ’Early Fuji’ (Figure 3 and Supplementary table S5). In total most of DEGs were clustered into 65 different GO terms: 23 for biological process, 20 for cell component and 21 for molecular function (Figure 3). In terms of biological process the highest number of DEGs was assigned to categories: biological regulation, cellular process, metabolic process, response to stimulus and single-organism process. In terms of cellular component, the highest number of DEGs was assigned to the categories: cell and cell part, membrane and membrane part integration and organelle and organelle part. In terms of molecular function, most of genes were assigned into two categories: binding and catalytic activity (Figure 3).

Next, DEGs from top 30 GO enriched categories were analyzed (Supplementary Figure S2). In the main category of biological process, top up-regulated DEGs (Supplementary Figure S3) were classified in terms to: GO:0042445 – hormone metabolism process, GO:0008202 – steroid metabolic process, GO:0009850 – auxin metabolic process and GO:0009684 – indoleacetic acid biosynthesis. There were other terms additionally enriched: GO:0023052 – signaling, GO:0050896 – response to stimulus, GO:0044700 – single organism signaling, GO:0006952 – defense response, GO:0007165 – signal transduction, GO:0007154 – cell communication, GO:0034050 – host programmed cell death, and GO:0009626 – plant hypo-sensitivity. In the category of cell component genes, two top GO enrichment terms were extracted: GO:0071944 - cell periphery and GO:0005886 – plasma membrane (cell membrane integration) (Supplementary Figure S4). In the category of molecular function, top GO enrichment terms were assigned to: GO:0043531 – ADP binding and GO:0016491 - oxidoreductase activity (Supplementary Figure S5).

2.3.2. KEGG Enrichment Analysis of DEGs

Analysis of significantly enriched pathways according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed that the highest numbers of DEGs were classified in six main groups: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, organismal systems (Figure 4a). The biggest group included genes involved in metabolic pathway regulation (in total 812 DEGs; blue blocks) and second - in genetic information processing (in total 205 DEGs; green blocks). In addition, the third group of genes was classified in pathways of plant metabolism in diverse environments (in total 99 DEGs; brown blocks). The subgroups of DEGs, classified in metabolism were assigned to: carbohydrate metabolism (170), amino acid metabolism (110), energy metabolism (105), and lipid metabolism (99). From the other groups the highest number of DEGs (90), was classified to the signal transduction (Figure 4b). Extracted up and down regulated genes recovered from Red Love® ‘General’ vs. ‘Early Fuji’ compartments are presented in Supplementary Figure S6. From this grouping, 30 top GO terms uncovered, three most abundant one from biological process were assigned to the galacturonate biosynthetic process and signal transduction; from cellular components – plasma membrane and extracellular space and form the group of molecular functions – genes involved in ADP binding and representing chalcone isomerase activity (Figure 4b)

Detailed analysis of the KEGG pathways enrichment classification showed that the group of DEGs most significantly up-regulated in red-fleshed Red Love® ‘General’ was assigned to plant hormone signal transduction pathway (mdm04075) (Figure 5). In the next, ten up-regulated DEGs from this group, involved in auxin, cytokinine, gibberellin, abscisic acid, brassinosteroids, jasmonic acid and salicylic acid hormonal signaling and regulation, were selected for further study (Figure 6).

The functions and cell localization of selected genes are presented in Table 3 and the selected genes were addressed for detailed expression profiling in the apple cultivars characterized by red-fleshed fruits.

2.4. Validation of Activity of Genes from Plant Hormone and Signal Transduction Pathway by RT-qPCR

The expression profiles of selected genes were analyzed in immature and ripe fruits of seven red-fleshed apple cultivars in comparsion to the immature and ripe white-fleshed fruits of ‘Free Redstar’ (Figure 7).

Overall, in comparison to ‘Free Redstar’, significant down-regulation of five, out of ten selected genes: DELLA (involved in gibberelin pathway), ABF (involved in abscisic acid pathway), BZR1 and TCH4 (brassinosteroids pathway), and MYC2 (responsible fo jasmonic acid transduction) was observed in immature fruits of: ‘Trinity’, ‘Alex Red’, Red Love® ‘Era’, Red Love® ‘Circe’, ‘Roxana’ and ‘Red Love® ‘Sirena’ (Figure 8). Interestly the most significant up regulation of those genes was evaluated for M. sieviersii f. niedzwetzkyana (characterized for red coloration of leaves, flower petals, core, stem, seeds, skin and flesh), representing wild apple variety, and concidered as the donor of target genes. Simultinously relevant overexpression of the same set of genes was noted for control immature white fruit samples of ‘Free Redstar’,

In case of ripe fruits of red fleshed cultivars all genes applied in this study showed significant variation in expression profiles. The activities of GH3, SnRK2, BRI1 and TCH4 genes were significanlty higher in fruit samples of ‘Alex Red’, Red Love® ‘Circe’ and Red Love® ‘Sirena’. In case of ‘Trinity’, Red Love® ‘Era’ and ‘Roxana’ cultivars four genes: B-ARR (cytokinine pathway), DELLA (gibberellin pathway), BZR1 (brassinosteroids pathway), and MYC2 (jasmonic acid pathway) showed significant down regulation. High expression level of ABF gene (abscisic acid pathway) (ranged between 10x to 50x number of transcript fold change) in comparsion to ‘Free Redsar’ white ripe fruit control, was noted for all red fleshed apple cultivars evaluated. Overall, the activity of analyzed genes mantained overexpressed in immature and ripe fruits of M. sieviersii f. niedzwetzkyana, meanwhile inactivation of all selected genes was observed in ripe fruits of white-fleshed ‘Free Redstar’ (Figure 9).

2.5. The activity of the Structural Genes in Red-Fleshed Apple Cultivars

Six genes, from polyphenol and anthocyanin biosynthesis pathway were analyzed to confirm their expression in red fleshed apple cultivars. Significantly high level of expression of all structural genes was noted in immature fruits of M. sieviersii f. niedzwetzkyana (i.e. showing about 18x fold changes for F3H) and ‘Free Redstar’ white-fleshed apple fruits (Figure 10). In case of DFR and UFGT negligible activity was noted for analyzed immature fruit flesh samples (data not shown).

Additionally, CHS, F3H and UFGT genes were significantly up-regulated in ripe fruits of ‘Trinity’, Red Love® ‘Circe’ and Red Love® ‘Sirena’. In general, significant up-regulation (fold change ranged from 10 to 60x) of PAL, CHS, ANS and F3H (encoding enzymes for anthocyanin precursors synthesis) was observed for all ripe fruits of red-fleshed apples. As we could expect, the activity of tested genes was significantly lower in ripe fruits of white-fleshed ‘Free Redstar’ (Figure 11). Simultaneously, the light red fruit flesh color of ’Roxana’ cv. may be the result of significantly higher expression of the PAL (60 fold change) precursor gene and relative inhibition of other structural genes, probably causing the breakdown of intensity of anthocyanin biosynthesis in the genome of this cultivar. Interestly, probably similar mechanism of this regulation was noted for white-fleshed ’Free Redstar’ (Figure 11). Generally, in red-fleshed apple cultivars four of structural genes (PAL, CHS, F3H, ANS), were significantly up-regulated, while high activity of DFR and UFGT genes was noted only for ‘Alex Red’, Red love® ‘Circe’ and ‘Red Love® ‘Sirena’.

Additionally, the lower activity of all tested genes was observed for ripe, in comparison to immature, fruits of M. sieversii f. niedzwetzkyana.

2.6. Genes from Signal Transduction Pathway Shows Inconsistencies in Regard Expression Profiles in Comparison to RNA-seq Experiment but Significant Intergenic Correlation

Genes from the plant signal transduction pathway selected on the basis of the RNA-seq experiment were expressed at high level in the fruits of the Red Love® ’General’ cultivated in China. They were chosen for further expression profiling with regard to transcriptome analysis. However, in the validation studies performed on the red-fleshed apple collection, grown in Poland, variable and inconsistent expression profiles of those genes were observed.

Gene-to-gene correlation coefficient between structural and uncovered genes was evaluated individually for a set of immature and ripe fruits of analyzed apple cultivars. The assessment of interaction between structural genes and those selected from the hormone signal transduction pathway in immature fruits of all tested cultivars, showed significant correlation between, i.e.: DFR F3H, UFGT and B-ARR, NPR1, GH3, SnRK2, BRI1 (Table 4A). In ripe fruits, the correlation was noted between the all identified genes and the structural UFGT gene (encoding flavonoid 3-O-glucosyl transferase, the gene responsible for final glycosides transformation and their pigmentation). This may explain that the most of the structural genes thus from anthocyanin pathway as well as plant hormone signal transduction pathway are potentially activated at early fruit developmental stage and shows significant intergenic correlation.

Generally, in ripe fruits significant correlation was observed between chalcone synthase gene (CHS) and B-AAR, GH3, SnRK2 and TCH4 (respectively involved in pathways of: cytokinins, auxin, salicylic acid and brassinosteroids), between gene encoding phenylalanine ammonialyase (PAL) and the same set of genes, as well as between UFGT and all genes selected form hormone transduction pathway (Table 4B.). The B-ARR, GH3, SnRK2 and TCH4 genes (from cytokinine, auxin, abscisic acid and brassinosteroids pathways respectively) showed correlation with CHS, PAL and UFGT genes and could be considered as possible indicators of regulation of glucoside color transformation in apple red fleshed fruits ready for harvesting.

3. Discussion

In the presented study, utilizing transcriptome comparisons and gene expression analysis of two apple cultivars producing red- and white-fleshed fruits, we have revealed new genes putatively related with the regulation of anthocyanin biosynthesis. Among of over 40 Gbase of filtered sequence data, two major groups of genes mapped on the reference apple genome were recovered. One group of sequence reads (43%) were assigned to be involved in metabolic pathways and second (31%) - in secondary metabolite biosynthesis. Based on differentially expressed genes (DEGs) form top 20 KEGG Gene Ontology (GO) term classifications (Figure 5; Results), we selected a set of 10 genes from plant hormone signal transduction pathway, involved in auxin (AUX), cytokinine (CK), gibberellin (GA), abscisic acid (ABA), brassinosteroids (BR), jasmonic acid (JA) and salicylic acid (SA) biosynthesis. Following RT-qPCR analysis of the genes of interests we have observed their variable activity with up-regulation noted for GH3, SnRK2, ABF and TCH4, and down-regulation for B-ARR, DELLA, BZR1, BRI1, MYC2 and NPR1. However, although all the selected genes showed up-regulation in red-fleshed Red Love^® ‘General’ cultivar used in RNA-seq experiment, the results obtained for the RT-qPCR, performed for seven cultivars studied, were not consistently uniform. A similar observation of inconsistencies between RNA-seq and RT-qPCR comprehensive analysis was reported by Everaert et al. [40]. The authors concluded that, depending on the analysis workflow, 15–20% of genes usualy are considered as ‘non-concordant’ in regard to the results obtained with RNA-seq and RT-qPCR. Those ‘non-concordant’ genes are defined when both approaches yield differential expression in opposing directions or when one method shows differential expression while the other does not [40].

The analysis carried out for representative set of seven red-fleshed apple cultivars collected at two fruit developmental stages, and compared to white-fleshed cv. ’Free Redstar’, allowed to show, that the changes in gene expression are closely related to the fruit ripening process [34,41,42]. All structural genes and the newly revealed GH3 (auxin response), SnRK2 and ABF (both involved in ABA signaling) as well as TCH4 (BR pathway) were activated in the ripe red-fleshed fruits. In accordance to fold change calculation, DELLA (GA response), B-ARR (CK response), BZR1 and BRI1 (BR receptors), as well as MYC2 and NPR1 (genes from JA and SA signaling pathways respectively), were defined as to be significantly inhibited in the ripe red-fleshed fruits. The putative mechanism of gene regulation of the plant hormone and transduction pathway is presented on the scheme below (Figure 12).

In accordance to presented scheme, similar observation was described by Su and coworkers [37]. They postulated that the structural genes from anthocyanin biosynthesis pathway were up-regulated in apple calli, suppressing expression of brassinosteroids resistant gene (BZR1) (BR accumulation = anthocyanin pathway activation). They also have reported that interaction of BZR1 with MdJa2 modulates the target genes transcription, and generally regulate plant response to brassinosteroids [37]. Based on observation of apple calli development, Zheng and coworkers also underlined the anthocyanin biosynthesis regulation by plant hormones, including brassinosteroids, which strictly influenced their accumulation [43]. Simultaneously, as it was reported by Su et al., application of external BR also promoted fruit maturation and pro-anthocyanin synthesis in apple tissue [37,43]. Our observations led to the comparable conclusion that the inactivation of BZR1 is positively correlated with promotion of anthocyanin accumulation in red-fleshed ripe apple fruits.

Brassinosteroids signaling is perceived by receptors BRI1 and BZR1 in cell membrane. As the recent research explains it is one of the major factors (similarly to bHLH) playing key role in regulation of brassinosteroids gene expression [43,44. Both, brassinosteroid receptor genes, analyzed in presented study, showed significantly low expression in the series of red-fleshed fruit cultivars, probably leading to brassinosteroids accumulation, thus promoting an anthocyanin synthesis. This selecting mechanism was described in previous studies of tomato, cucumber, strawberry and grape, indicates the key role of brassinosteroids in fleshy fruits development and ripening [42,45,46,47,48].

In parallel, based on our observations, we have found the significant inactivation of genes NPR1 and MYC2, respectively responsible for salicylic acid (SA) and jasmonic acid (JA) signal transduction, influencing in negative correlation in anthocyanin regulation in apple fruits. Therefore, our observation gives new insight in anthocyanin biosynthesis dependence on salicylic and jasmonic acids signal inhibition, not reported so far.

In this concern, previous research of Sasaki et al. and summarized by Wang and Chang explained the role of JA in regulation of secondary metabolites, plant defense response and organ development [49,2]. Then, later studies of Qi, An, and coworkers showed, that the activation of MYB and bHLH, significantly correlated with anthocyanin accumulation in A. thaliana and it is promoted by complex mechanism, including degradation of repressor protein (so called JA ZIM DOMAIN) [50,51]. Recent reports of Wang et al. also underline the role of jasmonic acid as signaling factor influencing the MYB24L gene and affecting the anthocyanin biosynthesis [30].

In conducted research, the newly studied DEGs such as DELLA (gibberellin pathway) and B-ARR (cytokinin pathway), showed significant down-regulation in the red-fleshed ripe fruits and up-regulation in white fruit flesh of ‘Free Redstar’. The observations are similar to the results obtained by Nawaz and coworkers, which confirmed up regulation of histidine kinase and DELLA genes in ripe fruits of white-fleshed apple ’Hanfu’ and its mutant. The authors explained, the level of gene transcript number is genotype-depended, resulting in apple flesh fruit pigmentation intensity; from flush flesh, pinky flesh, to dark red flesh [42,52].

Another founding of our work underlines negative correlation between gene activity and auxin and cytokinin production (gene activation = less auxins or gene inactivation = more cytokinin production) expressed by up-regulation of auxin response gene (GH3) and down regulation of (B-AAR) in red-fleshed cultivars (Figure 10, Results). This gene activity and metabolite production relationship was highlighted by Ji, Dikeman and coworkers who observed that increment of auxins together with decreasing of cytokinin in the cells can inhibit anthocyanin synthesis [53,54]. This phenomenon was observed in callus tissue of both A. thaliana and apple and underlined by Nguyen, Ji and coworkers as they have elucidated that the increased concentration of cytokinin and decreasing of auxin can significantly promote the expression of MYB transcription factors, leading to anthocyanin accumulation [55,56].

Additionally, we have noted up-regulation of ABF and SnRK2 genes, (involved in abscisic acid pathway) in red-fleshed ripe fruits. Notably, this was observed in all cultivars except ‘Roxana’ (pink fruit flesh) and M. sieversii f. niedzwetzkyana (a wild apple variety). This observation suggests a distinct regulation of fruit flesh coloration in these two cultivars. The gene activity in this context is genotype-dependent, highlighting the complexity of regulation of apple fruit flesh coloration. As it was described by Jiang and Joyce ABA plays the kay role in strawberry fruit color development (generally through PAL activity stimulation), as well as in grapes and litchi promotes the fruit coloration stage [57,58,59].

Through the effect of regulating the number of gene transcripts, we have noted the changes in targeted hormonal pathways, which, by regulating the content of hormonal proteins, influence: the mechanism of anthocyanins synthesis and accumulation in apple fruits. For the DELLA, B-ARR, BZR1, BRI1, MYC and NPR1 genes, we have observed a negative regulatory mechanism, by blocking the gene activity and thus resulted in hampering recognition or degradation of hormonal proteins.

The uncovered genes were correlated with the functional gene UFGT (responsible for final anthocyanin color transformation). To date, only limited studies have examined the negative regulatory mechanisms underlying anthocyanin synthesis in apple. The novel technologies applied in our research allowed understanding the mechanism of apple fruit color development. Uncovered genes being the candidates for molecular markers can accelerate the breeding of high quality red-fleshed apple cultivars. Application of such research provides the knowledge of regulation of fruit color development and anthocyanin biosynthesis, accumulation and degradation.

4. Materials and Methods

4.1. Plant Material

In RNA-seq experiment red fleshed Red Love® ‘General’ and white fleshed ‘Early Fuji’ apple cultivars were used. Apple fruits were collected from trees grown in experimental orchard of the Yantai Academy of Apple Research in China. Fruits (min. 3) of both cultivars were picked up at the harvest time (middle of September). From each fruit peeled flesh discs with the diameter of 2 cm and 1 cm of depth, cut with sterile blade, were collected. The samples were placed directly into liquid nitrogen and stored in -80℃ until RNA extraction.

In validation of gene expression experiments a set of apple cultivars characterized by different patterns of fruit pigmentation was used. This set included: ‘Alex Red’, ‘Roxana’, ‘Trinity’, Red Love® ‘Circe’ (‘Lureone’), Red Love® ‘Era’, Red love® ‘Sirena’ and Malus sieviersii f. niedzwetzkyana and a control white fruit flesh ‘Free Redstar’ (flesh pigmentation control). Tissue samples were collected, from the apple trees grown in the experimental orchard maintained at the National Institute of Horticultural Research in Skierniewice, Poland. For each cultivar three fruits were collected at two developmental stages: immature fruits sampled 100 days after full bloom (DAFB) and ripe fruits collected at harvest time. Then, the apple fruit flesh discs (as describe above) were prepared and immediately placed in liquid nitrogen for the RNA isolation. The samples were used for Reverse Transcriptase quantitative PCR (RT-qPCR) gene expression analysis of target structural genes as well as for differentially expressed genes revealed by RNA-seq.

4.2. RNA Extraction

For RNA-seq experiment, total RNA was isolated using MINIBEST Plant RNA Extraction Kit (TaKaRa, China) according to manufacture protocol. The concentration and quality of RNA was evaluated after electrophoresis in 1% agarose gel stained with GelRed Dye (Biotiun, China). High quality RNA preparations were sent to the commercial company for transcriptome sequencing.

For RT-qPCR total RNA was isolated using method described by Zeng and Yang [60] with minor modifications. The procedure was started with sample incubation for 20 minutes in CTAB buffer heated to 65℃. After centrifugation in a mixture of chloroform and isoamyl alcohol (24:1, v/v) RNA was resolved from the upper solution fraction and precipitated with 10 M lithium chloride overnight at 4°C. Then, the precipitate was dissolved in RNase-free water. The quality, degree of integration and concentration of RNA samples were assessed by using the Agilent 2100 Bioanalyzer and Expert 2100 software (Agilent Technologies, CA, USA).

4.3. RNA-seq Analysis

The library construction and sequencing were completed by Shanghai Ouyi Biomedical Technology Co., Ltd. The transcriptome sequencing by the Illumina technology was conducted by OE biotech Co., Ltd. (Shanghai, China). To perform quality check of the raw data and summarize the number of reads Trimmomatic software [61] was used. To map sequence reads on apple reference genome Htseq-count software [62] was used. As a reference ‘Golden Delicious’ apple genome assembly ASM211411v1 was applied: https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_002114115.1/. The DESeq [63] allowed standardize the number of counts of each single gene, calculate the multiple of difference and perform significant difference tests. To obtain the number of reads compared to protein-coding genes in each sample cufflinks [64] was applied. Finally, the gene expression level was calculated as FPKM (Fragments Per kb Per Million Reads), giving the number of fragments per kilobase length of a protein-coding gene per million fragments sequenced. The raw sequence data were stored in FASTQ file format. Evaluation of differentially expressed genes (DEGs) between Red Love® ‘General’ and ‘Early Fuji’ was performed with use of DEseq2 R package software. For this purpose RNA-seq data were compared and analyzed weather the same gene is differentially expressed within two samples. For this purpose, LogFC >1 (logarithm of fold change in the expression level of the same gene in two analyzed samples) and filtering data condition p-value<0.05 criteria were applied.

4.4. GO Enrichment Analysis

The gene function enrichment was performed based on compilation of the list of all protein coding genes and differential protein-coding genes (counted based on FPKM), using hypergeometric distribution test to calculate the representative GO function set. Than p-value of significant enrichment in the list of differential protein-coding genes was corrected by Benjamini & Hochberg’s multiple test to obtain fails discovery rate (FDR). The procedure translates experimental results into a series of regulated gene lists at multiple false discovery rate FDR cutoffs, and computes the p-value of the overrepresentation of the genes.

The top30 GO enrichments analysis was performed based on screened GO entries corresponding to the number of differential genes greater than 2 in the three categories (Biological Process (BP), Cell Components (CC), and Molecular Function (MF)), and sorted 10 entries according to the -log10Pvalue corresponding to each entry from large to small layout.

The fisher algorithm was used to conduct Cellular Components (CC), Biological Process (BP), and Molecular Function (MF) enrichment analysis on the differential genes between samples [65], and topGO was used to draw a directed acyclic graph for the enriched term to visually display the GO term of differentially expressed gene enrichment and their hierarchical (more specific) relationship.

4.5. KEGG Enrichment Analysis of Differential Genes

KEGG (Kyoto Encyclopedia of Genes and Genomes) database was screened for pathway analysis on differential protein coding genes. The hypergeometric distribution test was used to calculate the significance p-value of differential gene enrichment in each pathway entry. KEGG enrichment analysis of top20 pathway entries corresponding to the number of differential genes greater than 2, were sorted by the -log10Pvalue and presented by bubble chart (Figure 5; Results).

4.6. cDNA Synthesis and RT-qPCR

Total RNA (1µg) was reverse-transcribed into cDNA using the AffinityScript QPCR cDNASynthesis Kit (Agilent, CA, USA). The reverse transcription reaction was carried out with universal oligo-dT primer and reverse transcriptase (RT) in optimized thermal conditions: 25°C for 5 min., 42°C for 5 min. (oligo-dT annealing), 55°C for 15 min. (reverse transcription), 95°C for 5 min. (enzyme inactivation) in Biometra Basic thermocycler (Biometra, Germany). For differentially expressed genes uncovered in this study oligonucleotides were designed using Primer3plus software (https://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi). For structural genes, encoding key enzymes of anthocyanin biosynthesis pathway, oligonucleotides published by Kondo et al. [24] were used. Gene encoding ACTIN [39] was used as a RT-qPCR data normalizer (Table 5).

RT-qPCR was performed with Kapa SYBR qPCR kit in the presence of SYBR Green fluorescent dye (KapaBiosystems, USA) using RotorGen 6000 thermal cycler (Corbett Research, Australia). In each experimental setup, two pairs of specific primers, complementary to the evaluated structural gene and the selected target DEGs, were used in analogous reactions. The cDNA template was prepared in dilutions of known concentrations, enabling the preparation of a standard amplification reaction curve. The thermal profile of the RT-qPCR reaction was as follow: 95°C for 5 min. (polymerase activation), than 40 cycles including steps: 95°C for 15 s (denaturation), 60°C for 20 s (oligonucleotide annealing), 72°C for 20 s (fluorescence level detection). Relative expression (fold change) normalized in regard to ACTIN was determined on the basis of single data points derived from real-time PCR amplification curve threshold cycles (Ct values) (2^-∆∆Ct) described by Livak and Schmittgen [66]. For this purpose Rotor-Gene 6000 Series Software 1.7 was used (Corbett Research, Australia). The average of relative expression was normalized to the white flesh fruit control ‘Free Redstar’. Standard error of the mean ± SEM and t-significance at the levels of p-value> 0.05 (*), 0.01 (**), 0.001 (***) between ‘Free Redstar’ and red fleshed cvs were calculated separately (GraphPad Prism 10.0.3). Relative fold changes diagrams for each gene were drawn using GraphPad Prism 10.0.3. The gene-to-gen correlation matrix between structural and newly uncovered gens was calculated using Pearson correlation coefficient (r=1, -1) with p-value<0,05 (GraphPad Prism 10.0.3 software).

5. Conclusions

In this research we have carried out the evaluation of the expression profiles of genes related to the plant hormonal pathway in apple fruits. We have defined the set of genes showing diverse activity in immature and ripe fruits of red fleshed apple cultivars. The group of selected genes, ABF, SnRK2, as well as MYC2, NPR1 and DELLA seem to be involved in positive or negative regulation of plant signaling pathways and may play crucial role, mediating the abscisic acid, jasmonic acid, salicylic acid and gibberellin pathways respectively, thus influencing in anthocyanin transformation, accumulation and biosynthesis. All investigated genes showed significant correlation with the UFGT gene in ripe red-fleshed fruits. So far, there are many reports devoted to research on the involvement of auxins, cytokinine and brassinosteroids in apple fruit ripening and fruit flesh pigmentation. The molecular basis for the regulation of other hormonal transduction factors and their impact on apple red fruit flesh phenotype constitution have not been deeply investigated so far. Despite the preliminary research conducted, further analyzes, including important fruit quality parameters, such as fruits vitamin C, polyphenols, and anthocyanins value, are necessary to determine the suitability of uncovered genes to be used for early selection of red-fleshed apple varieties by molecular marker assisted selection (MAS) procedure.