In our previous work, we planted more perilla varieties in light incubator and found that the plant color of ‘PF899’ variety was more sensitive to light intensity. The variety is extremely sensitive to changes in light intensity, with green leaves in weak light and dark red leaves in increased light intensity. Hence, the light-sensitive perilla variety was selected in the study. After germination, the plants were grown in a greenhouse at the Guangzhou University of Chinese Medicine in Guangdong province, China (23°3’N, 113°23’E) under a weak light intensity (WL; 15 μmol/m²) until the three-leaf stage. The plants were transferred to moderate light (ML; 180 μmol/m²) for 6 days, followed by transfer to weak light (15 μmol/m²) for 6 days (recovery light conditions, RL). The plants were grown under a 12 h/12 h light/dark photoperiod and fertilized according to standard practices. Each treatment group was treated with three pots, each pot had five plants. Three biological replicates were performed for each treatment group. The leaves from three treatments were collected, immediately frozen in liquid nitrogen, and stored at –80°C for metabolomic and transcriptomic analysis. And in transcriptome analysis, a group of materials at moderate light (ML; 180 μmol/m2) for 12 h was added.

Leaves from different treatment groups were weighed and pulverized. To measure total anthocyanins, each 0.1 g leaf tissue sample was combined with 0.1% hydrochloric acid methanol reagent and incubated at 40°C for 40 min to extract total anthocyanins. The optical densities at 530, 620, and 650 nm were measured with a UV spectrophotometer; 0.1 mol/L hydrochloric acid methanol solution was used as a blank control. Total anthocyanin content was calculated according to Greey’s formula (Lee and Ohnishi, 2003; Nitta et al., 2005b). To measure total flavonoids, perilla leaf powder was combined with an equal volume of 70% ethanol, and 1 mL of the extract and 1 mL of 5% sodium nitrite solution were added to a volumetric bottle. After shaking well for 6 min, 1 mL of 10% aluminum nitrate test solution was added to the flask, which was then mixed for 6 min. After adding 10 mL of sodium hydroxide test solution, the sample was mixed. Finally, 70% ethanol was added to bring the volume to 25 ml. After 15 min, the absorbance of the sample at 500 nm was measured by UV-vis spectrophotometry with 70% ethanol as a blank control. The rutin was used as the standard curve to calculate the content of total flavonoids.

The freeze-dried samples were crushed, and 100 mg of powdered tissue was combined with 1.2 mL of 70% methanol. After centrifuging at 12,000 rpm for 10 min, the extracts were filtered (SCAA-104, 0.22 μm pore size), and the metabolites in the extracts were analyzed used an UPLC-ESI-MS/MS system (UPLC, SHIMADZU Nexera X2; MS, Applied Biosystems 4500 Q TRAP). The conditions were as follows: UPLC: column, Agilent SB-C18 (1.8 μm, 2.1 mm * 100 mm). The mobile phase consisted of solvent A (pure water with 0.1% formic acid) and solvent B (acetonitrile with 0.1% formic acid). Sample measurements were performed with a gradient program. The column temperature was set to 40°C and the injection volume was 4 μL. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. LIT, and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (Q TRAP), AB4500 Q TRAP UPLC/MS/MS System, equipped with an ESI Turbo Ion-Spray interface, operating in positive and negative ion mode, and controlled by Analyst 1.6.3 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature 550°C; ion spray voltage (IS) 5500 V (positive ion mode)/–4500 V (negative ion mode); ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set at 50, 60, and 25.0 psi, respectively; the collision-activated dissociation (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solution in QQQ and LIT mode, respectively. QQQ scans were acquired as multiple reaction monitoring (MRM) experiments with collision gas (nitrogen) set to moderate. Positive ion mode optimizes cluster removal voltage (DP) and Collision electric (CE) for individual MRM transitions were performed by further DP and CE optimization. A specific set of MRM transitions was monitored for each period based on the metabolites eluted within that period.

A Tiangen RNA Prep Pure Kit for Plants was used to extract total RNA from P. frutescens leaves. A NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States) was used to measure total RNA concentration and quality. RNA samples that met the requirements were sent to Wuhan Bena Technology Service Co., Ltd. (Wuhan, China) for transcriptomic sequencing. RNA-seq libraries (P. frutescens leaves collected at WL, ML12 h, ML6d, and RL × 3 replicates) were constructed and sequenced using the Illumina HiSeq4000 platform, and clean reads were mapped to the perilla reference genome using HISAT with default parameters.

The metabolome data in this study was identified with chemical library MWDB (Metware Database), which were established by Metware Metabolism Company. Materials was quantitative based on secondary spectrum information, and isotope signals were removed during analysis, which contained K + ions and Na + ionsRepeated signals of NH4 + ions and fragments of other substances with larger molecular weight were identified. Metabolites with significantly different levels between groups were determined based on Variable Importance in Projection (VIP) ≥ 1 and absolute log2FC (fold change) ≥ 1. VIP values were extracted from the OPLS-DA results, and score plots and permutation plots were generated using the R package MetaboAnalystR. The data were log transformed (log2) and subjected to mean centering before OPLS-DA. To avoid overfitting, a permutation test (200 permutations) was performed. Gene expression levels were calculated as fragments per kilobase of transcript per million mapped reads (FPKM). Differentially expressed genes (DEGs) were identified from normalized read count data using DESeq2.34. Genes with | log2 (fold change) | ≥ 1 and p < 0.01 were considered to be DEGs.

The HCA (hierarchical cluster analysis) results of samples and metabolites were presented as heatmaps with dendrograms, while PCC (Pearson correlation coefficients) between samples were calculated by the core function in R and were presented only as heatmaps. Both HCA and PCC were carried out with the R package pheatmap. For HCA, normalized signal intensities of metabolites (unit variance scaling) were visualized as a color spectrum. Unsupervised PCA (principal component analysis) was performed using the statistics function prcomp within R¹. The data were unit variance scaled prior to unsupervised PCA.

The identified metabolites were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Compound database². Gene function annotation was performed using four databases: National Center for Biotechnology Information (NCBI³), non-redundant protein sequences (NR⁴), the Swiss-Prot protein sequence database, and Gene Ontology (GO⁵). The annotated metabolites and genes were then mapped to the KEGG Pathway database⁶. Pathways with significantly regulated metabolites were subjected to MSEA (metabolite sets enrichment analysis), and their significance was determined based on p-values from hypergeometric tests. Finally, the DEGs in the WL, ML12h, ML6d, and RL groups were subjected to GO and KEGG pathway enrichment analysis.

The weighted gene co-expression network analysis (WGCNA) package was used to generate co-expression network modules between metabolites and genes. Using the automatic network construction function (blockwise Modules) with default parameters, co-expression modules were obtained based on TOM (topological overlap measure). The initial clusters were merged on eigengenes. The eigengene value was calculated for each module and used to search for associations with key anthocyanin substances. The transcriptional regulatory networks were generated by combining the Pearson correlation coefficient (PCC > 0.85) between genes and transcription factors, and the cis-element binding sites were predicted in the promoter regions of key anthocyanin genes in the same module. The networks were visualized by CYTOSCAPE (v.3.7.2, United States) (Kohl et al., 2011).

The conservative domain of the candidate MYBs and bZIPs from perilla and other plants were used to construct phylogenetic. Accession numbers of MYBs and bZIPs from other plants are shown in Supplementary Table S8. The neighbor-joining (NJ) tree was constructed using the poison model with MEGA-X with 1000 bootstrap replicates. The phylogenetic trees were visualized using the iTOL web tool⁷.

Quantitative reverse-transcription PCR was performed as previously described (Livak and Schmittgen, 2001). Total RNA was isolated from the samples using a MAGEN RNA Extraction Kit (MAGEN, United States), and reverse transcription was performed using an Evo M-WLV RT Kit (Accurate, Hunan, China). Primers for each gene (ANS, CHI, DFR, and so on) were designed using Primer 5.0; the primers are listed in Supplementary Table S9. The specificity of the primers was verified by agarose gel electrophoresis. PCR amplification was performed using a LightCycler 480 II REAL-TIME PCR system (Roche, Basel, Switzerland). The qPCR cycling conditions were 1 cycle of pre-denaturation at 95°C, 30 s; 40 cycles of 95°C for 5 s and 60°C for 30 s; dissolution at 95°C for 5 s; 60°C for 1 min. All genes were amplified with three replicates, and the perilla Actin gene was used as an internal reference gene.

To identify the leaves color and major pigment components under different light, we observed and measured the total anthocyanin and total flavonoid contents during different light intensity treatments. The perilla leaves color under different light intensities were obviously affected. The leaves and seedlings of the light-sensitive red perilla variety were green when cultured under weak light intensity (15 μmol/m², WL). After transferring to moderate light intensity (180 μmol/m², ML) for 6 days, the leaves and whole plants turned dark red. When returned to weak light intensity for 6 days (15 μmol/m², RL), the leaf color gradually returned to partially green (Figure 1A). Like the changes of the leaf color, total anthocyanin and total flavone levels significantly increased under ML treatment and then decreased in RL compared with WL (Figures 1B,C). These results demonstrate that light intensity significantly affects the leaf color and the accumulation of total anthocyanins and total flavones in leaves of the light-sensitive perilla variety.

To assess anthocyanin components under different light treatments, we performed metabolomic analysis using UPLC-ESI-MS/MS. The samples were clearly clustered by PCA (Supplementary Figure S1). In total, 293 flavonoid metabolites were identified, including 103 flavonols, 72 flavones, 34 anthocyanins, and 84 other derivatives (Figure 2A; Supplementary Table S1). There were 90, 55, and 81 differentially expressed metabolites (DEM) in WL vs. ML, ML vs. RL, and WL vs. RL, respectively (Figure 2B).

Quantitative analysis showed that approximately 140 (48%) flavonoid metabolites were affected under different light intensity treatments in perilla leaves. Among these, 62 metabolites were significantly induced under increased light intensity, including 23 flavonols, 19 flavones, and 5 anthocyanins (Supplementary Table S2). Among them, the contents of malonylshisonin increased 1.7-fold under ML treatment compared with WL and decreased 0.8-fold under RL treatment compared with ML (Supplementary Table S3). Naringin chalcone and dihydroquercetin (early intermediates in the anthocyanin biosynthesis pathway), as well as cyanidin-3,5-O-diglucoside, cyanidin-3-O-(6”-O-p-coumaryl) glucoside-5-O-glucoside, and cyanidin-3-O-(6”-O-p-coumaryl) glucoside (important precursors of malonylshisonin) were all detected in this study. Interestingly, the levels of naringin chalcone, dihydroquercetin, and cyanidin-3-O-(6”-O-p-coumaryl) glucoside increased by 4.3-, 5.8-, and 3.2-fold (Supplementary Table S3), respectively, under ML treatment compared with WL. By contrast, the cyanidin-3,5-O-diglucoside contents decreased by 0.6-fold under ML treatment compared with WL. These results suggest that light intensity affects accumulation of the metabolites in anthocyanin biosynthesis pathway in perilla.

According to the change trends of flavonoid metabolites under different treatments, they were divided into eight categories. The levels of subclass 4 and subclass 7 members initially increased, then decreased, possessed consistent pattern with total anthocyanin contents. Both malonylshisonin and naringin chalcone belong to subclass 7. The levels of subclass 2 and subclass 6 members continuously increased under ML and RL treatment (Supplementary Figure S2). The cyanidin-3-O-(6”-O-P-coumaryl) glucoside-5-O-glucoside and dihydroquercetin belong to subclass 2 and subclass 4, respectively. The four subclasses attracted our main attention.

To explore the molecular regulatory mechanism of the perilla anthocyanin metabolites in response to light intensity, we performed transcriptomic sequence. Considering the rapid response of transcriptomic level to light treatment, we added a short treatment (12 h) under 180 μmol/m² of light (ML12h). After RNA sequencing and data filtering, 701.64 million clean reads with 105 Gb clean data were obtained. The Q20 reached 97%, with an average GC content of 48%. When we mapped the sequencing reads to perilla reference genomes, we obtained an average unique mapping ratio of 82.94% and a multiple mapping ratio of 10.19% (Supplementary Table S4).

Subsequently, we identified DEGs among all samples. And most of DEGs were responsed to light intensity showing to the results of GO analysis (Supplementary Figure S3). Compared to the WL group, 3290 and 3026 genes were upregulated in ML12h and ML6d, respectively. Compared to the RL group, 4043 and 3219 genes were upregulated in ML12h and ML6d, respectively (Figure 3A). Among the upregulated genes under ML treatment, 983 genes showed increased expression under short (12 h) and long (6 days) treatments compared with WL (Figure 3B). These DEGs are mainly involved in metabolic pathways, phenylalanine and tyrosine tryptophan biosynthesis and phenylalanine metabolism in the KEGG enrichment (Figure 3C). By contrast, 1478 genes were downregulated in RL compared with ML (Figure 3D). The down-regulated genes were enriched in biosynthesis of secondary metabolites, phenylalanine and tyrosine tryptophan biosynthesis, tyrosine metabolism, phenylalanine metabolism, and flavonoid biosynthesis in the KEGG enrichments (Figure 3E). Based on the expression pattern, we divided all DEGs into 26 categories (Supplementary Figure S4). The 983 and 1478 DEGs most belong to Profile 20 (upregulated in ML12h and ML6d compared with WL) and Profile 23 (downregulated in RL compared with ML12h and ML6d), respectively. These photo-response genes attracted our main attention.

We identified metabolites and genes involved in the anthocyanin and flavones biosynthetic pathways (Figure 4). Malonylshisonin biosynthesis involves 10 catalytic reactions, containing the initial cyanidin biosynthesis and the latter cyanidin derivatization pathway (Saito and Yamazaki, 2002; Jiang et al., 2020). Naringin chalcone, dihydroquercetin, cyanidin-3,5-O-diglucoside, cyanidin-3-O-(6”-O-p-Coumaroyl) glucoside, cyanidin-3-O-(6”-O-p-Coumaroyl) glucoside-5-O-diglucoside, and malonylshisonin were identified in the metabolomic analysis. Except for cyanidin 3,5-O-glucoside, these metabolites obviously accumulated under increased light intensity. Then, we identified 95 genes encoding the 10 anthocyanin synthases in the malonyshisonin biosynthetic pathway (Figure 4 and Supplementary Table S5). Furthermore, 32 genes were upregulated after light treatment, which encode the key rate-limiting enzymes in the basic anthocyanin biosynthetic pathway. Among these, 18 genes (including PAL, 4CL, CHS, F3’H, ANS, 5-GT, and MAT) were more strongly upregulated under 12 h ML treatment, and five genes (including F3H, 3-GT, ACT, and CHI) had higher expression levels under 6 days ML treatment. Only one C4H gene was expressed at high levels under both 12 h and 6 days ML treatments. These results suggest that light intensity remarkably influences the gene expression and metabolite accumulation involving in anthocyanin biosynthesis in perilla leaves.

Moreover, we identified 14 light-regulated flavones based on analysis of the flavonoid biosynthesis pathway. Most of these compounds accumulated under light treatment. Specifically, the contents of luteolin, apigenin, and baicalin increased by 5.0-, 5.4-, and 11.6-fold, respectively, under ML treatment compared with WL. Comparatively, the levels of apigenin-6,8-di-C-glucoside, luteolin-7-O-glucoside, and chrysin decreased under ML treatment. Because many genes involved in the accumulation of these substances are unknown, we only detected FSII, F6H, and UBGT in the flavone biosynthetic pathway (Supplementary Table S5). Like the malonylshisonin biosynthesis genes, both FSII and UBGT were upregulated under ML treatment compared with WL. In view of luteolin, apigenin and more flavones possess pharmacological effects (Lei et al., 2020; Liang et al., 2020; Rizzo et al., 2021), the results show that strongly accumulated flavones under elevated light treatment, with provided potential application for further research.

To gain further insight into the regulation of anthocyanin biosynthesis in response to light intensity, we performed WGCNA to identify co-expression networks between anthocyanin metabolites and gene expression levels. Ten co-expression modules were identified based on their similar expression patterns (Figure 5A). The turquoise module is positively correlated with anthocyanins, malonylshisonin, and other metabolites (Figure 5B). We therefore focused on those genes in the turquoise module. We identified anthocyanin biosynthesis and transcription factor genes in this module and used gene pairs with correlation coefficients more than 0.85 to construct a regulatory network (Figure 5C and Supplementary Table S7). Five important genes involved in anthocyanin biosynthesis were identified in the turquoise module, including ANS, 4CL, DFR, and CHI, which were positively correlated with total anthocyanin and malonylshisonin contents. The 147 transcription factors were also identified, mainly encoding MYB, bHLH, bZIP, ERF, and NAC in the module. Notably, MYB5, MYB12, bZIP2, bZIP6, bHLH1, bHLH7, and bHLH9 had high connectivity with malonylshisonin and five anthocyanin biosynthesis genes. The MYB_related2, MYB_related4 had high positive correlation with malonylshisonin. Moreover, the expression level of most identified MYB and bZIP genes were strongly upregulated under ML treatment compared with WL and downregulated under RL treatment compared with ML (Figure 5D). By contrast, MYB22, bHLH1, bHLH7, and bHLH9 were negatively correlated with malonylshisonin content. MYB22 and bHLH1 were also negatively correlated with total anthocyanin contents (Supplementary Table S6). Hence, these transcription factors were speculated to play important roles in regulating malonylshisonin accumulation under different light intensities.

MYB transcription factors are major members of MBW ternary transcriptional complex, which is an important regulator of anthocyanin metabolism in many plants (Lloyd et al., 2017; Li et al., 2020; Wang et al., 2020; Zhang et al., 2020; Ning et al., 2021). In WGCNA analysis, we identified 38 MYB transcription factors, containing 37 R2R3-MYB subfamily members and one 1R-MYB (MYB10). These MYBs showed high correlation coefficient with malonylshisonin and five anthocyanin biosynthesis genes (Supplementary Figure S5). Among them, 22 R2R3-MYB members were highly induced under ML treatment, while 15 were downregulated in response to light treatment (Figure 5D).

To examine the MYB transcription factors involved in regulating anthocyanin accumulation, we constructed a phylogenetic tree combined the MYBs that we identified in perilla and were previously reported to regulate anthocyanin accumulation in other plants. The phylogenetic tree was divided into three subclasses (Figure 6A). Subclass I contains perilla MYB12, MYB2 and key anthocyanin-regulated MYBs from Arabidopsis and apple. Subclass II includes 11 perilla MYBs and three Arabidopsis R2R3-MYB family members. Among them, MYB8 and MYB14 showed peak expression under ML12h treatment. Subclass III includes 25 MYBs from perilla and IbMYB340. Among them, MYB4 and MYB37 were expressed at high levels and responded to light treatment.

For further explore the relationships between the MYBs and anthocyanin biosynthesis, we performed correlation analysis of 15 selected MYBs with anthocyanin metabolites and 17 key anthocyanin biosynthesis genes. Among these, 9 MYBs showed positive correlations with anthocyanin metabolites and 6 showed positive correlations with key anthocyanin biosynthesis genes (Figure 6B). Interestingly, MYB12 and MYB19 showed significantly positive correlations with most anthocyanin biosynthesis genes (Figure 6C). The MYB12 had been reported in malonylshisonin formation, which was found in responsive to light treatment in the study.

Furthermore, the cis-acting elements of four key anthocyanin biosynthesis genes, including DFR, CHI1, CHI2, and ANS, were analyzed. Several light-response elements and hormone-response elements were identified in the four anthocyanin genes (Supplementary Figure S7). Interestingly, we identified CCAAT-box and MYB-recognition element (MRE) cis-acting elements in the promoter regions of perilla CHI1 and ANS, respectively (Supplementary Figure S7). Combining WGCNA analysis, phylogenetic analysis, and cis-acting element analysis, more MYB transcription factors were predicted to regulate malonylshisonin accumulation under the light treatment in perilla.

The bZIP transcription factors were recently reported to regulate anthocyanin biosynthesis (Hu et al., 2019; Liu et al., 2019; Wang et al., 2020). In the turquoise module, we also identified 6 bZIP transcription factors (Supplementary Figure S6). These perilla bZIPs were divided into two subclasses in the evolutionary tree (Figure 7A). Subclass I consist of bZIP4 and bZIP3 from perilla, along with HY5 transcription factors from other plants. Subclass II consists of 4 bZIP transcription factors from perilla, along with MdbZIP44, LcABF3, and MdABI5. Co-expression analysis showed that bZIP2, bZIP5, and bZIP6 share positive regulatory relationships with key anthocyanin biosynthesis genes (CHI1, CHI2, ANS, and DFR), while bZIP1, bZIP3, and bZIP4 share negative regulatory relationships with these genes (Figure 7B). Interestingly, bZIP2, bZIP5, and bZIP6 showed significant positive correlations with total anthocyanin and malonylshisonin contents and most anthocyanin biosynthesis genes (Figure 7C). The promoter sequences of CHI1, CHI2, ANS, and DFR all contain G-boxes, a light-responsive binding site that is recognized by bZIP transcription factors (Supplementary Figure S7). The bZIP2 and bZIP6 were identified to encode G-box binding factor 3 in perilla, suggesting that they may play key roles in the light response and the anthocyanin regulation.

We measured the expression levels of anthocyanin biosynthesis genes (CHI1, CHI2, ANS, and DFR) and important transcription factor genes (MYB, bZIP, and bHLH) involved in anthocyanin regulation by qRT-PCR (Figure 8). Combined the transcriptomic and qRT-PCR data, the candidate genes showed similar expression patterns. The anthocyanin biosynthesis genes (CHI1, CHI2, ANS, and DFR), MYB genes, and bZIP genes were upregulated under increased light intensity. These genes were all upregulated after a short (12 h) moderate light intensity treatment and downregulated after a long treatment (6 days), and their expression continued to decline under low light intensity. The response of the anthocyanin biosynthesis genes to light intensity treatment might be the direct cause of accumulation of anthocyanin metabolites. The MYB transcription factor genes were more sensitive to short vs. long light treatment, pointing to their important regulatory functions in anthocyanin metabolite accumulation. The bZIP2 and bZIP6 are homologous to the important light-response transcription factor gene HY5. Both were upregulated after a short light treatment. We also identified bHLH1 and bHLH7 transcription factor genes in perilla. Interestingly, both genes were downregulated under light treatment. The bHLH1 and bHLH7 encode PIF3 and PIF7, respectively. The PIFs were negatively regulation factors in response to light in other plants. The downregulated trend of bHLH1 and bHLH7 under light treatment is consistent with their functions. The expression patterns of key anthocyanin biosynthesis genes and transcription factors in response to light points to their important roles in regulating anthocyanin biosynthesis.