Comparative transcriptome analysis identifies genes associated with chlorophyll levels and reveals photosynthesis in green flesh of radish taproot

The flesh of the taproot of Raphanus sativus L. is rich in chlorophyll (Chl) throughout the developmental process, which is why the flesh is green. However, little is known about which genes are associated with Chl accumulation in this non-foliar, internal green tissue and whether the green flesh can perform photosynthesis. To determine these aspects, we measured the Chl content, examined Chl fluorescence, and carried out comparative transcriptome analyses of taproot flesh between green-fleshed “Cuishuai” and white-fleshed “Zhedachang” across five developmental stages. Numerous genes involved in the Chl metabolic pathway were identified. It was found that Chl accumulation in radish green flesh may be due to the low expression of Chl degradation genes and high expression of Chl biosynthesis genes, especially those associated with Part Ⅳ (from Protoporphyrin Ⅸ to Chl a). Bioinformatics analysis revealed that differentially expressed genes between “Cuishuai” and “Zhedachang” were significantly enriched in photosynthesis-related pathways, such as photosynthesis, antenna proteins, porphyrin and Chl metabolism, carbon fixation, and photorespiration. Twenty-five genes involved in the Calvin cycle were highly expressed in “Cuishuai”. These findings suggested that photosynthesis occurred in the radish green flesh, which was also supported by the results of Chl fluorescence. Our study provides transcriptome data on radish taproots and provides new information on the formation and function of radish green flesh.


Introduction
Non-foliar plant organs such as the fruit of tomatoes, oranges, eggplants, and bananas, are commonly green in youth and lose chlorophyll (Chl) at maturity [1]. In contrast, the skin and flesh of the green radish is green throughout the developmental stages, remaining green even after harvesting for up to 6 months at low-temperature (-0.8 to -1˚C) [2]. Therefore, the greenfleshed radish is an ideal subject for investigating the formation, retention, and function of the interior green part in non-foliar organs. The green color of the plant is due to the presence of Chls, which are pigments essential for photosynthesis. The Chl metabolic pathway has been well studied in plant leaves and contains three different phases: Chl a biosynthesis, Chl cycle, and Chl degradation. Although Chl synthesis occurs in plastids, the enzymes involved in Chl a biosynthesis are encoded by nuclear genes in higher plants [3]. The Chl a synthesis pathway can be divided into four parts [4]. Part Ⅰ is the reaction of glutamate to 5-aminolevulinic acid (ALA, the precursor of Chls), with three enzymes (EARS, hemA, and hemL) involved in this process [5,6]. Part Ⅱ involves reactions from ALA to uroporphyrinogen Ⅲ, which are catalyzed by hemB, hemC, and hemD [7]. Part Ⅲ involves the reactions from uroporphyrinogen Ⅲ to protoporphyrin Ⅸ, with the aid of three enzymes-hemE, hemF, and hemY [7]. Part Ⅳ involves the reactions from protoporphyrin Ⅸ to Chl a. The first step in Part Ⅳ is the insertion of Mg 2+ into protoporphyrin IX. This step is catalyzed by Mg-chelatase, which possesses three distinct subunits, chlD, chlI, and chlH [8]. In addition to Mg-chelatase, at least six other enzymes are also responsible for the formation of Chl a in part Ⅳ, including chlM, chlE, POR, DVR, chlG, and chlP [7]. The interconversion between Chl a and Chl b is known as the Chl cycle, which is the second phase of Chl metabolism. Five enzymes are involved in the Chl cycle, CAO, chlG, CLH, NOL (NYC1), and HCAR [9,10]. Chl degradation is the third phase of Chl metabolism, with SGR, CLH, PAO, PPD, and RCCR participating in this phase [11].
Excessive Chls or Chl intermediates can produce reactive oxygen species and harm biomolecules, while too few Chls will decrease the efficiency of photosynthesis. Thus, the Chl metabolism pathway is precisely regulated in plants, and the expression of genes involved in Chl metabolism may differ greatly depending on the organs, developmental stages, time, and environmental factors relating to each plant [12]. For example, Chl metabolic genes encoding hemA, chlH, chlE, and POR were found to be closely related to Chl content in the petals and leaves of chrysanthemum [13], while the gene encoding Mg-chelatase played an important role in the green petals of the carnation [12]. Compared with that of the golden variant of the fruit, the green kiwifruit exhibits a lower expression of SGR2 during the maturation stages [1]. In tea leaves, Li et al. [14] found that the expression levels of PORA were significantly correlated with Chl a content, and the expression of two NYC1 genes showed a significant correlation with Chl a and Chl b content. A recent study on the skin of "Xinlimei" radish taproot showed that several genes associated with Chl metabolism presented significant higher expressions during the taproot skin greening process, such as hemA, chlH, POR, chlG, HCAR, CAO, NOL, CLH, and RCCR [15].
Green plant tissues containing Chls are known to use solar energy for photosynthesis. To date, various strategies have been developed to measure and detect photosynthesis, among which Chl fluorescence is a popular method and is an in vivo photosynthesic probe [16]. Chl fluorescence is effective for monitoring photosynthesis occurring in plant internal tissues, where photosynthesis is difficult to measure accurately [15,17]. Furthermore, the genes involved in photosynthesis are often used as markers for plant photosynthesis, such as RBCS, RBCL, PEPC, and LHCB [1,18]. Although the photosynthetic performance of non-foliar organs and tissues has been widely assessed, limited research is available on the green flesh of radish taproots.
In this study, RNA sequencing (RNA-seq) was performed to compare gene expression profiles between the flesh of green-fleshed "Cuishuai" (GF) and white-fleshed "Zhedachang" (WF) radish varieties during taproot development. We comprehensively explored the differentially expressed genes (DEGs) in green and white flesh tissues during radish taproot development, targeting Chl metabolism and photosynthesis-related pathways. Chl contents, Chl fluorescence, and light response curves of the taproots of these two varieties were also measured. Our results form a sound basis for further studies on the formation, retention, and function of green flesh in radish taproots at the transcriptional level.

RNA extraction, library construction and sequencing
Total RNA was isolated using the Trizol reagent (Promega, Madison, WI, USA), after which a Nano Photometer1 spectrophotometer (IMPLEN, CA, USA) was used to check for the purity of the RNA. Next, RNA concentration was measured using the Qubit1 RNA Assay Kit in a Qubit1 2.0 Fluorometer (Life Technologies, CA, USA), and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total of 3 μg RNA per sample was used as the input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext1 Ultra™ Directional RNA Library Prep Kit for Illumina1 (NEB, Ipswich, MA, USA) according to manufacturer recommendations. Finally, 30 RNA libraries were constructed with the library quality assessed using the Agilent Bioanalyzer 2100 system. Each cDNA library was sequenced on an Illumina Hiseqxten-PE150 platform. Raw RNA-seq data have been deposited in the Sequence Read Archive (SRA) in NCBI. The BioProject number is PRJNA684971.

Bioinformatics analysis
Raw RNA-seq reads were first filtered by removing adapter sequences, low-quality sequences (Sanger base quality�20), and sequences with N content greater than 5%. The clean reads obtained were then mapped to the Radish reference genome (ftp://ftp.kazusa.or.jp/pub/radish) using HISAT2 software [20]. Subsequently, the featureCounts software was utilized to assign sequence reads to the genomic features [21]. DEGs between different samples were determined based on the DEGSeq R package [22] with a threshold of |log 2 FoldChange|�2 and adjusted P-value�0.01. The Benjamini and Hochberg [23] method was applied for the P-value adjustment.
All sequences were aligned to protein databases (UniProtKB/Swiss-Prot, NR, GO, KEGG) using BLASTX with an E-value�10 −5 to acquire their functional annotations. Gene Ontology (GO) analysis was performed using Goseq [24] based on three groups: molecular functions, biological processes, and cellular components. The GO terms with adjusted P-value�0.05 were considered as significantly enriched terms. For Kyoto Encyclopedia of Genes and Genomes (KEGG)-based annotation analysis, the R/ClusterProfiler package [25] was used to test the statistical enrichment of DEGs of pathways with adjusted P-value�0.05.

Quantitative real-time PCR analysis
Gene-specific primers were designed using the Primer Premier 3.0 plus. The actin gene was used as an internal control to standardize the results. Primer pairs used in the qRT-PCR are shown in S1 Table. The qRT-PCR analysis was carried out using SuperRealPreMix Plus (SYBR Green) on a Bio-Rad iQ5 Real-Time PCR platform. The qRT-PCR reaction was initially performed at 95˚C for 15 min, followed by 40 cycles of 10 s at 95˚C and 30s at 60˚C. Melting curves were used to verify primer specificity, and all reactions were performed in four replicates. The comparative Ct method (2 −ΔΔCT ) was used to calculate the relative expression levels of the target genes.

Auto-fluorescence, photosynthetic activity, and Chl content in radish flesh
GF is a green radish with dark green skin and green flesh (Fig 1A), while the WF has both white skin and flesh ( Fig 1E). Chl auto-fluorescence was observed in the GF (Fig 1C), however, it was absent in the WF (Fig 1G and 1g). Furthermore, the FV/FM, FPSII and ETR of the GF were measured using the IMAGING-PAM. The FV/FM value was 0.694 for the GF (Fig 2A). With an increase in photosynthetically active radiation (PAR), FPSII significantly decreased and was completely suppressed at 308 μmol�m -2 �s -1 , while ETR increased and reached a maximum value of 3.575 at 161 μmol�m -2 �s -1 , before decreasing ( Fig 2B).
The flesh tissues of GF and WF were collected from each of the five developmental stages (S1-S5), and Chl content was measured using the spectroscopic method. Only trace amounts of Chls were detected in the flesh tissues of WF, whereas numerous Chls were identified in the GF across all five stages of development (Fig 3). In the GF, Chl a, Chl b and total Chl content remained at a high level across all stages. In the WF, no significant differences were observed for Chl a, Chl b, and total Chl content across the five stages. It is worth noting that Chl a/b ratio did not change significantly during the developmental period. It was 0.98-1.48 in GF and 0.44-0.57 in WF.

RNA-seq analysis
Radish flesh tissues of five corresponding stages from GF and WF were used to construct 30 independent libraries, which were then sequenced using the Illumina Hiseqx-ten-PE150 platform. The sequencing results are presented in Table 1   Validation of the RNA-seq data by qRT-PCR Ten DEGs were randomly selected to confirm the RNA-seq results using qRT-PCR. The expression patterns of selected genes were similar to the results of RNA-Seq (Fig 4A-4J) and a significant correlation was found between them (R 2 = 0.728, P = 3.591e-15) ( Fig 4K). This showed that the results of RNA-seq analysis were consistent with those of qRT-PCR analysis. The qRT-PCR and RNA-seq values for the target genes are listed in S2 Table.

DEGs analysis
To obtain the DEGs between GF and WF, the GF was compared to the WF at the corresponding stage (GFS1 vs. WFS1, GFS2 vs. WFS2, GFS3 vs. WFS3, GFS4 vs. WFS4, and GFS5 vs. WFS5). In summary, 8416 DEGs between GF and WF were detected at five stages (S3 Table, Fig 5). A total of 4161, 3887, 4183, 4189 and 4074 DEGs were obtained for S1, S2, S3, S4 and S5, respectively ( Fig 5A). It should be noted that a total of 1367 DEGs were differentially expressed across the five stages (Fig 5B), of which 609 genes were up-regulated and 742 genes were down-regulated ( Fig 5C).
We performed the GO functional enrichment analysis using GOseq software to further understand the functions of these DEGs. To summarize, 180 enriched GO terms were identified across the five stages (Fig 6 and S4 Table). At S1, the number of enriched GO terms was 118, which was the highest across five stages. A total of 108 GO terms were identified at S2. At S3, only 83 GO terms were identified. At S5 and S6, 91 and 106 were identified, respectively. More than 32% of the GO terms [58] were simultaneously enriched across all five stages ( Fig  6A). Photosynthetic-related terms, such as photosynthesis, photosystem I, photosystem II, light harvesting in photosystem I, chloroplast, chloroplast membrane, and chloroplast thylakoid, were included in these 58 terms (Fig 6B). The GO terms of 609 constantly up-regulated and 742 persistently down-regulated DEGs across the five stages were also analyzed. With regard to DEGs constantly up-regulated in GF compared with WF, chloroplast, chloroplast thylakoid membrane, chloroplast envelope, thylakoid, chloroplast thylakoid, photosynthesis and chloroplast stroma were the seven most enriched GO terms; and chloroplast, chloroplast stroma, chloroplast envelope, chloroplast thylakoid membrane, protein-chromophore linkage, response to red light and response to far-red light were the seven most enriched GO terms among DEGs persistently down-regulated in GF vs. WF (S2 Fig). The results of GO enrichment showed that many DEGs could function in the chloroplast or on the chloroplast envelope and may play important roles in the photosynthetic process.
Additionally, we performed KEGG pathway enrichment analysis of these DEGs. In summary, 273 DEGs were significantly enriched in ten KEGG pathways among the different stages (P value�0.05) (S5 and S6 Tables). There were five pathways enriched at all five stages, including photosynthesis (ko00195), photosynthesis-antenna proteins (ko00196), glyoxylate and dicarboxylate metabolism (ko00630), carbon fixation in photosynthetic organisms (ko00710), and porphyrin and Chl metabolism (ko00860). The results of KEGG pathway enrichment analysis were similar to those of the GO enrichment.

Chl metabolic pathway
Many genes that participate in the Chl metabolic pathway have already been identified and characterized in higher plants. We identified a total of 74 genes involved in this pathway in the present study, among which 34 genes showed significantly differential expression levels between GF and WF (Fig 7). The transcription levels of 12 genes were continuously higher in GF than those in WF from S1 to S5, including five chlP genes (Rsa1.0_00104.1_g00052.    01, S3 Fig). The 12 DEGs with continuously high expression levels in GF from S1 to S5 were all significantly positively correlated with the Chl content. These genes may play important roles resulting in the differences in Chl content between GF and WF.

Carbon fixation pathway
A total of 114 genes participating in the carbon fixation pathway were identified in this study, which included 41 DEGs (Fig 8 and S7 Table). Carbon reactions of photosynthesis, including the C 3 carbon cycle, C 4 carbon cycle and crassulacean acid metabolism, convert CO 2 into carbohydrates. The C 3 carbon or Calvin cycle is the most essential pathway for CO 2 assimilation. We found that 29 DEGs encoding 12 key enzymes were enriched in the C 3 carbon cycle. Among them, there were 25 significantly up-regulated genes, including the orthologs of RBCS, PGK, GAPA, ALDO, FBP, PRK and RPE. In addition, 12 DEGs were enriched in the C 4 carbon

PLOS ONE
cycle. Of these, six genes had higher expression levels in the GF than those in the WF at least one of the five stages. Three PPDK genes and six PEPC genes (key genes of the C 4 pathway) were also identified in our transcriptome data. All PPDK genes and two PEPC genes were expressed at extremely low levels in the GF. Although another four PEPC genes were expressed at the middle levels, no difference was detected in their expression levels between GF and WF tissues. Therefore, it is unlikely that the C 4 pathway exists in the radish green flesh.

Photosynthesis-related Pathways
A total of 64 genes specifically involved in the photosynthesis pathway were differentially expressed between GF and WF tissues in the present study, including 23 genes associated with the PSII core, 23 genes associated with the PSI core, two genes associated with the cytochrome b6/f complex, 13 genes associated with photosynthetic electron transport and three genes associated with F-type ATPase (Fig 9A). Among these 64 genes, a vast majority of the genes (up to 75%) were significantly up-regulated at all five stages in GF tissues.
The antenna system of photosynthesis consists of proteins encoded by light-harvesting Chl a/b binding (LHC) genes. These LHC genes are divided into two sub-classes: LHCA associated with photosynthesis I (PS I) and LHCB associated with photosynthesis II (PS II). A total of 36 LHC genes, including 11 LHCA genes (three LHCA1, four LHCA2, three LHCA3 and one LHCA4) and 25 LHCB genes (eight LHCB1, five LHCB2, two LHCB3, three LHCB4, three LHCB5 and four LHCB6), were differentially expressed between GF and WF during the process of taproot development (Fig 9B). Thirty-three of the above 36 genes showed remarkably higher expression levels across different stages in GF than those in WF tissues; Rsa1.0_02580.1_ g00002.1 (LHCA2) mRNA content in GF was significantly higher than that in WF tissues at S4, and Rsa1.0_01789.1_g00003.1 (LHCA2) and Rsa1.0_02221.1_g00007.1 (LHCB4) were significantly up-regulated at all five stages except S4.
The circadian rhythm pathway was also significantly enriched according to the comparative transcriptome data of GF and WF. We identified 16 DEGs involved in this pathway, of which eight genes belonging to two families were transcription factors (Fig 9D), such as six HY5 genes (bZIP transcription factors) were identified. Among them, Rsa1.0_00016.1_g00067.1, Rsa1.0_00195.1_g00004.1, Rsa1.0_02208.1_g00006.1, and Rsa1.0_12036.1_g00001.1 showed significantly lower expression levels in GF than those in WF during the stages from S2 to S5, while Rsa1.0_03795.1_g00001.1 expression was significantly lower in GF than that in WF only at S4. The expression of the other HY5 (Rsa1.0_02116.1_g00001.1) was complicated, and was higher in GF than that in WF at S1 but lower at S2, S3, S4 and S5. One TCP transcription factor gene (Rsa1.0_00058.1_g00008.1) was expressed at higher level in GF than that in WF at S2 and S4, and the other (Rsa1.0_05545.1_g00002.1) only at S4. Besides the seven transcription factor genes, eight functional DEGs relevant to the circadian rhythm pathway were also identified. Of these, five genes including PRR7, CRY2, CHS and CSNK2B (Rsa1.0_30426.1_g00001.1 and Rsa1.0_50931.1_g00001.1) were expressed at higher levels in GF than those in WF. However, PHYA, PHYB and CSNK2B (Rsa1.0_05900.1_g00002.1) showed lower expression levels in the GF.

Validation of candidate genes in different radish cultivars
Based on the correlation analysis (S3 Fig) and expression levels of DEGs (Fig 7), we selected five genes that participated in part Ⅳ of Chl biosynthesis as candidate genes for further validation in ten white-and green-fleshed cultivars (S8 Table) using qRT-PCR. These five genes were all highly expressed in green-fleshed cultivars, but showed low expression in whitefleshed cultivars (Fig 10). Therefore, our results confirmed that chlP, POR, chlM and chlE are potential candidate genes responsible for the higher Chl content in green flesh than that in white flesh samples from radish taproots.

Discussion
To understand the molecular mechanism of green flesh formation in radish, we extracted RNA from taproots of two representative cultivars, GF (green-fleshed radish) and WF (white-fleshed radish), across five developmental stages (S1-S5). RNA sequence analysis was subsequently carried out. Furthermore, differences in gene expression levels were analyzed between GF and WF at the corresponding stages. GO annotation showed that genes related to chloroplast, chloroplast membrane, chloroplast thylakoid, Chl biosynthetic process, and photosynthetic electron transport were highly enriched among the up-regulated genes in GF. The KEGG pathway analysis also supported the results of GO enrichment, which showed that five photosynthesis-related pathways (photosynthesis, photosynthesis-antenna proteins, glyoxylate and dicarboxylate metabolism, carbon fixation, and porphyrin and Chl metabolism) were enriched at all five stages. In addition, we observed clear Chl auto-fluorescence in GF tissues. These results indicate that the green flesh of taproots may contain chloroplasts, which may be active in GF tissues.
Chl is a major photosynthetic pigment responsible for the green color of the plants. Most of the genes participating in this pathway have been identified and characterized in the leaves of higher plants [4]. Some studies have also been carried out in non-foliar tissues, such as kiwi flesh [1], tomato fruits [26], green flower petals [13], radish taproot skin [15], barley developing spikelets [27], and litchi pericarps [28]. In the current study, we found 74 genes, including 34 DEGs, involved in Chl metabolism in radish flesh (Fig 7). Pearson's correlation analysis showed that the expression levels of approximately 80% of above 34 DEGs were significantly correlated with the Chl content (P�0.01) (S3 Fig). However, only 12 genes, chlI (1), chlM (1), chlE (2), POR (2), chlP (5), and CAO (1), showed higher expression in GF than those in WF from S1 to S5. Interestingly, 11 of these 12 genes were involved in Part Ⅳ in Chl biosynthesis, from protoporphyrin Ⅸ to Chl a. Five chlP genes were identified in radish flesh, and their expression levels were consistently high in GF. ChlP genes encoding geranyl-geranyl reductase provide a hydrophobic alcohol moiety for Chl synthesis [29]. The chlP antisense tobacco plants showed a reduction in Chl content, which resulted in a pale or variegated phenotype [30], and the rice OschlP mutant, 502ys, accumulated lesser Chl compared with that by the wild type [29]. POR catalyzes protochlorophyllide to form the chlorophyllide. There are three POR isoforms (AtPORA, AtPORB, and AtPORC) in Arabidopsis [31]. Two POR isoforms have been identified in barley, tobacco, and rice [31][32][33]. However, only one POR gene has been detected in cucumber and pea [34,35]. Five POR genes with different expression patterns were detected in the transcriptome data of radish flesh. Rsa1.0_08969.1_g00001.1 and Rsa1.0_17789.1_g00001.1, orthologs of AtPORA, showed very low expression levels. Rsa1.0_00630.1_g00021.1 (ortholog of AtPORB), Rsa1.0_00161.1_g00008.1 (ortholog of AtPORC) and Rsa1.0_00556.1_g00016.1 (orthologs of AtPORC) presented higher expression in GF than that in WF flesh (S6 Table). These results indicate that three POR isoforms may exist in radish flesh, and PORB and PORC are essential for Chl formation in radish green flesh. Previous studies also found that POR was closely related to Chl content in chrysanthemum petals and peony leaves [13,36].
Genes such as CAO, chlG, NOL, CLH, and HCAR were expressed in the Chl cycle. CAO and chlG catalyze the formation of Chl b from chlorophyllide a and hydroxychloro-phyllide a, while NOL and CLH are responsible for the reverse progress, i.e., from Chl b to hydroxychlorophyllide a [4]. The high expression of CAO and chlG and the low expression of NOL and CLH in GF suggest that Chl b synthesis was more active than its degradation in GF. This is one of the possible reasons for the accumulation of Chl b in radish green flesh. In addition, HCAR encoding hydroxymethyl Chl a reductase was not detected in the current study. This result differed from previous studies, which showed that a high expression level of HCAR was detected in plant green tissues [13,15]. This difference may be due to the specific tissue (radish green flesh) used in our study. Furthermore, only three DEGs involved in Chl degradation were identified in this study, and the expression of these genes was up-regulated in GF at only one stage, which may indicate that Chl degradation was not active in GF at most stages.
In GF, the contents of Chl a, Chl b and total Chl decreased at S2 compared to those in other stages (Fig 3). A previous study on the cyanobacterium Synechocystis [37] reported that the amount of Chl a decreased with increasing growth rate. We found that the taproot relative growth rate of GF between S2 and S1 was 343.31%, which was much faster than that between other stages (S4 Fig). This decrease in S2 may therefore be caused by the rapid growth of GF taproots from S1 to S2. In addition, we examined the expression trends of the 34 DEGs involved in Chl metabolism in GF across five stages (especially at S2), and found that the expression levels of these genes did not significantly decrease at S2 compared to those in other stages (Fig 7). We further analyzed the correlation between the expression levels of these 34 DEGs and Chl content only in GF. Unlike in both GF and WF, no obvious positive correlation was apparent in GF alone (S9 Table). Previous studies have shown that leaf Chl content is associated with post-transcriptional regulation [38,39], which may also have an effect on the changes in Chl content at different developmental stages in GF. Future experiments are required to explore this possibility.
Leaf photosynthesis is a well-established and well-researched process [40]. However, leaves are not the only part of the plant where photosynthesis can occur. Photosynthesis also occurs in many non-foliar organs, such as in cucumber fruit, Arabidopsis silique walls, and cotton stems, bracts and capsule walls [18,41,42]. The functions of the electron transport chain have been detected even in the pith of the tree stem [17]. In our study, numerous genes associated with photosynthesis were identified, and the vast majority of these genes showed high expression levels in GF (S6 Table). For example, four RBCS genes were highly expressed in GF, however three of them showed low expression in WF. RBCS genes encoding the small subunit of RUBISCO are regarded as key markers of photosynthesis. ALDO and SBPase are two essential enzymes involved in the Calvin cycle. Previous studies have shown that over expression of SBPase and/or ALDO in transgenic plants can result in enhanced photosynthesis [43][44][45][46]. We found six ALDO genes and one SBPase gene that were differentially expressed between GF and WF. Without exception, all these genes exhibited high expression in the GF. Additionally, many key genes participating in photosynthetic electron transport (e.g., petC, petH, psaH, psaK, atpH, atpG, Lhca, Lhcb) and photorespiration (e.g., PGP, HAO, gcvH, glyA, glnA, purU) also showed high expression in GF, although the expression levels of Rsa1.0_05506.1_g00002.1 (psbR) and Rsa1.0_00125.1_g00015.1 (psaH) decreased in GF (Fig 9), probably because of their redundant functions. Chl auto-fluorescence was observed, and high values of FV/FM, FPSII, and ETR were also detected in the GF. These results strongly suggest that photosynthesis occurs in the green flesh of GF taproots. Furthermore, while the radish is a typical C 3 plant [47], some DEGs were arranged in the C 4 pathway according to the KEGG analysis (Fig 8 and  S7 Table). Genes involved in the C 4 pathway are also expressed in C 3 plants, which may perform different functions in C 4 and C 3 plants [48,49]. Further, no difference was detected in the expression levels of PEPC and PPDK (key genes of the C 4 pathway) between GF and WF tissues in our transcriptome data (Fig 8 and S7 Table). Therefore, it is unlikely that the radish green flesh could perform the C 4 photosynthetic pathway.
In this study, the circadian rhythm pathway was enriched by KEGG pathway enrichment analysis. The circadian clock is considered a developmental manager in plants [50]. Many studies have shown that circadian rhythms are closely integrated with photosynthesis [51]. Some genes involved in circadian rhythm (such as CRY, PHY, and HY5) also play important roles in photosynthesis. CRY and PHY are blue-light and red/far-red-light receptors, respectively [52]. Chaves et al. [53] reported that CRY2 may serve to enhance blue-light sensing under conditions where light is restricted, although CRY2 has overlapping functions with CRY1. We observed that CRY2 showed high expression level in GF during the entire period of taproot development, while PHYA and PHYB genes exhibited barely detectable expression levels in GF. The environmental light must pass through the outer tissue to reach the inner part, the process of which reduces the light quantity and quality [17]. Sun et al. [54] reported that the inner tissues of the stem and root were exposed to an internal light environment enriched in far-red light. However, the expression levels of CRY and PHY might suggest that the radish green flesh may use weak blue-light, rather than far-red light, to activate photosynthesis, however further studies are needed to confirm this. HY5 is a bZIP-type transcription factor [55], and it can positively regulate Chl synthesis and chloroplast biogenesis [56,57]. Therefore, our study also focused on the expression of HY5. Six HY5 genes with differential expression levels between GF and WF were identified in the current study. Compared with WF, all HY5 genes showed down-regulated expression from S2 to S5 in GF (S6 Table). The decreased expression of HY5 genes contradicted the increased concentration of Chl in GF. This could be caused by the highly diverse functions of HY5 [58,59].

Conclusions
In this study, RNA-Seq technology was used to compare the gene expression patterns between the flesh of green-fleshed GF and white-fleshed WF during the different stages of taproot development. Through a combined analysis of transcriptome data and Chl content in GF and WF, we can infer that Chl biosynthesis is likely a dynamic, regulated process in GF. Compared to that in the WF, Chl accumulation in GF may be due to the higher expression of Chl biosynthesis genes, especially those associated with Part Ⅳ, and the lower expression of Chl degradation genes, which was also supported by the expression profiling of five candidate genes in diverse cultivars. Additionally, a large number of photosynthesis-related genes were highly expressed in GF, where Chl auto-fluorescence was evident and photosynthetic activity was also observed. These results strongly suggest that the green flesh of GF can perform photosynthesis. Further studies are needed to identify the key regulated genes involved in Chl biosynthesis and to clarify how photosynthesis is initiated and then progresses in the green flesh of radish taproot.    Table. Correlation coefficients between expression levels of 34 DEGs involved in Chl metabolism and contents of Chl a, Chl b and total Chl only in GF. (XLSX) S10 Table. The raw data of photochemical efficiency and ETR for Fig 2. (XLSX) S11 Table. The raw data of Chl content for Fig 3. (XLSX) S12