Comparative RNA-Seq analysis to understand anthocyanin biosynthesis and regulations in Curcuma alismatifolia

Rhizomes from C. alismatifolia cv. ‘Dutch Red’ and ‘Chocolate’ were obtained from a horticultural company (Jinluan) in China. Rhizomes were grown in a greenhouse within the campus of Minnan Normal University. During the experiment, the plants were grown at 25 °C under 500 μmol · m⁻² · s⁻¹. Colourful bracts from ‘Dutch Red’ were harvested at three continuous developmental stages [early flowering small flowers (SF), half-flowering (HF) and blossomed flowering (BF)]. The yellow bracts from ‘Chocolate’ were harvested at the peak flowering stage yellow flowers (YF). For each sample, three biological replicates (3 g each) were collected. These tissue samples were immediately subjected to RNA extraction.

High-quality total RNA from the tissue samples was isolated using the RNeasy mini Kit (Qiagen, Germany), following the manufacturer's instructions. gDNAase (provided by the Kit) was added to remove the genomic DNA. The purity and concentration of isolated RNA were checked by a 2100 Bioanalyser (Agilent Technologies, Inc., Santa Clara, CA, USA) and quantified using an ND-2000 (NanoDrop Thermo Scientific, Wilmington, DE, USA). Only high-quality RNA samples (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 8.0, 28S:18S ≥ 1.0, >2 μg) were used to construct the sequencing library.

RNA purification, reverse transcription, library construction and sequencing were performed at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) according to the manufacturer's instructions (Illumina, San Diego, CA, USA). The RNA-Seq transcriptome libraries were prepared using an Illumina TruSeqTM RNA sample preparation Kit (San Diego, CA, USA). Poly(A) mRNA was purified from total RNA using oligo-dT-attached magnetic beads and then fragmented by fragmentation buffer. Taking these short fragments as templates, double-stranded cDNA was synthesised using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA, USA) with random hexamer primers (Illumina). Then the synthesised cDNA was subjected to end-repair, phosphorylation and ‘A’ base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 200–300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (New England Biolabs, Boston, MA, USA) for 15 PCR cycles. After quantifying by TBS380, two RNA-Seq libraries were sequenced in single lane on an Illumina Hiseq NovaSeq 6000 sequencer (Illumina, San Diego, CA, USA) for 2 × 150 bp paired-end reads.

The raw reads obtained from HiSeq2000 sequencing were filtered by fastp software (v0.19.6): reads containing adaptors, reads with more than 5% unknown nucleotides and low-quality reads with more than 20% of bases with a quality value ≤20 were deleted (Wang, 2007). Then all clean reads were assembled with Trinity (release 2013-02-25) with the following settings: inchworm: -K –L 25; Chrysalis: -min_glue 2 -glue_factor 0.05 -min_iso_ratio 0.05 -kk 48 –strand -report_welds -max_mem_reads 1,000,000 (Grabherr et al., 2011). We then filtered out any transcripts with less than 1% of the per-component expression level using a script bundled with Trinity. ORFs were extracted using the Perl script bundled with Trinity (transcripts_to_best_scoring_ORFs.pl). To estimate the number of full length transcripts that had been assembled in our data sets we used the tool benchmarking universal single-copy orthologs (BUSCO) (v5.2.2), which is based on evolutionarily informed expectations of gene content, with default settings. HMMER is then run to assign a score to the candidate amino acid sequence (Seppey et al., 2019). Finally, Trinity connects the transcripts and obtains sequences that cannot be extended on either end. Such sequences are defined as unigenes. When multiple samples from the same species are sequenced, the unigenes from each sample's assembly can be further analysed through sequence splicing and the removal of redundancy using sequence clustering software to acquire non-redundant unigenes with the longest length possible. In the final step, transcripts and unigenes were separately searched against the National Center for Biotechnology Information (NCBI) protein non-redundant (NR), the Clusters of Orthologous Groups of Proteins (COG) and Swiss-Prot databases using DIAMOND (v0.8.37.99) to perform function annotations with a significant threshold E-value of 10⁻⁵. BLAST2GO (v2.5.0) (http://www.blast2go.com/b2ghome) program was used to get GO annotations of unigenes for describing biological processes, molecular functions and cellular components (Conesa et al., 2005). Metabolic pathway analysis was performed using KOBAS (v2.1.1) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Ogata et al., 1999).

To identify DEGs between two different samples, the gene expression level was calculated using fragments per kilobase of exon per million mapped reads (FPKM). RNA-Seq by expectation-maximization (RSEM) (http://deweylab.biostat.wisc.edu/rsem/) (Dewey and Li, 2011) was used to quantify gene and isoform abundances. R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R, http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) (Robinson et al., 2010) was utilised for DEGs between two samples (Anders et al., 2013). The threshold for the p-value was determined by the false-discovery rate (FDR). Unigenes with FDR ≤0.001 and ratio of FPKMs of the two samples larger than 2 (genes for which FPKM <1 were filtered) were considered to be DEGs in this study. In addition, functional-enrichment analysis including GO and KEGG were performed to identify those DEGs that were significantly enriched in GO terms and metabolic pathways at Bonferroni-corrected p-value ≤0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) (Xie et al., 2011).

Primers used for qRT-PCR assays were designed with the online Primer3 (v0.4.0) software and they amplified PCR products that varied from 80 bp to 200 bp. The housekeeping gene ACTIN (Li et al., 2021b) was used as internal reference control for normalisation. qRTPCR assays were performed as described previously (Li et al., 2019). The relative expression of the genes was calculated using the 2^−ΔCt method.

To obtain a more detailed understanding of the anthocyanin regulatory network within different varieties and different flowering developments, four samples, derived from the red bracts from the SF, HF and BF stages of ‘Dutch Red’, and the yellow bracts from BF stage of ‘Chocolate’ (YF), were used as materials for RNA sequencing (Figure 1A). Each sample was performed with three biological replicates; thus a total of 12 sequenced libraries were constructed.

In total, approximately 576.45 Mb of raw read data were generated from Illumina HiSeq sequencing ranging from 45.24 Mb to 52.98 Mb in each sample (Table S1 in Supplementary Materials). After raw read filtering, a total of 566,269,472 clean reads (the number of total clean nucleotides was 83,999,856,857 nt; the Q20 and GC percentages were 98.31% and 49.25%, respectively) were obtained from the 12th libraries (SF-1, SF-2, SF-3; HF-1, HF-2, HF-3; BF-1, BF-2, BF-3; YF-1, YF-2, YF-3). Then, clean reads were assembled into 159,687 transcripts and 69,453 unigenes, reaching an averagely mapped ratio of 73% in samples. The assembled transcript length ranged from 500 bp to over 4,500 bp with an average of 1,236 bp and N50 of 1,830 bp (Table 1). Other statistical results obtained from the C. alismatifolia transcriptome sequencing and assemblies are also listed in Table 1. The accession numbers of all raw data in the Short Read Archive (SRA) Sequence Database in the NCBI are listed in Table S1 in Supplementary Materials.

These transcripts offer a potential source for the identification of functional genes. Subsequently, reads with an FPKM value <0.5 were removed, and 95,625 (mean from three replicates, the same below), 88,479, 89,429 and 78,537 transcripts were found in SF, HF, BF and YF, respectively. On average, 63% of the transcripts were in the 0.5–5 FPKM range, and 35% of the transcripts were in the 5–100 FPKM range (Figure 1B; Table S2 in Supplementary Materials). A principal component analysis showed that there were highly correlated transcriptome characteristics between the biological replicates of each sample (Figure 1C). A total of 44,536 transcripts were shared between samples (Figure 1D).

The transcript and unigene sequences were first aligned against the NR and Swiss-Prot databases using BlastN and BlastX searches with an E-value less than 10⁻⁵. This analysis indicated that 104,303 transcripts (65.32% of the total transcripts) and 34,189 unigenes (49.23% of the total unigenes) were matched to the NR database (Table S2 in Supplementary Materials), while 82,193 transcripts (51.47%) and 26,589 unigenes (38.28%) were matched to the Swiss-Prot database (Table S3 in Supplementary Materials). The annotated transcripts and unigenes in the COG, GO and KEGG databases are also listed in Table S3 in Supplementary Materials. However, 53,694 transcripts (33.62%) and 34,270 unigenes (49.34%) could not be matched to any functions included in these five databases.

Comparison of gene expression levels in varieties and different development stages revealed a total of 15,449 DEGs (BF vs. YF, 11,000 DEGs; BF vs. HF, 4,761 DEGs; BF vs. SF, 7,742 DEGs; HF vs. SF, 1,535 DEGs), and only 1% (215 DEGs) of the DEGs between varieties and different development stages are the same (Figure 2A). Interestingly, the number of down-regulated DEGs was greater than the number of up-regulated DEGs in red bracts of blooming flowers of ‘Dutch Red’, as indicated from the comparisons of BF vs. HF (1,643 DEGs up-regulated; 3,118 DEGs down-regulated) and BF vs. SF (3,591 DEGs up-regulated; 4,151 DEGs down-regulated); however, more genes were significantly up-regulated in the bracts of ‘half flowering’ of ‘Dutch Red’, as indicated from the comparisons of HF vs. SF (1,198 DEGs up-regulated; 337 DEGs down-regulated) (Figure 2A).

Figure 2

KEGG analysis of the identified DEGs showed most enrichment in the following pathways: ‘Phenylpropanoid biosynthesis’ (map00940), ‘Plant-pathogen interaction’ (map04626), ‘Plant hormone signal transduction’ (map04075), ‘Flavonoid biosynthesis’ (map00941), ‘Starch and sucrose metabolism’ (map00500), ‘Flavone and flavonol biosynthesis’ (map00944) and ‘Terpenoid backbone biosynthesis’ (map00900) (Figure 2B). We noted that among the DEGs associated with flavonoid, those involved in anthocyanin biosynthesis, and transcripts of shikimate O-hydroxycinnamoyltransferase (HCT, EC2.3.1.133), the enzyme which is of great importance in plant secondary metabolite production, are more enriched in the initial blossoming stage when compared with the blossomed stage of ‘Dutch Red’. Those proteins, involved in stress resistance and antioxidant defence, were more enriched in the half and fully blossomed stages when compared with the initial blossoming stage. In addition, most of the transcripts that participated in photosynthesis were more enriched in the yellow bracts of ‘Chocolate’ when compared with those in ‘Dutch Red’.

Heat maps of the clustering analysis results also suggested that genes related to the ‘Flavonoid biosynthesis’, ‘Flavone and flavonol biosynthesis’ and ‘Phenylpropanoid biosynthesis’ pathways were significantly up-regulated in the early stage bracts of ‘Dutch Red’ when compared with those in the half and fully blossomed stages (Figure 3A). Analysis of DEGs related to ‘Signal transduction’ revealed that auxin-responsive proteins were screened as the most variable gene family, and auxin-responsive small auxin up RNA (SAUR) genes were important in the early development stage of red bracts. Meanwhile, IAA transcripts were mainly produced in both the half and fully blossomed stages when compared with the initial blossoming stage. Besides, the transcripts of 3-ketoacyl-CoA synthase (XP_009405939.1), a key enzyme in fatty acid biosynthesis, were observed exclusively enriched in an early development stage compared with the half and fully blossomed stages (Figure 3A). From the clustering analysis of blossoming stages between ‘Chocolate’ and ‘Dutch Red’, DEGs related to ‘flavonoid biosynthesis’ and ‘Photosynthesis’ were more up-regulated in the yellow bracts of ‘Chocolate’ (Figure 3B). DEGs in the ‘Porphyrin and chlorophyll metabolism’ pathway were significantly enriched in ‘Chocolate’. DEGs in ‘Phenylproanoid biosynthesis’ showed that peroxidase transcripts were significantly enriched in the red bracts of ‘Dutch Red’, which was suggestive of a stress response in the blossomed bracts of ‘Dutch Red’ (Figure 3B).

Figure 3

By annotation in the public database, more than 30 genes affiliated to nine structural gene families were predicted to participate in the anthocyanin biosynthesis pathway (Figure 4). Anthocyanin biosynthesis started with the conversion of phenylalanine to coumarate-CoA by phenylalanine ammonia lyase (PAL), cinnamate-4-hydroxylase (C4H) and 4-coumarate: CoA ligase (4CL). There were seven unigenes annotated as PAL and 12 unigenes annotated as 4CL, among which DN13924_c2 and DN18026_c0 had a significantly lower expression in the blossomed stage of the ‘Dutch Red’. Subsequently, dihydroflavonol was formatted by coumarate-CoA and malonyl-CoA catalysed by CHS, CHI, flavanone-3-hydroxylase (F3H) and flavonoid 3′,5′-hydroxylase (F3′5′H), which were the key enzymes in the metabolism of anthocyanin. The formation of various anthocyanidins by dihydroflavonols was then catalysed by DFR and ANS. Among them, DN14172_c2, DN16691_c4 and DN17923_ c1 encode F3H, DFR and CHS, respectively (Table 2). The first two were the same genes as that submitted previously in the C. alismatifolia by Kriangphan et al. (2009), while DN17923_c1 was closely related to CHS (AEU17693.1) gene in C. longa (Resmi and Soniya, 2012). They all showed significantly higher expression in the bracts of ‘Chocolate’ (FPKM as 1440.4, 570.8 and 2667.3, respectively) (Figure 4; Table 2). Another three anthocyanin structural genes, DN10857_c0, DN16064_c1 and DN12990_c1, have 70.9%, 84.2% and 84.0%, respectively, nucleic acid similarity with CHI (XP_009384766.1), F3′5′H (XP_009411862.1) and ANR (XP_009412668.1) in Musa acuminate. DN10857_c0 was expressed significantly lower in the blossomed stage of ‘Dutch Red’ (FPKM as 38.0) than those in the bracts of SF, HF and YF (FPKM as 501.8, 344.7 and 561.0, respectively). DN16064_c1 was expressed significantly lower in the bracts of YF than those of ‘Dutch Red’. A total of 14 genes were identified in the clade with glucosyltransferase-like family with different patterns of expression (Figure 4). Among them, DN12119_c3 and DN13592_c1 encode the anthocyanidin/anthocyanin 3-O-glucosyltransferases, which catalyse unstable anthocyanidin into anthocyanin. They are the last key enzymes in the anthocyanin biosynthetic pathway. DN12119_c3, DN16301_c1 and DN13065_c1 had significantly lower expressions in BF than that in SF and HF, while other unigenes showed no significant divergence in different stages or cultivars (Table 2; Figure 4).

Figure 4

In addition to the structural genes in the anthocyanin biosynthetic pathway, TFs also play important roles in flower colour development through regulating the temporal and spatial expression of structural genes. Many DEGs identified in the KEGG database were annotated as TFs, including MYB, AR2/ERF, C2C2, bHLH and WD. Their abundance levels and detailed information are shown in Figure 5A and Table S4 in Supplementary Materials. MYB superfamily and WD-repeat transcripts were the most abundant found to regulate the flower development and flower formation (Figure 5B). DN17014_c0 was annotated as expansin-A15, which belongs to the MYB superfamily, and DN13707_c0 was annotated as WD-repeat–containing protein. They were the two important TFs that expressed highly during the whole flower developments (Table 3). Further comparison showed that DN17014_c0 shares sequence similarity with the TF MYB105-like gene, and the average value was about 95%. DN11279_c3 and DN19744_c1 annotated as MYB-related proteins are highly expressed during initial flower development in ‘Dutch Red’ (FPKM as 105.1 and 160.6, respectively), but during the blossomed stages they were significantly decreased. This suggests that DN11279_c3 and DN19744_c1 were more likely the floral activators of ‘Dutch Red’. Another MYB protein, DN18449_c0, annotated as trihelix TF GT-3a, was exclusively expressed highly in the blossomed bracts of ‘Dutch Red’ (FPKM as 116.7), but was significantly lower in the blossom bracts of ‘Chocolate’ (Table 3). A comparison of the MYB TFs between ‘Dutch Red’ and ‘Chocolate’ showed that DN21970_c4, annotated as divaricata-like MYB gene, had significantly high expression in the blossomed bracts of both ‘Dutch Red’ and ‘Chocolate’ (Table 3), indicating that DN21970_c4 was possibly the common TF in regulating later floral development in C. alismatifolia.

Figure 5

To evaluate the reliability of RNA-Seq analysis using Illumina sequencing, 14 candidate genes in the anthocyanin biosynthesis and regulations were selected for qRT-PCR test in the four bracts (SF, HF, BF and YF). The oligonucleotide primer sequences used to amplify these transcripts for each gene and related unigene are listed in Table S5 in Supplementary Materials. In general, our qRT-PCR results show a high degree of consistency with the RNA-Seq results (Figure 6). Expression patterns of six (42.85%) genes (F3H1, F3H2, MYB4, WD2, MYB5 and DFR) had almost similar expression patterns but with very small partial inconsistencies compared with the RNA-Seq results. Seven (50%) genes (CHS2, CHI, F35H, MYB1, MYB2, MYB3 and WD1) fit well with the RNA-Seq results across all four samples. The gene CHS1 expression was up-regulated significantly in the yellow ‘Chocolate’ bracts, which was inconsistent with our FPKM database. It is rational and acceptable that there are certain differences in direct comparisons between the RNA-Seq and qRT-PCR results due to different normalisation methods and other technical biases.

Figure 6