Four advanced inbred radish lines “NAU-YH”, “NAU-XLM”, “NAU-XBC”, and “NAU-YZH” were used. The cultivars exhibit red skin-white flesh, green skin-pink-purple flesh, white skin-white flesh and red skin-red flesh, respectively (Figure 1). Seeds were selected and surface sterilized before being germinated on moist filter paper in darkness for 3 days. They were then transplanted into plastic pots containing 1:1 mixture of sterilized soil and peat substrate, and cultured in the greenhouse. The growth conditions included a 14 h light/10 h darkness photoperiod with an average temperature of 18°C. The inability of radish cortex cells to undergo division and expansion results in splitting, an occurrence that is vital for the initiation of taproot thickening. The development of cortex splitting is an important signal of the initiation of taproot thickening growth in radish due to the inability of the cortex cells to undergo division and expansion (Wang et al., 2013). Cortex splitting has been found to begin at around 12 days after sowing (DAS), thus 10 DAS is the pre-cortex splitting, while the peak of root cortex splitting is at 30 DAS. The maximum taproot thickening is achieved at 50 DAS. Subsamples of leaf and taproot issues were collected at 10 DAS (pre-cortex splitting stage), 30 DAS (cortex splitting stage) and 50 DAS (taproot thickening stage). At maturity stage, prior to experiments, radishes were briefly, manually peeled to separate the skin and the flesh, which were then cut into small cubes. Samples in three biological replicates were separated into different batches for anthocyanin and total RNA extraction. Samples for anthocyanin were used immediately, while those for RNA extraction were frozen in liquid nitrogen and stored at −80°C until use.

Individual samples were ground into fine powder in the presence of liquid nitrogen before anthocyanin extraction. The anthocyanin content was measured with the modified method (Mehrtens et al., 2005; Yin et al., 2015). Totally, 2.00 g of the ground powder were dissolved in 6 mL of an extraction solution [concentrated HCl + 80% (v/v) ethanol, 1:99], then extracted by shaking on a mechanical shaker at 110 rpm at room temperature, while shielding from light for 24 h. Then 1 mL of the mature radish extract or 2 mL of seedling (10 DAS) extract was filtered and diluted to 10 mL with 0.4 M sodium acetate (pH 4.5) and 0.025 M hydrochloric acid (pH 1.0) buffers. The absorbance was then observed at 530 and 657 nm on a UV-VIS spectrophotometer (ND752, SPSIC, Shanghai, China). The anthocyanin content (Q) detected was calculated as Q =A530 - (0.25×A657)g FW, where g FW is the fresh mass of the sample in grams. Values, representing means from three independent experiments were subjected to analysis of variance (ANOVA) using PROC GLM code of SAS version 9 (2005). Significant means were separated using the Tukey's Honestly Significant Difference Test at P ≤ 0.05, while Microsoft Office Excel version 2013 was used to generate figures.

The cetyltrimethylammonium bromide CTAB method (Liu et al., 2003) was used for extraction of genomic DNA from mature leaf and root tissues of “NAU-YH”. After being digested with RNase, it was re-extracted with phenol/chloroform/isoamyl alcohol (25:24:1) and chloroform/isoamyl alcohol (24:1). Total RNA from the leaf and root tissues of all the four cultivars at the three previously mentioned developmental stages was extracted according to the previous protocol (Xu et al., 2013). Integrity analysis of DNA and RNA was performed by electrophoresis on a 1% ethidium bromide stained agarose with gel. cDNA was synthesized from 2 μg of RNA using a PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa Bio Inc., Dalian, China). Prior to reverse transcription, DNase was used to remove contaminating DNA according to the manufacturer's instructions. Total RNA was reverse-transcribed using M-MLV reverse transcriptase (Promega, Madison, WI, USA) and an oligo d (T) 18 primer.

For RNA Sequencing Library Construction and Illumina Sequencing, radish (Raphanus sativus L.) advanced inbred line “NAU-YH” was used. Sequencing was done on Illumina HiSeq™ 2500 platform at the Beijing Genomics Institute (BGI, Shenzhen, China). The construction of the library and Illumina sequencing were performed according to a method previously described (Cheng et al., 2013). Radish unigene sequences from the de novo transcriptome data, deposited in the NCBI Sequence Read Archive repository: SRX707630 (Yu et al., 2016) were analyzed to identify and isolate the genes associated with anthocyanin biosynthesis. Using nucleotide sequences from these genes as queries, BLASTx was done against the NCBI Gene Bank (http://www.ncbi.nlm.nih.gov/genbank), radish genome (http://www.nodai-genome-d.org/) (Kitashiba et al., 2014) and “NAU-YH” transcriptome databases. For validation of the transcriptome data, fifteen candidate genes including RsPAL1, Rs4CL3, RsCHS3, RsCHI, RsF3H, RsF3′H1, RsDFR, RsANS, RsANR, RsUFGT, RsTT12, RsSAM, RsOMT, RsGSTU5, and RsGSTF10 were isolated from “NAU-YH”. Gene-specific primers (GSPs) were designed based on the genomic nucleotides and used to amplify the ABGs using mixed root and leaf DNA as the template. Primers used for cloning are shown in Table S1. Each PCR reaction was carried out in a total volume of 20 μL containing 2.0 mM Mg²⁺, 0.2 mM dNTP, 1.0 μM of gene specific primer, 0.8 μL Taq DNA polymerase (TaKaRa Bio Inc., Dalian, China) and 20 ng of diluted DNA template. The conditions for PCR were as follows: 94°C for 3 min; 35 cycles of 50 s at 94°C, 50 s at annealing temperature (Tm) and 90 s at 72°C, and finally extension at 72°C for 10 min. PCR products were purified and cloned into a pMD19-T vector (TaKaRa Bio Inc., Dalian, China). Three independent positive clones from each isolated gene were sequenced on an ABI3730 sequencer (Applied Bio systems, USA). These cloned sequences were aligned with the corresponding unigenes from the transcriptome (https://www.ncbi.nlm.nih.gov/sra/SRX707630/), to access the coverage. The chromosomal location and related predicted sequences were obtained through a manual search from (http://www.nodai-genome-d.org/). ORF Finder and the BLAST (https://www.ncbi.nlm.nih.gov/orffinder/) programs were employed to analyze the cDNA and amino acid sequences. Pfam (http://www.pfam.xfam.org/) was used to identify conserved domains. Coding sequences (CDS) were aligned to DNA sequences and schematics generated using Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/).

Following gene validation, a spatiotemporal analysis of the respective anthocyanin biosynthetic genes' transcript levels, was performed using quantitative real time PCR (RT-qPCR). Sequence-specific primers used for RT-qPCR were designed using Beacon Designer v7.0 (Premier Biosoft International, USA; Table S2). An iQ™ 5 Multicolor Real-Time PCR Detection System (Bio-Rad Laboratories, Berkeley, CA, USA) was used to perform PCR. Each reaction (20 μL) contained 10 μL of SYBR® Premix Ex Taq (Takara), 2.0 μL of cDNA and 0.2 μM of each primer. PCR was carried out under the program of 95°C for 3 min, 40 cycles of 95°C for 5 s, 58°C for 30 s, and 72°C for 10 s. The data were analyzed using iQ™ 5 Optical System Software (version 2.1, Bio-Rad) and expression levels of the gene normalized to RsActin gene (Xu et al., 2012). Relative fold expression changes were calculated using the 2^−ΔΔC_T method (Livak and Schmittgen, 2001). The ANOVA was performed using PROC GLM code of SAS version 9 (2005) and means separated using Tukey's Honestly Significant Difference Test at P ≤ 0.05, while Microsoft Office Excel version 2013 was used to generate figures. Specifically, point analysis was done, depicting the significant variations between the leaf and root tissue, within each gene. Pearson's correlation coefficient, MetaboAnalyst 3.0 (http://www.metaboanalyst.ca/faces/Secure/upload/StatUpload View.xhtml) was used for correlation analysis between gene expression levels and accumulation of total anthocyanin.

Differences in the quantity of anthocyanin were observed in the leaf, young root, flesh and skin of different radish genotypes at different growth stages (Figure 2). At all sampling points, there were no detectable amounts of anthocyanin in the leaf and root tissues of “NAU-XBC”. At 10 DAS significant anthocyanin amounts were only detected in the leaf and root of “NAU-YZH”. There was a steady global increase in anthocyanin content at 30 DAS, wherein, anthocyanin was detected in the leaves of “NAU-YH” and “NAU-YZH”, being highest in “NAU-YZH” while insignificant in “NAU-XLM”. In the root tissues at 30 DAS, anthocyanin was detected in “NAU-YZH”, “NAU-YH”, and “NAU-XLM” with consistently high levels in “NAU-YZH” and lower but equally significant amounts recorded in “NAU-YH” and “NAU-XLM”.

Among the three developmental stages, the maximum value of anthocyanin content was recorded at 50 DAS. At this stage, the root and the skin were separated into the respective skin and flesh components. Consistent with other stages, “NAU-YZH” accumulated considerable amounts of anthocyanin in both tissues. Anthocyanin was only detected in the leaf of “NAU-YZH” and “NAU-YH” being significantly high in the latter cultivar. Sufficiently high amount of anthocyanin was recorded in the flesh of “NAU-YZH” closely followed by “NAU-XLM” but barely detectable in the “NAU-YH” flesh. In the skin at 50 DAS, significant anthocyanin amounts were detected in “NAU-YZH” and “NAU-YH” (NAU-YZH>NAU-YH) (Figure 2).

To identify the genes associated with anthocyanin biosynthesis, the anthocyanin biosynthetic pathway was analyzed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG), which divided this pathway into three distinct phases: the general phenylpropanoid biosynthesis pathway (ko00940), the flavonoid biosynthesis pathway (ko00941) and later the specific anthocyanin biosynthetic pathway (ko00942). According to de novo assembled transcriptome (SRX707630), using local blast search and sequences functionally annotated by the KEGG pathway analysis, 20 genes encoding enzymes of the anthocyanin biosynthesis pathway were identified. These included three PAL (K10775, EC: 4.3.1.24) syntenic genes (14 unigenes), five 4CL (K01904, EC: 6.2.1.12) syntenic genes (17 unigenes), four CHS (K00660, EC:2.3.1.74) syntenic genes (10 unigenes), two F3′H (K00475, EC:1.14.11.9) syntenic genes (11 unigenes) and one gene each for C4H, CHI (K01859, EC:5.5.1.6), DFR (K13082, EC: 1.1.1.219 1.1.1.234), ANS (K05277, EC: 1.14.11.19), F3H (K05280, EC 1.14.13.21) (8, 6, 6, 11 and 11 unigenes, respectively), and 129 unigenes were found to correspond to genes involved in methylation, glucosylation and glycosylation: UFGT (K12338, EC: 2.4.1.298), UGAT (K12937, EC: 2.3.1.254), AT (K12930, EC: 2.4.1.115), MT (K12931, EC: 2.3.1.171), GT1 (K12938, EC: 2.4.1.-) (Table 1).

A total of 102 assembled unigenes were annotated as those that correspond to enzymes involved in the upstream anthocyanin biosynthesis pathway. Specifically, 39 assembled unigenes were found to correspond to the genes in the general phenylpropanoid pathway. From malonyl-CoA to the colored unstable flavonoids, 55 unigene sequences related to six enzymes involved in the flavonoid biosynthetic pathways were isolated. One-twenty nine unigenes were found to be associated with the glycosylation of different anthocyanin aglycons in the specific pathway (Table S4). In the downstream steps of the anthocyanin biosynthesis pathway, stable anthocyanins are trafficked to the anthocyanic vacuole via the non-covalent activity of glutathione S-transferase (GST, EC: 2.5.1.18); 61 unigenes were found to correspond to different GSTs.

Partial DNA fragments or full-length sequences of 15 genes encoding anthocyanin biosynthesis-related enzymes were isolated through cloning. Full length genomic fragments of nine genes (RsCHS3, RsCHI, RsF3H, RsF3′H1, RsANS, RsANR, RsUFGT78D2, RsGSTU5, and RsGSTU10), and partial fragments of six genes (RsPAL1, Rs4CL3, RsSAM, RsOMT, RsTT12, and RsDFR), were isolated. The sequences of these genes were submitted to the National Center for Biotechnology Information (NCBI)/GenBank database under the following accession numbers: RsCHS3 (MF182893), RsCHI (MF182892), RsF3H (MF182895), RsF3′H1 (MF182896), RsANS (MF182899), RsANR (MF182891), RsUFGT78D2 (MF183115), RsTT12 (MF182901), RsGSTU5 (MF182897), RsGSTU17 (MF182898), RsPAL1 (MF285801), Rs4CL3 (MF285800), RsSAM (MF182902), RsOMT (MF182900), and RsDFR (MF182894). The length of these gene sequences varied from 900 to 1,579 bp, and the unigene coverage to corresponding genomic sequences was above 93% (Table 2).

Genomic sequence analysis depicted the different domains of each of the ABGs, which bear the catalytic sites that enable the unique reactions performed by each gene. Intron classification revealed that most of the anthocyanin biosynthetic genes bear phase one introns with predominantly asymmetrical exons, except RsSAM which has one symmetrical (1-1) exon (Figure S2). The analysis also revealed unique features in the flavanone hydroxylase genes; the RsF3′H1 gene located on scaffold 45 is intronless, while RsF3H bears exclusively phase 0 introns. The RsSAM ABG has the highest number of introns (10) and the isoelectric points of radish ABGs are between 4.589 (RsCHI) and 11.87 (RsPAL1), which categorizes the radish ABGs as stable (Table 3).

RT-qPCR was utilized to analyze transcript levels of 14 radish genes associated with anthocyanin biosynthesis, relative to the cultivar, developmental stage and plant tissue. Four color variant radish genotypes (“NAU-YZH”, “NAU-YH”, “NAU-XBC” and “NAU-XLM”) at three developmental stages (10, 30, and 50 DAS) and the leaf, young root, taproot flesh and skin were used. The initial steps of the flavonoid pathway from 4-coumaroyl CoA through chalcone and naringenin to dihydroflavonol are catalyzed by CHS, CHI, F3H, and F3′H (Tanaka et al., 2008). RsCHS3 and RsCHI abundantly expressed in the leaf and root tissues of all genotypes under study, with the highest expression in “NAU-YZH” leaf and “NAU-XLM” root at 10 DAS (Figure 3). The latter gene was 12-folds up regulated at 30 DAS characterized by high levels in “NAU-YZH” leaf and “NAU-YH” root, and non-traceable in all tissues of “NAU-XBC” while there was a general dramatic decline of RsCHI in the root (Figure 3). At 50 DAS, however, the expression levels of RsCHS3 increased by 2-folds with consistent high expression in “NAU-YZH” and “NAU-YH” leaf, flesh and skin, respectively (Figure 3). RsCHI on the other hand generally showed no change in transcript levels with significant high expression in the leaf of “NAU-YH” and the skin and root of “NAU-YZH”. Notably, while the lowest transcripts of RsCHI significantly were expressed in “NAU-XLM” leaf, this gene was increased by 4-folds in the flesh at 50 DAS (Figure 3).

The transcripts of RsF3H were equally expressed in both leaf and root but with significantly higher expression of RsF3′H1 in “NAU-YZH” and “Nau-YH” at 10 DAS (Figure 3). However, RsF3′H peaked in the leaf and was barely detected in the root at 30 DAS, while RsF3H remained consistently up regulated with high expression in “NAU-YZH” leaf and “NAU-YH” root (Figure 3). There was no change in transcript levels of RsF3H but RsF3′H1 globally increased by 10-folds, and both genes showed a consistent high expression in the leaf and skin of “NAU-YZH” and “NAU-YH” respectively, at 50 DAS (Figure 3). Transcripts of RsF3H were notably significantly repressed in the leaf and skin but increased by 5-folds in the flesh of “NAU-XLM”.

In the synthesis of dihydroflavanols, which are intermediates in anthocyanin biosynthesis, the genes RsDFR and RsANS play critical synthesis functions, while RsANR acts as an inhibitory gene to RsANS. The RsANS gene was down-regulated in the leaf of “NAU-XBC”, and up-regulated in “NAU-YZH” and “NAU-YH” with no significant variation in the root tissues at 10 DAS (Figure 3). Transcript levels increased in the leaves, against traceable levels in the root at 30 DAS and continued to increase at 50 DAS while remaining highly expressed in “NAU-YZH” leaf, flesh and skin, “NAU-XLM” flesh and “NAU-YH” skin (Figure 3). In all genotypes under study, transcript levels of RsDFR were not variable in the 10 DAS and 30 DAS leaves and 10 DAS root but with the exception of “NAU-XBC”, decreased in the root at 30 DAS (Figure 3). However, RsDFR expression steadily increased at 50 DAS with the highest levels depicted in the “NAU-YZH” and “NAU-XLM” flesh and “NAU-YZH” and “NAU-YH” skin (Figure 3). There was a high expression of RsANR in “NAU-YH” leaf and over 10-fold expression in “NAU-XBC” leaves. Over 6-fold increased expression was recorded in the leaf of “NAU-XLM” but silenced in the root tissues at 50 DAS (Figure 3).

RsSAM, RsOMT, and RsUFGT are involved in glucosylation and methylation of anthocyanidins, resulting in stable compounds. Transcripts of RsSAM were highly expressed in tissues of “NAU-YZH” and the colored tissues of “NAU-YH” and “NAU-XLM”, and undetectable in all tissues of “NAU-XBC” at 30 and 50 DAS (Figure 4). RsOMT on the other hand was significantly expressed in the “NAU-YZH” and “NAU-XLM” with higher expression in the root and the leaf at 10 and 30 DAS, respectively. Notably, transcripts of RsOMT were significantly elevated at 50 DAS in the white colored “NAU-XBC” (Figure 4). Nevertheless, RsUFGT consistently increased across the three developmental stages, with consistent high expression in the colored radish cultivars, although the expression dramatically decreased in the root of “NAU-YH” at 30 DAS (Figure 4). Generally, the transcripts of RsSAM, RsOMT, and RsUFGT were upgraded in tissues expressing high anthocyanin content.

There was no distinct variation in the glutathione transferase transcripts at 10 DAS, but their levels heightened at 30 DAS in the leaf and were barely detectable in the root at this stage (Figure 4). However, RsGSTF11 was over 15-folds up-regulated in the leaf, with a constant up-regulation in all tissues of the “NAU-YZH” at 50 DAS. On the other hand, RsGSTU5 and RsGSTF10 transcript levels increased in the skin with an up-regulation in the “NAU-XBC” and “NAU-YH” flesh and skin, respectively (Figure 4).

The genotype “NAU-YZH” was used to analyze the relationship between transcript levels and total anthocyanin content. Three sampling points, three biological replicates and the leaf, young root (10 and 30 DAS), leaf, taproot skin and flesh tissues (50 DAS) were regarded as independent factors in the transcript pairwise comparisons of total anthocyanin, yielding 21 transcriptional data points for each gene (Figure 5). Pearson's correlation coefficient was used for the correlation analysis.

Transcripts of RsF3H, RsUFGT, RsANS, and RsF3′H1, showed a significant correlation (r > 0.85) with total anthocyanin content (Table S3).The co-expression patterns between the gene expression and anthocyanin contents were generated and are depicted in a heat map of clustered correlations (Figure 5). The set of genes from number 8–11 had the strongest positive correlation with total anthocyanin, including RsF3H (r = 0.92), RsUFGT (r = 0.92), RsF3′H1 (r = 0.84), and RsANS (r = 0.89) while those in lines 1–4 exhibited the strongest negative correlation: RsGSTU5 (r = 0.62), RsSAM (r = 0.5), RsGSTF10 (r = −0.47), and RsANR (r = −0.47) (Table S3). Linear regressions performed between the cumulative transcription of each of the 14 genes and the corresponding anthocyanin contents were shown in Figure S3. These data indicated that RsF3H, RsUFGT, RsANS, and RsF3′H might be candidate genes ascribed to the red pigment of the “NAU-YZH” root.

The anthocyanin concentrations of radish under study increased with increase in days after sowing, resulting in higher concentrations in mature taproots than sprouts of colored tissues. There were variations in the levels of anthocyanin in the skin and flesh organs of the four different radish cultivars. The increase in anthocyanin accumulation as the crop advances in growth could be due to increased accumulation of anthocyanin content in tissues, resulting from increase in size and hence the capacity of the plant's biosynthetic machinery. It is speculated that this occurrence is associated with the coordination of plant metabolism because of the ability of plants to regulate the metabolism in respective organs at various developmental stages (Majdi et al., 2014). These results concur with those found in yam, cherry, grape, and chili (Aza-Gonzalez et al., 2013; Liu et al., 2013; Xie et al., 2015; Yin et al., 2015). There was also a difference in the accumulation of anthocyanin in tissues of different radish at each developmental stage. For instance in “NAU-YH” anthocyanin in the skin was higher than in the flesh, while “NAU-XLM” recorded higher anthocyanin content in the flesh and significantly low amounts in its skin, but these contents increased progressively with growth, though at varying rates.

The length of cloned gene sequences from our study varied from 900 to 1,579 bp, and the unigene coverage to corresponding genomic sequences was above 93%, indicating that the radish anthocyanin unigenes from RNA-seq were successfully assembled and were viable for further investigation.

Gene structure analysis predicted the significant domains in each gene, which help to explain their mechanism of action. For instance, the distinct substrate specificity of the NAD_binding _4 domain of RsDFR with ASN-active catalytic site explains the preferential accumulation of pelargonidin and cyanidin derivatives instead of malvidin or petunidin (Shimada et al., 2005). The analysis also revealed that RsF3′H is intronless and that most of the ABGs bear phase one introns with predominantly asymmetrical exons, except RsSAM which has one symmetrical (1-1) exon. Transcript abundance of RsF3′H was found to be consistently lower when compared to RsF3H. Introns have been found to contribute to increased protein abundance while elimination of introns renders protein product undetectable. This impeccable intronic function is referred to as intron-mediated enhancement (Akua et al., 2010). Extentive analyses further confirmed these inferences where genes bearing introns in yeast were found to have enhanced mRNA and protein levels than intronless genes (Juneau et al., 2006). The pattern of spliceosomal introns has a significant relationship with the conservation of splice signals sequence of exons, therefore, the relatively dominated phase 1 introns in the ABGs in radish indicates a relative reduction in conservation of these genes (Long and Deutsch, 1999).

Despite the consistent occurrence of all ABGs in any studied variety, a species-specific control of the genes of the basic pathway and at key branching points is assumed to contribute to the differences in anthocyanin content and the shift from lighter to darker hues as the crop develops (Ageorges et al., 2006). In the present study, variations in anthocyanin coloration across tissues of four radish genotypes are ascribable to alterations in the unique expression patterns of the overall set of anthocyanin genes.

The transcripts of ABGs portrayed a commensurate expression for most of the genes in the leaf and root tissues of the four radish cultivars at 10 DAS, revealing that ABGs are coordinately expressed at early stages of development. However, the expression of anthocyanin pathway genes is consistent with the accumulation of anthocyanins: as the plant begins to pigment, the expression of key anthocyanin biosynthesis genes is detected to obviously increase.

Chalcone synthase (CHS) and Chalcone isomerase (CHI) are the first genes in the flavonoid branch of the anthocyanin biosynthesis pathway. RsCHS3 recorded the highest level of transcripts when compared to all the other genes and consistently increased throughout plant development. Its transcript abundance correlated positively with total anthocyanin, signifying that it is a key gene for anthocyanin synthesis, findings that echo previous research (Park et al., 2011; Chen et al., 2016). The loss of the RsCHS3 gene from the 30 DAS stage, through 50 DAS, in the “NAU-XBC” could perhaps be the single most factor contributing to the loss of color in this cultivar. A knockout mutation of the gene impeded anthocyanin accumulation in seeds resulting in a transparent testa in Arabidopsis (Shirley et al., 1995). Although RsCHI correlated positively with anthocyanin content, its transcription levels were markedly down- regulated in the root at 50 DAS, being non-consistent with previous reports (Xu et al., 2014). This result indicates that a single gene is not responsible for anthocyanin accumulation and that anthocyanin biosynthesis involves the coordinated mechanism of many genes (Walker et al., 2007; Yu et al., 2012).

The transcription of RsF3′H1 seems to have been developmentally activated after pre-cortex splitting (30 DAS) in “NAU-YZH”, “NAU-YH” skin, and “NAU-XLM” flesh, but with lower transcript levels than RsF3H and was barely transcribed in “NAU-XBC”, which is white colored. However, the expression of RsF3H and RsF3′H1 greatly increased in all cultivars at maturity except for the non-pigmented, taproot skin of “NAU-XLM”, taproot flesh of “NAU-YH” and tissues of “NAU-XBC”. In the latter cultivar, the levels of transcripts remained low but detectable. Striking variations were also observed in the regulation of the flavonoid hydroxylase genes. The expression profile of RsF3H was relatively high even before 50 DAS and transcripts of this gene were present at maturity in all cultivars, including the white cultivar “NAU-XBC”.

The reciprocal expression levels of these flavonoid hydroxylase genes most likely contribute to different color hues in radish. Previous studies also proposed a plausible cause of color transition from pelargonidin type in “NAU-YZH” to cyanidin type in “NAU-XLM” which could be due to a reduced RsF3′H1 mRNA level (Mudalige-Jayawickrama et al., 2005) resulting in the inefficient production of dihydroquercitin and increased dihydrokaempferol and subsequently, pink-purple coloration. Differential expression of transcripts encoding flavonoid-hydroxylase was also reported in mulberry (Morus alba L), cauliflower (Brassica oleraceae var botrytis), and Chinese bayberry (Myrica rubra) (Chiu et al., 2010; Niu et al., 2010; Li et al., 2014). Notably, RsF3H and RsCHI were clearly expressed in the white-flesh radish “NAU-XBC”, but anthocyanin accumulation was not detected. These findings suggest that RsF3H and RsCHI are likely more highly regulated than the other anthocyanin biosynthesis structural genes. Additionally, findings from other studies indicate that CHI and F3H are early genes in the anthocyanin biosynthetic pathway and are coordinately expressed, with increased transcript levels toward plant maturity (Ravaglia et al., 2013).

Dihydroflavonol 4-Reductase (DFR) encoded by a single gene (DFR) converts dihydroflavonols to leucoanthocyanidins. Besides the high expression in colored tissues, this gene was also found to be well expressed in the white tissues of “NAU-XBC”, and is, therefore, likely to participate in the synthesis of other secondary metabolites like free-auxins (Shen et al., 2013). These results relate to those obtained in litchi and kiwifruit (Montefiori et al., 2011; Wei et al., 2011).

Anthocyanidin synthase/leucoanthocyanidin oxidase (LDOX) catalyzes the conversion reaction of leucoanthocyanidins to colored anthocyanidins. However in this study, in spite of the high expression of RsANS in the early developing stage, the leaf color of three colored radishes is green. This could result from the catalysis of Anthocyanidin synthase (ANS) substrate into flavonol or/and the inhibitory role of the RsANR, which converts the colored anthocyanidins to epicatechin, resulting in redirection of anthocyanin pathway into proanthocyanidin pathway (Jaakola et al., 2002). Our study is consistent with findings from bilberry (Jaakola, 2013). However in the later stages ANS is highly expressed in colored tissues and at stages coinciding with elevated anthocyanin amounts, consistent with previous findings (Zhang et al., 2015).

Down-stream genes such as glucosyltransferases greatly influence the direction of anthocyanin synthesis through the regulation of anthocyanidin glucosyltransferase in thickening radish taproot.

It was found that the transcription level of RsUFGT was much lower in the white-flesh cultivar “NAU-XBC” than in the red and pink colored radishes. Similar to RsCHS3, our results also showed that the mRNA levels encoding RsUFGT, the specific gene for anthocyanin biosynthesis, increased proportionally to the anthocyanin content across all the three developmental stages, suggesting that these two genes are under a different regulatory regime, in comparison to the other ABGs in radish, and that the biosynthesis of anthocyanin is controlled at an earlier stage as reported in previous studies (Kobayashi et al., 2001). The expression of O-Methyltransferase (RsOMT) in the red colored “NAU-YZH” was higher than that of RsUFGT, which was also well expressed in all genotypes, unlike RsOMT. The purple-pinkish colored “NAU-XLM” expressed comparatively lower levels of OMT when compared to RsUFGT. It has been proposed that the methylation of B-ring hydroxyl groups causes a shift toward deeper red colors (Tanaka et al., 2008). Therefore, it can be inferred that the relative abundance of RsF3H to RsF3′H1 and RsOMT to RsUFGT could per se explain to a greater extent the phenotypic variation of anthocyanin content, color hue and color intensity in radish. Secondly, the key regulation point for quantitative anthocyanin variation is in the downstream pathway at the RsUFGT level, but the qualitative differences are precisely controlled upstream of RsUFGT at the flavonoid hydroxylases' level and at RsOMT which is downstream of RsUFGT.

The sequestration of anthocyanin from the cytoplasm to the vacuole is poorly understood, and various mechanisms have been put across to explain the process, including transport proteins like GSTs and MATE transporters. The GSTs encoded by a group of Glutathione-S-transferase genes, whose specific functions remain to be elucidated in radish, have been implicated in the transport of anthocyanin in other crops (He et al., 2010; Gomez et al., 2011). In this study, the featured GSTs were found to positively coincide with total anthocyanin and exhibit similar expression patterns to the RsUFGT, although the correlation between cumulative transcription and total anthocyanin throughout crop development was lower. The abundance of GST transcripts in the present study exhibited genotypic specificity, as RsGSTF10 was found to correlate with elevated anthocyanin in “NAU-XLM”, RsGSTU5 with “NAU-YZH” and RsGSTF11 with “NAU-YZH”, “NAU-YH” and “NAU-XLM”. RsGSTU5 may contribute to spatial differential accumulation in red radish, owing to its elevated transcript levels in the root against the suppressed transcript levels of key ABGs at 30 DAS. In anthocyanin transport, it has been suggested through mutant analysis that the anthocyanin defective mutants were unable to accumulate anthocyanins into the vacuoles (Conn et al., 2008) implying that glutathione transferases are possible anthocyanin transporters.

At 30 DAS, the red colored radish accumulated significant amounts of anthocyanin in the root despite the down regulation of major ABGs in this tissue. It was also found that TT12, a MATE transporter was consistently up-regulated in the red colored radish but low amounts in the white colored tissues. TT12 was reported to mediate anthocyanin transportation in Arabidopsis (Marinova et al., 2007).

To our knowledge, this is the first report describing the spatial-temporal expression patterns of anthocyanin biosynthetic genes (ABGs) in radish. Our results demonstrate the coordinated expression of ABGs in relation to anthocyanin accumulation in radish tissues and that there may be a common regulatory mechanism governing the coordinated expression of related genes. Furthermore, it appears that the major control point to anthocyanin biosynthesis in radish is UFGT. Globally, the correlation of anthocyanin content with coordinated gene regulation would be the key contributing factor to phenotypic and spatial-temporal anthocyanin accumulation in radish.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.