In this study, Malus spp. ‘Flame’, a green-fruited cultivar, was used. Eight-year-old trees were grafted onto Malus hupehensis and planted at the Crabapple Germplasm Resources Nursery at the Beijing University of Agriculture (40.l°N, 116.6°E). Three trees showing similar growth were used, and fruit samples were collected from annual branches growing at the edge in the southeast direction. The fruits were collected 20, 40, 60, 80, and 100 days after budding (S1-S5, ‘S’ represents ‘stage’). All flesh samples were frozen in liquid nitrogen upon collection and stored at −80°C prior to high-pressure liquid chromatography (HPLC) analysis or RNA extraction.

‘Flame’ tissue culture plants were harvested from one-year-old branches before spring bud germination, and the culturing conditions were as previously described (Tian et al., 2017).

RNA degradation and contamination were visualized on 1% agarose gels. RNA purity was confirmed using a Nano Photometer^® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit^® RNA Assay Kit in a Qubit^® 2.0 fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using an RNA Nano 6000 Assay Kit and a Bioanalyzer 2100 system (Agilent Technologies, CA, USA) (Yang et al., 2018).

A total of 3 µg of RNA per sample was used as input material for the RNA sample preparation. Sequencing libraries were generated using an NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s recommendations, and index codes were added to label each sample. To preferentially select the 150 to 200 bp cDNA fragments, an AMPure XP system was used to purify the library fragments. High-fidelity DNA polymerase, Universal PCR primers and index (X) primers were used in the PCRs. An Agilent Bioanalyzer 2100 system was used to assess the quality of the library.

An index of the reference genome was built using Bowtie v2.2.3, and paired-end clean reads were aligned to the apple (Malus domestica) reference genome using TopHat v2.0.12 (Trapnell et al., 2009; Riccardo et al., 2010). HTSeq v0.6.1 (https://pypi.python.org/pypi/HTSeq) was used to count the read numbers mapped to each gene (Anders et al., 2015).

Differential gene expression analysis of the five groups (three biological replicates per group) was performed using the DESeq R software package (1.18.0) (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) (Benjamini and Hochberg, 1995). The resulting P-values were adjusted using the Benjamini and Hochberg approach for determining the false discovery rate (Benjamini and Hochberg, 1995). Genes with an adjusted P-values < 0.05 found by DESeq were considered differentially expressed (Anders and Huber, 2010).

Blast2GO software was used to identify enriched GO terms. GO terms with corrected P < 0.05 were considered significantly enriched for DEGs (Conesa et al., 2005). KOBAS software was used to test the statistical enrichment of DEGs in the KEGG pathways (http://www.genome.jp/kegg/) (Mao et al., 2005).

The R WGCNA package was used to identify modules of highly correlated genes based on fragments per kilobase of transcript per million mapped reads (FPKM) data (Zhang and Horvath, 2005). The WGCNA analysis was performed according to established methods (Zhan et al., 2015). Genes with the highest degree of connectivity within a module are referred to as intramodular hub genes (Langfelder and Horvath, 2008). The gene annotation information was taken from the KOBAS 2.0 annotation results.

HPLC and qRT-PCR analyses were performed according to previously published methods (Tian et al., 2011). Frozen samples (approximately 0.8–1.0 g fresh weight) were extracted with 10 mL extraction solution (methanol: water: formic acid: trifluoroacetic acid = 70: 27: 2: 1) at 4°C in the dark for 72 h. The supernatant was passed through filter paper and then through a 0.22-μm Millipore™ filter (Billerica, MA, USA). Trifluoroacetic acid: formic acid: water (0.1: 2: 97.9) was used as mobile phase A, and trifluoroacetic acid: formic acid: acetonitrile: water (0.1: 2: 48: 49.9) was used as mobile phase B for HPLC analysis. The gradients used were as follows: 0 min, 30% B; 10 min, 40% B; 50 min, 55% B; 70 min, 60% B; 30 min, 80% B. Detection was performed at 520 nm for anthocyanins and at 280 nm for PAs (Revilla and Ryan, 2000). All samples analyzed consisted of three biological triplicates.

The expression levels of related genes were analyzed using qRT-PCR and SYBR Green qPCR Mix (TaKaRa, Ohtsu, Japan) with a Bio-Rad CFX96 Real-Time PCR system (BIO-RAD, USA) according to the manufacturers’ instructions. The PCR primers were designed using NCBI Primer BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and are listed in Table S1 . qRT-PCR analysis was carried out in a total volume of 20 μl containing 9 μl of 2×SYBR Green qPCR Mix (TaKaRa, Ohtsu, Japan), specific primers at 0.1 μM each, and 100 ng of template cDNA. The reaction mixtures were heated to 95°C for 30 s, followed by 39 cycles at 95°C for 10 s, 50 to 59°C for 15 s, and 72°C for 30 s. A melting curve was generated for each sample at the end of each run to ensure the purity of the amplified products. The transcript levels were normalized using the Malus 18S ribosomal RNA gene (GenBank ID DQ341382, for crabapple) as the internal control and calculated using the 2^ ^(−ΔΔCt) analysis method. All samples analyzed consisted of three biological replicates extracted from three different batches of fruits.

PA accumulation in crabapple leaves was visualized via infiltration with DMACA stain. Samples were stained via incubation overnight with DMACA solution (0.2% DMACA w/v in methanol: 6 M HCL, v/v = 1:1). The PA content was determined as previously described (Wang et al., 2017).

Full-length RAP2-4 (MD15G1365500) and RAV1 (MD13G1046100) constructs were PCR-amplified from a cDNA library derived from Malus crabapple leaves (cv. “Flame”) using gene-specific primers and Taq DNA polymerase (TaKaRa, Ohtsu, Japan) according to the manufacturer’s instructions. Full-length RAP2-4 and RAV1 were cloned into a modified pBI101 vector using seamless cloning at the NdeI and KpnI sites. The PCR primers used are shown in Table S1 (Tian et al., 2016).

A. tumefaciens cells were grown, collected, and resuspended to a final optical density of 1.5 at 600 nm in a solution of 10 mM MES, 10 mM MgCl₂, and 200 mM acetosyringone and then incubated at room temperature for 3 to 4 h without shaking. The infiltration protocol and culture methods for transient expression assays in crabapple plantlets and apple fruits were adapted as previously described (Tian et al., 2016). All samples were analyzed from at least three biological replicates.

The open reading frames of RAP2-4 and RAV1 were cloned into the EcoRI and SacI sites of pGADT7 (Clontech, Palo Alto, CA, USA) under the control of the galactokinase 4 (GAL4) promoter to yield the effector constructs. The promoter fragments of McCHS (MD13G1285100), McCHI (MD01G1117800), McF3H (MD02G1132200), McDFR (MD03G1214100), McANS (MD01G1153600), McUFGT (MD09G1141700), McLAR1 (MD16G1048500), McLAR2 (MD13G1046900), McANR1 (MD10G1311100), and McANR2 (MD05G1335600) were ligated into the pHIS2 plasmid, the sites are located upstream of the LacZ reporter gene (BD Biosciences, Shanghai, China). The background of the pHIS2 vectors was suppressed using 3-amino-1,2,4-triazole (3-AT). The yeast one-hybrid assay methods were as previously described (Wang N. et al., 2018). The primers used for the yeast one-hybrid assays are shown in Table S1 .

Raw sequencing data in this manuscript have been deposited in National Center for Biotechnology Information Sequence Read Archive under accession number PRJNA546094.

To compare the variation in flavonoid content during crabapple fruit development, we selected fruit flesh for high-performance liquid chromatography (HPLC) from three biological replicates of ‘Flame’ fruit at five different developmental stages (35, 60, 95, 120, and 150 days after full bloom) ( Figure 1A ). The main PA compounds procyanidin B1, procyanidin B2, and epicatechin were detected by HPLC, and the overall PA content (procyanidin B1, procyanidin B2, epicatechin) decreased during fruit development ( Figure 1B ). Moreover, anthocyanin accumulation was only detected at stage 5, and phloridzin accumulated in the early stages of fruit development.

We used RNA-seq analysis to profile the transcriptomes of ‘Flame’ fruit at five representative developmental stages ( Figure 1A ) (three biological replicates). The total numbers of clean reads in the RNA-seq libraries ranged from 10,664,693 to 13,963,577, and > 83% of paired reads were mapped to the apple genome ( Table 1 ). Pearson correlation analysis showed that the three biological replicates had highly consistent transcriptome profiles across all developmental stages (r ² = 0.858 to 0.986; Figure S1 ). The percentages of exonic sequences ranged from 42.95% to 56.98%, and the percentages of intronic sequences ranged from 2.19% to 5.42% ( Figure S2 ).

We also observed that many flavonoid biosynthetic genes were highly expressed during development ( Figure S3 ). Flavonoid pathway biosynthetic genes, including CHS (chalcone synthase) (MD04G1003300, MD13G1285100), CHI (chalcone isomerase) (MD01G1167300, MD07G1186300, MD07G1233400), flavanone 3 beta-hydroxylase (F3H, MD02G1132200, MD15G1246200), F3H (flavanone 3-hydrocylase) (MD02G1132200), dihydroflavonol 4-reductase/flavanone 4-reductase (DFR, MD15G1024100), flavonoid 3’-monooxygenase (F3’H, MD14G1210700, MD06G1201700), leucoanthocyanidin dioxygenase (LDOX, MD06G1071600, MD03G1001100), and flavonol synthase (FLS, MD08G1121600, MD15G1353800) showed > 11-fold differential expression ( Table S3 ). Notably, the PA-related biosynthetic genes LAR (MD13G1046900, MD16G1048500) and ANR (MD05G1335600, MD10G1311100), encoding enzymes associated with PA biosynthesis, exhibited a > 20-fold decrease in expression with the development of crabapple fruit, and these results were consistent with the trends in the accumulation of the PAs ( Figure 1C ).

To explore the molecular basis of the variations in PA content and obtain a more detailed understanding of the PA regulatory network during crabapple fruit development, the expression of each gene at the fifth developmental stage was compared to that over the consecutive developmental stages and then filtered using |log₂(fold-change)| > 1 or < −1 and a false discovery rate (FDR) < 0.05. The most DEGs were found for S1 vs. S5 (8,710), while the S1 vs. S2 comparison had the fewest DEGs (1,645) ( Figure S3 ).

Many DEGs related to ‘phenylalanine metabolism,’ ‘phenylpropanoid biosynthesis,’ and ‘flavonoid biosynthesis,’ which are associated with PA biosynthesis, were significantly enriched during fruit development, as shown by KEGG analysis ( Figures 2A, B ). Interestingly, several DEGs related to ‘plant hormone signal transduction’ were also enriched throughout fruit development. We also identified several plant hormone signal transduction and response proteins, including auxin-responsive genes, ethylene signal transduction pathway genes and ethylene-responsive genes, among the DEGs related to ‘plant hormone signal transduction’ during fruit development ( Figure 2C , Table S2 ). Interestingly, several auxin-responsive genes were significantly enriched throughout fruit development, which implied that auxin might be involved in PA biosynthesis during fruit flesh development in crabapple.

Notably, ethylene signal transduction pathway genes and ethylene-responsive genes, including ethylene-responsive transcription factors and the ethylene receptor, were upregulated in a stage-specific manner during the later stages of fruit flesh development. In addition, many DEGs related to ethylene biosynthesis and signal transduction, such as ERS (ethylene response sensor) (MD03G1292200), EIN3 (ETHYLENE INSENSITIVE 3) (MD07G1053800), and ERFs (ethylene response factors) (MD04G1009000, MD11G1306500, MD17G1209000, MD16G1140800), were enriched during fruit development, suggesting that ethylene may play an important role in regulating PA biosynthesis in the late stages of fruit flesh development.

To further identify the specific transcription factors involved in regulating PA biosynthesis during crabapple fruit development, a total of 9,471 DEGs were used in a WGCNA analysis, resulting in 17 distinct modules ( Figure 3A ). The MElightcyan module was the highest correlative module with procyanidin B2, and this module included 4,297 genes and had the highest correlation with PA accumulation (0.81) across all developmental stages ( Figure 3B ). In the MElightcyan module, 137 and 296 genes were related to ‘signal transduction mechanisms’ and ‘transcription,’ respectively ( Figure 3B , Tables S3 and S4 ). These analyses further showed that as the fruit developed, ethylene signal transduction pathway genes gradually increased, while ethylene-responsive genes gradually decreased. In addition, genes associated with ‘transcription’ were highly enriched for the MYB, bHLH, ERF, and ARF families. These data indicate that ethylene signal transduction pathway genes or ethylene-responsive genes play an important role in the regulatory network of PA biosynthesis and that MYBs, bHLHs, ERFs, and ARFs may be involved in the metabolism of PAs.

To further identify the specific TFs that participate in regulating PA biosynthesis during crabapple fruit development, 296 genes encoding transcription factors, including members of the MYB, bHLH, ERF, TCP, bZIP, WRKY, and WD40 families, were further analyzed via a correlation network ( Figure 3C ). The top 150 genes that showed the most connections in the network based on their high K_ME (eigengene connectivity) values were defined as hub genes, and the enrichment for MYB, bHLH, and ERF transcription factors were detected.

Through correlation networks and gene expression trends, we identified 12 transcription factors from the MYB, bHLH and ERF families as hub genes ( Figure 3C ). These genes included ethylene-response factors (ERF105, MD07G1248600; ERF023, MD01G1083000; RAP2-4, MD15G1365500; ERF1A, MD04G1058000; ERF5, MD06G1051900; RAV1, MD13G1046100; DREB1A, MD06G1072300; DREB1D, MD04G1067800; ERF1E, MD06G1072200; and ERF061, MD14G1127700), MYB44 (MD08G1107400) and bHLH13 (MD07G1151000) ( Figure 4A ).

To visualize the reliability of the RNA-seq data, the expression level of hub genes were detected by qRT-PCR analysis ( Figure 4B ). The results showed that the expression levels of these genes gradually decreased during fruit development. To better validate the reliability of the selected hub genes, we generated a heat map showing the correlation data. The RNA-seq and qRT-PCR data showed strong correlations (> 0.70) between PAs and the expression of related candidate PA regulators ( Figure 4C ).

The ERF family gene RAP2-4 (RELATED TO APETALA 2-4) participates in regulating plant development and stress resistance via light perception and ethylene signaling (Lin et al., 2008). By contrast, RAV1 (related to ABI3/VP1) is known to be a suppressor involved in flower development, growth, and stress responses (Hu et al., 2004). Here, the expression of RAP2-4 and RAV1 showed the greatest correlations with PAs accumulation among the 12 identified TFs during fruit development ( Figures 4A, C ).

To assess the role of RAP2-4 and RAV1 in PA biosynthesis, 1221 bp PAP2 and 1098 bp RAV1 cDNA sequences from ‘Royalty’ were cloned from fruit flesh, and they were predicted to encode 406 and 365 amino acids, respectively.

Subsequently, the vector constructs 35S::RAP2-4 and 35S::RAV1 (pBI101 for overexpression) were overexpressed in ‘Flame’ tissue culture plants, resulting in stronger DMACA staining and higher PA content in 35S::RAP2-4 expressed leaves than in the control ( Figure 5A ). We also observed weaker DMACA staining and lower PA content in RAV1-overexpressing leaves than in control leaves ( Figure 5A ). Gene expression analysis by qRT-PCR further indicated that the expression of the PA-related genes McLAR1 and McANR1 increased compared to those in the control in RAP2-4-overexpressing leaves, while McANR2 expression was lower in RAV1-overexpressing leaves than in the control ( Figures 5B, C ).

We also overexpressed the RAP2-4 and RAV1 genes in ‘Red Fuji’ fruits to further reveal the role of these two TFs. Agrobacterium tumefaciens cultures containing 35S::RAP2-4 or 35S::RAV1 were individually injected into apple fruits. The fruits infiltrated with 35S::RAP2-4 rapidly accumulated PAs, resulting in deep blue staining in the peels. In contrast, apple fruits infiltrated with the 35S::RAV1 showed weaker DMACA staining. ( Figure 5D ). By using qRT-PCR, the expression of McLAR1 strongly increased compared to control peels in RAP2-4-overexpressing fruit peels, we also observed a significant decrease in the expression of McLAR2 and McANR2 when RAV1 was overexpressed. ( Figures 5E, F ). These results suggest that RAP2-4 may act as an activator in PA biosynthesis, while RAV1 acts as a suppressor.

To verify the speculation that PA biosynthetic genes might be regulated by RAP2-4 and RAV1 in crabapple, a yeast one-hybrid assay was employed to test their ability to bind the promoters of McCHS, McCHI, McF3H, McDFR, McANS, McUFGT, McLAR1, McLAR2, McANR1, and McANR2. The results showed that RAP2-4 bound the promoter of McLAR1, and RAV1 bound the promoter of McANR2. From these results, we deduced that the PA biosynthetic genes McLAR1 and McANR2 might be candidate target genes of RAP2-4 and RAV1, respectively ( Figure 6 ).