Metabolome and transcriptome profiling reveal regulatory network and mechanism of flavonoid biosynthesis during color formation of Dioscorea cirrhosa L.

Plant material collection

A total of four tuber types with distinct colors were obtained from their natural habitat located in the hilly region of Shunde, Guangdong Province, China (112.65°E, 22.88°N). All samples were collected from the same location. The color gradient of the four tubers was easily visible (Fig. 1A). Therefore, four differently colored tubers with three independent biological replicates (12 samples in total) were used in this study. The tuber flesh was flash frozen in liquid nitrogen and kept at −80 °C for subsequent experiments.

Determination of total flavonoid content

The detection of total flavonoid content was carried out using the colorimetric aluminum nitrate method (Moreno et al., 2000). The tubers were dried and powdered, then 1.0 g of each sample was accurately taken, mixed with 30 ml of 60% ethanol, extracted ultrasonically for 40 min, followed by suction filtration and centrifugation. A total of 0.1 ml of extract from each sample was pipetted into a burette, then 0.5 ml of 5% NaNO₂ solution, 0.5 ml of 10% Al (NO₃)₃ solution, and 4 ml of 4% NaOH solution were then respectively added. The volume was adjusted with 60% ethanol to 10 ml, shaken, and then the absorbance value was measured at 510 nm. The total flavonoid content was calculated using rutin as the standard (Chen, Wang & Wan, 2010). A significance analysis of the total flavonoid content between groups was performed using the t. test (version 4.0.3) package in R.

Metabolite identification and quantification

Widely-targeted metabolite profiling analyses were performed by Metware Biotechnology Co., Ltd. (Wuhan, China). The identification and quantification of the metabolites were performed as described previously, using the LC-ESI-MS/MS system (UPLC, Ultra Performance Liquid Chromatography; Shimpack UFLC SHIMADZU CBM30A; MS/MS, Applied Biosystems 4500QTRAP) (Dong et al., 2014). Using the MWDB database (Metware Database, http://www.metware.cn), the mass spectrometry data of all samples were obtained and processed using Analyst (software version 1.6.3).

Metabolomic data analysis

A principal component analysis (PCA) was conducted to understand the overall differences between samples in each group as well as the variability between samples (Chen et al., 2009; Eriksson et al., 2006). The hierarchical cluster analysis (HCA) method was used to show the accumulation pattern of different metabolites in the four groups, and the R package heatmap was used to draw the cluster heatmap. The OPLS-DA method was used to maximize the metabolomic differences between the pairs of samples (Thévenot et al., 2015), and the variable importance in projection (VIP) value was used to calculate the relative importance of each metabolite. The PCA and OPLS-DA were both performed using R (software version 4.0.3; R Core Team, 2020). Metabolites with a fold change ≥ 2 and a fold change ≤ 0.5 were selected, and metabolites with a VIP > 1 were considered differential metabolites for group discrimination (Mizoi, Shinozaki & Yamaguchi-Shinozaki, 2012). The K-means analyses were performed in R.

RNA-sequencing, data processing, and annotation

The total RNA isolation was conducted, as previously described, with three biological replicates prepared for each tuber type. The mRNA library construction and sequencing of each sample were performed in the Illumina Hiseq platform using the paired-end reads method. Then, clean reads were obtained by processing the raw reads and removing the low-quality bases, undetermined nucleotides, and adapter sequences. Trinity (software version 2.6.6) (Grabherr, Haas & Yassour, 2011) was used to assemble clean reads, and the longest clusters were obtained using Corset hierarchical clustering (version 1.07) (Davidson & Oshlack, 2014), with soft clusters regarded as unigenes. The gene saturation was calculated and the gene expression level was measured in each sample using reads per kilobase per million (RPKM) value (Wang et al., 2012a; Wang et al., 2012b; Pertea et al., 2015).

To further investigate the functions of unigenes, the assembled unigenes were annotated by the following databases: non-redundant NCBI protein database (NR) (ftp://ftp.ncbi.nih.gov/blast/db) (Deng et al., 2006), Clusters of Orthologous Groups of proteins (KOG) (http://www.ncbi.nlm.nih.gov/COG) (Tatusov et al., 2001), a manually annotated and reviewed protein sequence database (SwissProt) (http://www.gpmaw.com/html/swiss-prot.html) (Boeckmann et al., 2003), Gene Ontology (GO) (http://geneontology.org/) (Ashburner et al., 2000), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp) (Minoru et al., 2008). BLAST (software version 2.10) (Altschul et al., 1990; Camacho et al., 2009) was used to align the sequences to the database and extract the pathways.

Differentially expressed genes (DEGs) analysis

The DEGs analysis was performed using the DEGseq (version 1.30.1) (Varet et al., 2016) package in R. The Benjamini–Hochberg method was used to perform a multiple hypothesis test correction on the hypothesis test probability p-value <0.05 to obtain the false discovery rate (FDR). Differential gene screening was based on the threshold of —log₂Fold Change— ≥ 1 and FDR < 0.05 (Love, Huber & Anders, 2014). To compare the differences in DEGs among samples, we calculated the correlation of samples and drew the expression heatmap based on the Pearson’s correlation coefficient. A gene counting analysis was implemented by featureCounts (version 1.6) (Liao, Smyth & Shi, 2013). Furthermore, the identified DEGs were subjected to an enrichment analysis using GO and KEGG annotations and a pathway analysis using the clusterProfile (version 1.6) R package.

Co-expression network construction

To investigate the gene modules associated with pigment accumulation and PA synthesis, a weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA (version 1.31) (Langfelder & Horvath, 2008)) package in R. Network construction and module identification were conducted using the topological overlap measure (TOM). Gene modules were calculated based on the default settings with minModuleSize = 30, mergeCutHeight = 0.25, and power = 14. The eigengene value, total connectivity and intramodular connectivity, KME (eigengene-based connectivity), and p-values were calculated for the genes in each module, and the correlation coefficients among modules were used to construct the expression network.

Comparison of tuber coloration

To investigate pigment accumulation in D. cirrhosa, we collected tubers at four color stages: LR (light red), RD (red), DR (dark red), and BR (brownish-red). Obvious color differences could be observed in these four tubers (Fig. 1A). Flavonoids are the key substances that affect the phenotypic color changes of plants (Wang et al., 2019a; Wang et al., 2019b). In this work, total flavonoid content was detected in the tubers on the basis of three biological replicates (Table S1). The results showed that the total flavonoid content reached the highest level in DR and then declined to the lowest level in BR. The t-test analysis revealed that 5 out of the 6 groups had significant differences in flavonoid content (Fig. 1B). The results showed that the content of flavonoids was significantly different in the pairwise comparison, indicating those samples could be used for further analysis.

Identification of differentially accumulated metabolites (DAMs) in the tubers

A total of 531 metabolites were detected based on the UPLC MS/MS detection platform, primarily divided into 16 categories including: amino acids and their derivatives, flavonoids, alkaloids, phenolic acids, lipids, nucleotides and their derivatives, fatty acids, lignans, and steroids (Table S2). Among these, 62 flavonoid metabolites (FMs) were identified in this study: 12 flavonoids, 27 flavonols, five dihydroflavones, six flavanols, three chalcones, two dihydroflavonols, and seven proanthocyanidins (Fig. 2C, Table S3). To compare the overall metabolite differences among the groups and the variability between the intra-group samples, the metabolomic dataset obtained was subjected to a principal component analysis (PCA). The results showed that these groups were clearly separated in score, with 48.9% of the first principal component (PC1) and 25.25% of the second principal component (PC2) (Fig. 2A). The heatmap hierarchical clustering showed that the metabolite contents in LR, RD, DR and BR varied greatly (Fig. 2B). These results indicate that the metabolites of the four groups differ significantly and that the metabolomic data were highly reliable.

The orthogonal partial least squares-discriminant analysis (OPLS-DA) method was used to find different metabolites between groups by calculating the correlation between components and by removing irrelevant difference information (Boccard & Rutledge, 2013). A volcano map was used to show the OPLS-DA results (Fig. S1). Based on the threshold of VIP (variable importance) >1 and fold change ≥ 2 or ≤ 0.5, a total of 467 DAMs were identified in the comparison groups. Among these, 234 DAMs were detected between LR and RD (124 up- and 110 down-regulated), 205 DAMs were detected between RD and DR (136 up- and 69 down-regulated), and 390 DAMs were detected between DR and BR (87 up- and 303 down-regulated), respectively (Fig. 2D, Table S4). The results indicate that the DR-BR period may be a crucial turning point for metabolite down-regulation. Conversely, the RD-DR period may be a pivotal point for metabolite up-regulation. Additionally, a KEGG analysis was also performed to explore the functional classification of the DAMs obtained from each comparison group. Except for the RD vs. DR group, the other five groups were found to possess the “flavonoid biosynthesis” KEGG pathway. “Phenylpropane biosynthesis” was enriched in the LR vs. RD, DR vs. DR, LR vs. DR, and RD vs.. BR groups, and “flavone and flavonol biosynthesis” was enriched in DR vs. BR, LR vs. BR, and RD vs. BR (Table S5, Fig. S2). These results indicate that flavonoid biosynthesis may play an important role in the color formation of tubers.

Dynamic trends and patterns of DAMs

By using a K-means analysis, 467 DAMs were divided into nine clusters with similar dynamic and change patterns (Fig. 3, Table S6). Among these, class 4 and class 9 harbored the greatest number of metabolites. The class 4 metabolites gradually increased from LR to DR, then decreased in BR, which was consistent with the change of total flavonoid content, primarily including phenolic acids, amino acids and derivatives, phenolamine, organic acids, saccharides, and flavonoids. Noticeably, only five flavonoids (pinocembroside, 3′-O-Methyl-(-)-epicatechin, gallocatechin, epicatechin-epiafzelechin, and persicoside) were assigned to class 4; these flavonoids may play an important role in tuber color formation. The metabolites in class 9 significantly accumulated in BR, including 12 flavonoids, 10 amino acids, and eight phenolic acids. These results showed that the flavonoids distributed in class 9 may be important metabolites that affect the color formation of BR tubers whereas metabolites in class 4 may affect the color formation of LR, RD and BR tubers.

FM difference in the tubers during color formation

To compare the FM changes in tubers, a heatmap was used to display the varying levels of these compounds. A considerable portion of FMs increased from LR to DR and occurred at high levels in DR, whereas they declined to low levels in BR, including phloretin-4′-O-glucoside, aromadendrin-7-O-glucoside, naringenin-7-O-glucoside, and quercetin-3-O-galactoside (Fig. 4). Only seven FMs were found in BR tubers: naringenin, tamarixetin, luteolin, kaempferol, quercetin-3,4′-dimethyl ether, quercetin, and azaleatin (Fig. 4, Table S3). These upstream substances all exhibited the highest levels in BR compared to the other three groups. In contrast, the downstream metabolites, primarily luteolin-4′-O-glucoside, luteolin-3′-O-glucoside, naringenin-7-O-glucoside, kaempferol-3-O-glucoside, kaempferol-3,7-O-dirhamnoside, quercetin-3-O-glucoside, and quercetin 3-galactoside, exhibited the lowest levels in BR (Fig. 4). These results indicate that BR has more upstream substrates to produce downstream flavonoids, and that these FMs may be associated with the brownish color formation in BR. This may be because the synthesis pathway of downstream metabolites in BR was blocked, resulting in the accumulation of upstream metabolites.

PAs are downstream metabolites of the flavonoid biosynthesis pathway (Taruscio, Barney & Exon, 2004). The PAs detected in this study included procyanidin A6, procyanidin B1, procyanidin B2, procyanidin B3, procyanidin B4, procyanidin C1, and procyanidin C2, catechin, epicatechin, and their derivatives (Fig. 4). It is noteworthy that both epicatechin and procyanidin B2 showed high levels in all groups, while procyanidin B3 only exhibited high levels in BR (Fig. 4). In addition, eight tannins were detected in total. The results showed that 1-O-Galloyl-D-glucose and 2-O-Galloyl-D-glucose were expressed at high levels in RD and DR, while arecatannin B1 showed the highest levels in BR (Fig. 4, Table S3). These findings suggest that these PAs may the key metabolites leading to color formation in tubers, and the tannins primarily accumulated in the later stages of tuber color formation.

Transcriptome analysis

A total of 12 cDNA libraries were generated, and 33.3 Gb of raw bases were produced. High-quality clean reads were generated by removing low-quality reads. The Q30 content of the reads was more than 90%, and the GC content was 45%. Further de novo assembly yielded a total of 131,334 unigenes for analysis, with a mean length of 1,172 bp and an N50 length of 1,863 bp (Table S7). The annotation results of seven public databases showed that 88,183 (67.14%) unigenes were annotated in NR, 70,409 (53.61%) were annotated in GO, 66,135 (50.36%) were annotated in KEGG, 55,987 (42.63%) were annotated in KOG, and 61,518 (46.84%) were annotated in SwissProt. These results indicate that the RNA-seq assembled in this study were reliable for further analysis. Prior to this study, little genomic information was available for D. cirrhosa. This work provides the genomic data for D. cirrhosa.

Functional annotation and enrichment analysis of DEGs

A total of 22,865 DEGs were identified in this study. Among these, 3,195 down-regulated DEGs and 4,487 up-regulated DEGs were identified in LR vs RD, 5,574 down-regulated DEGs and 4,698 up-regulated DEGs were identified in RD vs. DR, and 4,026 down-regulated DEGs and 8,360 up-regulated DEGs were identified in DR vs. BR (Fig. 5A, Table S8). The number of up-regulated genes in DR vs. BR was significantly higher than down-regulated genes, whereas RD vs DR harbored more down-regulated genes than up-regulated genes (Fig. 5B). Thus, we concluded that the RD-DR and DR-BR periods were two important stages for pigment accumulation and color change.

Furthermore, a KEGG analysis was carried out to map the functions of the unigenes between groups. The results showed that plant hormone signal transduction, and flavonoid and isoflavonoid biosynthesis were the dominant biological process categories in the LR vs. RD period. Plant hormone signal transduction, isoflavonoid biosynthesis, and flavonoid biosynthesis were the top categories in RD vs. DR. In DR vs. BR, spliceosome, plant hormone signal transduction, mRNA surveillance pathway, ribosome biogenesis in eukaryotes, and flavonoid biosynthesis were the dominant categories (Fig. 5C). Studies have shown that hormone- and transcription-related genes both perform functions in flavonoid accumulation (Rai et al., 2016; Zheng et al., 2021a; Zheng et al., 2021b). In this work, we found that plant hormone signal transduction and flavonoid biosynthesis may play important roles in the color formation of D. cirrhosa tubers, and there may be a close relationship between these two biological processes in regulating PA synthesis in D. cirrhosa.

Flavonoid biosynthesis pathways and related DEGs between the four groups

A total of 67 key candidate unigenes were found to be involved in PA biosynthesis by searching the KEGG database, including 13 enzymes that directly influence the flavonoid biosynthesis pathways (Table 1, Table S9). These regulated genes might be crucial for the synthesis of flavonoids in D. cirrhosa.

Phenylalanine ammonia-lyase gene (PAL) is the first enzyme that catalyzes the conversion of phenylalanine to cinnamic acid (Ferrer et al., 2008). CHS is effective in flower color regulation (Chen et al., 2012). Flavonol synthase (FLS) genes control flavonol synthesis (Holton, Brugliera & Tanaka, 2010). CHI and FLS can catalyze dihydroflavonol to produce kaempferide, luteolin, and quercetin (Winkel-Shirley, 2001). In this study, the heatmap was conducted to show the gene expression levels among tubers (Fig. 6). Results showed that some PA-related genes only exhibited significantly higher expression in BR, including Cluster 6992.42289 (CHS), Cluster 6992.51311 (FLS), and Cluster 6992.42986 (CHI). This may explain why some metabolites (e.g., luteolin, naringenin, and quercetin) were only found in BR and were not detected in other tubers (Fig. 4). Notably, two PAL genes (Cluster 6992.41423 and Cluster 6992.47474) and one CHS gene (Cluster 6992.41106) showed high expression levels in LR, RD, and DR, but low expression levels in BR, which was consistent with the changes in flavonoid metabolites (luteolin-3-O-glucoside, luteolin-4-O-glucoside, aromadendrin-7-O-rutinoside, hesperetin-7-O-glucoside, and tamarixetin) (Fig. 4). ANS and ANR genes are the key regulators that catalyze cyanidin and dihydromyricetin to PAs, respectively, and DFR is the key enzyme for the catalysis of flavanol to leucoanthocyanidin (Pang et al., 2007; Xie, Sharma & Paiva, 2003). In this work, we found that three PA genes, including the DFR gene (Cluster 6992.39010), ANS gene (Cluster 6992.40887), and ANR gene (Cluster 6992.50236), showed consistently high expression patterns. This is consistent with the high level of PA metabolites found in tubers (epicatechin and procyanidin B2), indicating that these are the key genes for PA synthesis. In contrast, the expression levels of other PA genes showed no significance between groups. These results indicate that these genes are important in the synthesis of flavonoids and PA metabolites and that these may function together to regulate metabolite biosynthesis in D. cirrhosa.

Co-expression network identified PA-related DEGs

WGCNA was used to reveal the correlation patterns of genes associated with PA biosynthesis. A co-expression network analysis was conducted based on the reads per kilobase per million (RPKM) value of 7,963 DEGs, resulting in 15 distinct modules, labeled with different colors, in which genes in the same modules had high correlation coefficients (Fig. 7A, Table S10). Furthermore, the total flavonoid content and the key genes related to PA biosynthesis were both used as trait data for the module-trait relationship analysis (Fig. 7B). Among the modules, MEblue had the highest correlation with PA biosynthesis and was related to ANS_Cluster-6992.41106 by the highest level of r = 0.81, p = 0.001. In contrast, MEblack had the lowest correlation with ANS_Cluster-6992.40887 by the lowest level of r = −0.81, p = 0.001. This indicated that a strong positive regulation and a negative regulation of PA synthesis genes were distributed in the MEblue and MEblack modules, respectively. Through calculating the eigengenes in these two modules, significantly high expression levels were exhibited in DR and BR (Fig. 7D). Further KEGG mapping of the modules with high correlation coefficients showed that MEblue harbored the most KEGG terms (Fig. 7C). It is worth noting that the MEblue and MEblack modules were both highly related to the biosynthesis of phenylpropanoids, flavonoid, amino acids, peroxisome, and the glucagon signaling pathway, while plant hormone signal transduction exhibited particularly high enrichment levels in each module, suggesting potential functions in the color formation process of D. cirrhosa tubers.

Identification of hub TFs involved in pigment accumulation

Aux/IAA transcription factors and AP2/ERF-ERF transcription factors belong to plant hormone gene families, which regulate auxin and ethylene biosynthesis, respectively. They are involved in various biological functions, such as plant development, fruit and seed maturation, and disease resistance (Carole et al., 2013; Ying et al., 2010; Mizoi, Shinozaki & Yamaguchi-Shinozaki, 2012). To further study the correlation between TFs and flavonoid pathway genes in D. cirrhosa, we examined the expression patterns of the 37 candidate TFs with high expression levels in the samples. Among these, 11 were plant hormone pathway TFs (Fig. 8A, Table S11). Results showed that almost all TFs were significantly expressed in DR and BR. In comparation, gene expression levels in LR and RD were low, indicating that flavonoid and hormone biosynthesis were greatly enhanced in the late stages of tuber color formation (DR and BR). Pearson correlations were calculated between TFs and flavonoid genes. We used Cytoscape (version 3.8.2) (Shannon et al., 2003) to draw the gene connection network between 24 transcription factors and 18 flavonoid genes with high correlation (Fig. 8B). Among the TFs, Cluster-6992.36828 (AP2 TF) harbored the most edges related to the flavonoid pathway and was highly expressed in BR. Cluster-6992.45020 (AP2 TF) and Cluster-6992.48662 (MYB TF) harbored the second-most edges and showed high expression levels in the LR and BR tubers, respectively. The Ethylene RESPONSE FACTOR (ERF) and Aux/IAA transcription factor family genes harbored a number of edges with flavonoid genes, including Cluster-6992.45555 (AUX), Cluster-6992.53098 (bZIP), Cluster-6992.43923 (AP2), Cluster-6992.51615 (bZIP), and Cluster-6992.43923 (AP2). Notably, the DFR and ANR genes were also found to be highly associated with TFs, suggesting they may be key regulators in PA biosynthesis (Fig. 8B). These results indicate that the flavonoid TFs and hormone TFs may contribute to flavonoid/PA biosynthesis in D. cirrhosa.