Transcriptomic dynamics changes related to anthocyanin accumulation in the fleshy roots of carmine radish (Raphanus sativus L.) characterized using RNA-Seq

Plant material and experimental design

Based on the quantification of anthocyanin levels of carmine radish “Hongxin 1” using HPLC analysis in our previous study (Gao et al., 2019b) (Fig. S1), young fleshy roots obtained from the development stage of carmine radish “Hongxin 1” were used for RNA-Seq in this study. The carmine radish “Hongxin 1” was selected and kept in a greenhouse at the experimental base of the Yangtze Normal University experiment in 2019. First, seeds of “Hongxin 1” were sown in sterile soil under normal growth conditions (23 °C, 16 h light/8 h dark) for 2 weeks. Subsequently, vernalization treatment was conducted, with the 2-week-old plants transferred to and maintained in a cold room (5 ± 1 °C, 12 h light/12 h dark) for 15 days. After vernalization treatment, the plants were grown in a growth room under normal growth conditions (23 °C, 16 h light/8 h dark). Finally, fleshy roots obtained from each stage of development (seedling stage (SS), initial expansion (IE), full expansion (FE), bolting stage (BS), initial flowering stage (IFS); full-bloom stage (FBS) and podding stage (PS)) of carmine radish “Hongxin 1” were collected for RNA-Seq, with three independent biological replicates for each stage and two technical replicates. All harvested tissue was immediately frozen in liquid nitrogen and stored at −80 °C for RNA-Seq analysis.

Sample preparation and library construction

Library construction was conducted with NEBNext Ultra RNA Library Prep Kits for Illumina (NEB, Ipswich, MA, USA), as follows: The mRNA was isolated with oligo(dT) magnetic beads from approximately 5 µg of total RNA and subsequently converted into short fragments by fragmentation buffer. Then, short fragments were converted into first-strand cDNAs and used as templates for second-strand synthesis with random hexamers. PCR amplification was conducted using the desired purified synthesized cDNA fragments (QiaQuick PCR kit). Ultimately, 200 bp paired-end reads were generated from the prepared library with 2 replicates using Illumina HiSeqTM 2000.

Read processing and differentially expressed gene (DEG) identification

Using Trimmomatic software, clean reads were obtained by filtering out low-quality reads (base quality ≤10) and adaptor-only reads and then trimming the remaining reads. After that, the small subunit (SSU) and large subunit (LSU) rRNA sequences were downloaded from the Silva database (Christian et al., 2013) and used for high-quality read alignment using BWA software (Li & Richard, 2009), and the mapped rRNA reads were removed by a homemade Perl script. Subsequently, we used Trinity software to de novo assemble the remaining clean reads into transcripts. In addition, BLASTx searches of the Swiss-Prot protein databases and nonredundant protein (Nr) database of the National Center for Biotechnology Information (NCBI) were used to annotate all unigenes to obtain the assembled sequences. The parameters were set as follows: E-value < 1e−10, identity > 70%, query coverage ≥ 80%, and default values for other parameters. After gene annotation, the Gene Ontology (GO) database (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.kegg.jp/) were used to annotate the functions of assembled unigenes. To assess the abundances of assembled transcripts, we first mapped the clean reads obtained from different fleshy root libraries using Bowtie2 on the de novo-assembled transcriptome. Then, quantification of the de novo assembly transcript was assessed with RSEM, and only transcripts with FPKM ≥ 1 were considered significantly expressed transcripts (Dewey & Li, 2011). Finally, differentially expressed genes (DEGs) with a corrected P-value < 0.05 between each set of compared samples were identified through screening by noiseqbio (Tarazona et al., 2012). Furthermore, DEGs related to anthocyanin biosynthesis in carmine radish were analyzed and plotted using neighbor-joining clustering with a homemade R script.

GO functional annotation and KEGG pathway analysis of comodulated differentially expressed genes (DEGs) in the growth stages of carmine radish

Co-modulated DEGs (DEGs common to all the dynamic growth stages of fleshy roots in carmine radish) were identified based on Venn diagram analysis. Subsequently, we conducted GO annotation of comodulated DEGs through the AgriGO website (http://systemsbiology.cau.edu.cn/agriGOv2/) and KEGG pathway enrichment analysis using KOBAS software (Xie et al., 2011). Their respective graphs were constructed using a homemade R scripts.

Validation of candidate DEGs involved in the growth stages of carmine radish using real-time qRT-PCR

To validate the results of the RNA-Seq analysis, 11 DEGs with substantial alterations related to the growth stages of carmine radish were chosen for validation by qRT-PCR. The primers used for qRT-PCR experiments were designed by Real-time PCR (TaqMan) Primer and Probes Design Tool (https://www.genscript.com/tools/real-time-pcr-taqman-primer-design-tool), and the radish ACTIN gene was used as a reference gene (Table S1). Following the standard protocol for the ABI 7500 system, the amplification of the candidate genes was performed using qRT-PCR in triplicate as described by Jian et al. (2015). The fold change in the expression levels of target genes was calculated using the relative quantitative method (2^−ΔΔCT) (Schefe et al., 2006), and variance analysis was conducted to evaluate differences in the relative expression levels of the genes between different samples, followed by multiple comparisons using Duncan’s least significant range (LSR) tests with SPSS statistical software. Different letters indicate significant differences at the p = 0.05 level. Of those, a, b, c, d, e, f, or g indicates a significant difference from a column with no superscript letter in common (P < 0.05); columns marked with the same letter are not significantly different at the 5% level by Duncan’s multiple range test.

DEGs related to the dynamic growth stages of fleshy roots in carmine radish

Here, transcriptome analysis of anthocyanin biosynthesis in the carmine radish “Hongxin 1” was conducted followed by HPLC analysis, and young fleshy roots of “Hongxin 1” at different developmental stages were selected for RNA-Seq. The cDNAs obtained from fleshy roots at seven growth phases (SS, IE, FE, BS, IFS, FBS and PS) were sequenced using Illumina sequencing technology. To harvest clean reads for further assembly, we filtered out low-quality reads (base quality ≤10) and adaptor-only reads and then trimmed the sequences to obtain high-quality reads. Subsequently, we aligned the high-quality reads to the small subunit (SSU) and large subunit (LSU) rRNA sequences and removed rRNA reads by a homemade Perl script. After removing rRNA sequences, the average percentage of clean read counts was 70% of raw tags in each library (Table 1).

To identify the candidate genes related to the dynamic growth stages of fleshy roots in “Hongxin 1” radish, we normalized the expression levels of all the globally expressed genes for further analysis, and the results showed that highly distinct gene expression profiles exist in the dynamic growth stages of fleshy roots (Fig. 1; Table S2). Cluster analysis according to the gene expression dynamics was performed by the k-means method, and differentially expressed genes were divided into 9 clusters associated with distinct expression profiles. More interestingly, we found that the putative candidate genes belonging to Cluster 8 were consistent with the dynamics of the anthocyanin profiles of fleshy roots during the development of carmine radish, while the putative candidate genes categorized into Cluster 9 had the opposite pattern. Furthermore, DEGs related to anthocyanin accumulation in carmine radish (Raphanus sativus L.) fleshy roots were identified between other developmental periods (“IE_root”, “FE_root”, “BS_root”, “IFS_root”, “FBS_root” and “PS_root”) and the “SS_root” period. Compared with “SS_root”, 1,629, 1,037, 1,385, 1,521, 1,574 and 917 DEGs were generated in each different periods of development, including both upregulated DEGs (878, 755, 718, 838, 852 and 555 transcripts) and downregulated DEGs (751, 282, 667, 683, 722 and 362 transcripts) (Fig. 2A). Of those, 126 comodulated DEGs were identified based on a Venn diagram (Fig. 2B), and we used a heatmap to display the expression pattern changes of comodulated DEGs (common DEGs involved in the dynamic growth stages of fleshy roots in carmine radish) with different colors (Fig. 2C). Moreover, we found that some comodulated DEGs showed similar expression trends in the dynamic growth stages of fleshy roots, which were again consistent with anthocyanin profiles of the fleshy roots, such as a series of functional enzymes that acted as important regulators in anthocyanin biosynthesis, including dihydroflavonol 4-reductase (RsDFR: Cluster_13775), flavonoid 3′-monooxygenase (RsF3′H: Cluster_4431), leucoanthocyanidin dioxygenase (RsLDOX: Cluster_3903) and chalcone synthase (RsCHS: Cluster_39833), as well as some regulatory enzymes, such as anthocyanidin 3-O-glucoside 2′″-O-xylosyltransferase (Rs3GGT1: Cluster_9270), coumaroyl-CoA:anthocyanidin 3-O-glucoside-6″-O-coumaroyltransferase 1-like (Rs3AT1: Cluster_46827), UDP-glycosyltransferase 75C1-like (RsUGT75C1: Cluster_2736) and UDP-glycosyltransferase 78D2-like (RsUGT78D2: Cluster_11854). In addition, transport proteins and transcription factors, namely, glutathione S-transferase F12 (RsGSTF12: Cluster_24268), MYB transcription factor (RsMYB2: Cluster_28373), and zinc finger, RING-type protein (RsRZEP: Cluster_7186), were identified (Fig. 2C; Table S3).

Functional annotation of DEGs related to the dynamic growth stages of fleshy roots in carmine radish

To explore the regulatory mechanisms of DEGs related to the dynamic growth stages of fleshy roots in carmine radish, GO annotation and KEGG pathway enrichment of those putative DEGs were conducted. The results indicated that “anthocyanin-containing compound biosynthetic process (GO:0009718)” and “anthocyanin-containing compound metabolic process (GO:0046283)” were commonly overrepresented in the dynamic growth stages of fleshy roots after the initial expansion (i.e., 40 days after planting), “flavonoid biosynthetic process (GO:0009813)” and “flavonoid metabolic process (GO:0009812)” were overrepresented in IFS, FBS and PS, but “pigment biosynthetic process (GO:0046148)” and “pigment metabolic process (GO:0042440)” were overrepresented in IFS and FBS. Moreover, we found that the GO terms “glucosinolate biosynthetic process (GO:0019758)” and “glucosinolate metabolic process (GO:0019757)” were overrepresented in IE (Fig. 3A; Table S4). KEGG pathway enrichment analysis identified five significantly enriched pathways related to anthocyanin synthesis traits in carmine radish, including flavonoid biosynthesis, flavone and flavonol biosynthesis, diterpenoid biosynthesis, anthocyanin biosynthesis, and benzoxazinoid biosynthesis (Fig. 3B; Table S5).

Validation of DEGs related to anthocyanin synthesis traits through qRT-PCR

To confirm the RNA-Seq results, eleven putative DEGs (Cluster_13775, Cluster_4431, Cluster_3903, Cluster_39833, Cluster_9270, Cluster_46827, Cluster_2736, Cluster_11854, Cluster_24268, Cluster_28373 and Cluster_7186) were selected and subjected to qRT-PCR analysis (Fig. 4). The changes in expression of these 11 putative DEGs related to anthocyanin synthesis were measured by log2-fold change, and the Pearson correlation analysis indicated excellent concordance between RNA-Seq and qRT-PCR (R2 = 0.76198) (Table S6; Fig. 5). Here, transcripts of dihydroflavonol 4-reductase (RsDFR: Cluster_13775), flavonoid 3′-monooxygenase (RsF3′H: Cluster_4431), chalcone synthase (RsCHS: Cluster_39833) and leucoanthocyanidin dioxygenase (RsLDOX: Cluster_3903) were identified and validated using qRT-PCR. The transcripts of RsDFR, RsF3′H and RsCHS, but not RsLDOX, were significantly upregulated in the different growth stages of fleshy roots among other different development periods from seedling stage to full-bloom stage but decreased in the podding stage, showing consistency with the dynamics of anthocyanin profiles of fleshy root in development of carmine radish (Figs. 4A, 4B, 4C and 4D); similar patterns were observed for UDP-glycosyltransferase 75C1-like (RsUGT75C1: Cluster_2736) and UDP-glycosyltransferase 78D2-like (RsUGT78D2: Cluster_11854), which have been reported to be involved in the methylation of anthocyanidins, resulting in stable compounds. We found that the expression of RsUGT75C1 was significantly different from that of RsUGT78D2 during the fleshy root growth development stages compared with the seedling stage (Figs. 4E and 4F). In addition, anthocyanidin 3-O-glucoside 2′″-O-xylosyltransferase (Rs3GGT1: Cluster_9270) and coumaroyl-CoA:anthocyanidin 3-O-glucoside-6″-O-coumaroyltransferase 1-like (Rs3AT2: Cluster_9270) might act as key regulatory enzymes for the formation of anthocyanins in carmine radish. We found that Rs3GGT1 was changed significantly compared with Rs3AT2 (Figs. 4G and 4H). Moreover, transport proteins and transcription factors such as glutathione S-transferase F12 (RsGSTF12: Cluster_24268), MYB transcription factor (RsMYB2: Cluster_28373), and zinc finger, RING-type protein (RsRZEP: Cluster_7186) were also identified. The transcripts of RsGSTF12 and RsRZEP were not significantly differentially expressed in the other growth stages of fleshy roots compared with the seedling stage (Figs. 4I and 4K), but we found that the transcript of RsMYB2 (Fig. 4J) was significantly upregulated in the fleshy root growth stage and decreased in podding stage compared with other development periods from seedling stage to full-bloom stage. This pattern was consistent with the dynamics of anthocyanidin profiles across the fleshy root development stages in the carmine radish.