Eight cultivars of A. andraeanum (Hort.) were grown in 20-cm plastic containers filled with a soilless substrate composed of 60% peat mixed with coconut chaff and perlite (STANLEY, Linyi, China) in a shaded greenhouse under a maximum photosynthetic radiation of 200 μmol m^-2 s^-1 in Guangzhou Flower Research Center, Guangzhou, Guangdong Province, China. Daytime and nighttime temperatures in the shaded greenhouse were maintained from 25 to 28°C and 19 to 21°C, respectively with a relative humidity ranging from 70 to 90%. Plants were fertilized and watered as described by Chen et al. (2012). Pinching of flower buds was used to promote vegetative growth in the initial flowering stage. Subsequently, plants entered their bloom stage.

Spathes of these cultivars varied in color ranging from red to white. For transcriptomic analysis of spathe colorations, three plants were randomly selected per cultivar, one spathe was collected from each plant at stage 7 of spathe development, i.e., peduncle fully extended but spathe showing full coloration without being open ( Supplementary Figure S1 ). Thus, there were three biological replicates per cultivar. To analyze anthocyanin, chlorophyll, and carotenoid levels in fully open spathes, spathes at their developmental stage of 10 ( Supplementary Figure S1 ), i.e., spathes were fully expanded but flowers were not dehisced (Collette et al., 2004), were collected (about 200 mg) from three randomly selected plants per cultivar. Collected spathes were immediately extracted with appropriate solutions mentioned below.

Total RNA was extracted from collected spathes. RNA concentrations were measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies,CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and also checked by electrophoresis using 1% agarose gel ( Supplemental Table S1 ). The mRNA was enriched using magnetic bead with oligo (dT) and used for cDNA synthesis. The first strand cDNA was synthesized using random hexamer primer and M-MuLV reverse transcriptase (Dingguo Changsheng, Inc. Beijing, China). The second strand cDNA synthesis was subsequently performed using DNA polymerase I with RNase H. cDNA fragments were preferentially selected from 150 to 200 bp in length. PCR amplification was performed for 13 cycles as recommended. The library preparation was performed according to the method described by Liu et al. (2021). The cDNA libraries were sequenced by Illumina Novaseq™ system at NextOmics Bio-Tech. Co. (Wuhan, China). Paired-end sequencing of the polyA enriched library was performed on an Illumina Novaseq 6000 platform generating 150 nt reads.

RNA-Seq data were analyzed to identify key structural genes and transcription factors involved in flavonoid biosynthesis and chlorophyll metabolism. Trinity software (version 2.4.0, Supplementary Table S2 ) was used to assemble the short reads (Grabherr et al., 2011). Trinity commands were as follows: Trinity –seqType fq –max_memory 120G –CPU 30 –min_kmer_cov 2 –full_cleanup –left clean.R1.fq.gz –right clean.R2.fq.gz –output trinity.out. Other parameters were default. Based on the sequence similarity principle, TIGR Gene Indices clustering tools (TGICL, version 2.1, Supplementary File S1 ) was further clustered to eliminate redundancy (Pertea et al., 2003). CDS and protein sequences were predicted by Transdecoder software (version 5.7.0, https://github.com/TransDecoder/TransDecoder/releases, default parameters); then, a unigene sequence set of 8 cultivars (24 samples) was obtained ( Supplementary Table S3 ). The analysis on the transcriptome assembly to assess completeness was conducted using BUSCO (version 5.4.2) (Manni et al., 2021). The lineage dataset is embryophyta_odb10.

Unigene functions were annotated based on the following database: the NCBI COG/KOG database (http://www.ncbi.nlm.nih.gov/COG/), Kyoto Encyclopedia of Gene and Genomes database (KEGG, https://www.genome.jp/kegg/), non-redundant protein database (NR, ftp://ftp.ncbi.nih.gov/blast/db/), SwissProt protein database (https://www.uniprot.org/), and Gene Ontology (GO, http://geneontology.org/). Identification of flavonoid biosynthesis genes was annotated using KIPEs (Knowledge-based Identification of Pathway Enzymes, https://pbb-tools.de/KIPEs/) based on transcript sequences (Pucker et al., 2020). MYB were also predicted with a dedicated tool MYB annotator (Automatic annotation of MYBs, https://pbb-tools.de/MYB_annotator/) (Pucker, 2022). bHLH were predicted with a dedicated tool bHLH annotator (Automatic annotation of the bHLH gene family in plants, https://pbb-tools.de/bHLH_annotator/) (Thoben and Pucker, 2023). The identifications of flavonoid biosynthesis genes, MYB and bHLH were listed in Supplementary File S2 .

Gene expression levels were estimated by RSEM software (version 1.1.12) and presented as FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) (Li and Dewey, 2011). Differentially expressed genes (DEGs) were analyzed between different groups by DESeq2 software (version 1.38.3) (Love et al., 2014). DEGs were detected and screened between the eight comparisons, using fold change ≥ 2 and false discovery rate (FDR) ≤ 0.05. The “pheatmap” R package (version 4.2.2, https://www.r-project.org/) was used to determine DEGs abundance, and the scale of heatmap for DEGs were converted from the FPKM was Log2 of the fold change. The enrichment analysis was performed in GO term and KEGG pathway.

For validation of candidate genes associated with the spathe coloration, red spathe cultivars Te Lun Sa (TLS), 2016, and A302; green spathe cultivars Kai Xin Guo (KXG), A166, and A086; and white spathe cultivars Bai Ma (BM), A231BAI, and A168 at stage 7 ( Supplementary Figure S1 ) were harvested, immediately immersed in liquid nitrogen, and stored at -80°C for qRT-PCR analysis. These nine cultivars were grown in the same conditions as the aforementioned eight cultivars.

Anthocyanin is red in acid solution. Its color is in proportion to the anthocyanin content. The absorption peak of anthocyanin acid solution is 530, with which molar extinction coefficient is 4.62 × 104. Thus, ultraviolet-visible (UV-Vis) spectrophotometer can be used to measure the anthocyanin content.

To analyze anthocyanin, collected 0.2 g spathes from eight cultivars were extracted in 10 mL of ethanol: hydrochloride (99:1, v/v) at 32°C for 4 h. The sample was centrifuged at 5,000 rpm for 5 min (Centrifuge 5804 R, Eppendorf, Shanghai, China), and the supernatant was stored at 4°C. The residue was extracted with the extraction solution 1 to 2 times until the supernatant turned colorless. Absorbencies were measured at 530 nm and 650 nm wavelengths with an UV-Vis spectrophotometer (UV4800, Unico, Shanghai, China). The relative content of anthocyanin was calculated according to the method (Yang et al., 2022):

For quantification of total chlorophyll and carotenoid, fresh spathe samples were ground to powder in a mortar with liquid nitrogen, extracted with 80% acetone, filtered into a test tube, and finally extracted with ethyl acetate. The suspension was centrifuged for 5 min at 2,000 rpm. The absorbance of the supernatant was recorded at 663 nm, 645 nm, and 440 nm, respectively. The total contents of chlorophyll (Chl) and carotenoid were calculated according to the following formula (Zhou et al., 2013; Vaculík et al., 2015):

The UV/Vis spectrophotometry was used to evaluate the proanthocyanidins (PA) content. About 200 mg spathe was used for extraction. Soluble PAs from fresh anthurium spathes were extracted and measured by Innovabio Testing Technology Co. (Nanjing, Jiangsu). After ethanol precipitation, the PA content was determined by normal butyl alcohol-hydrochloric acid method at 546 nm. Proanthocyanidins (Cas: 4852-22-6, Yuanye, Shanghai, China) was used as a standard for PA quantification. The content of PA was calculated based on the standard curve according to the method (Li X. et al., 2021).

About 100 mg spathe was ground to powder and extracted by ethyl alcohol: distilled water: hydrochloric acid (2:1:1, v/v/v) for supernatant preparation. SCIEX Qtrap6500 mass spectrometer System (AB, Massachusetts, USA) coupled to a LC-30AD high performance liquid chromatograph (Shimadzu, Kyoto, Japan) was used for liquid chromatography-tandem mass spectrometry (LC-MS/MS). Cyanidin chloride (Cas: 528-58-5), and pelargonidin chloride (Cas: 134-04-3, Yuanye, Shanghai, China) was used as a standard for cyanidin and pelargonidin quantification, respectively. The standard substance was dissolved in hydrochloric acid-methanol solution as single standard stock solution. LC-MS/MS analysis was performed by Innovabio Co. (Nanjing, Jiangsu) according to the method of Tomas (2022).

Spathes at stage 7 ( Supplementary Figure S1 ) harvested from the above-mentioned cultivars were ground to a powder with a mortar and a pestle in liquid nitrogen. RNA was extracted using the RaPure Plant RNA Kit (Magen, Guangzhou, China) following the instructions of the manufacturer. Total RNA was quantified using NanodropND1000 UV/Vis spectrophotometer (Thermo Scientific, Walthman, MA), and the purity was assessed by the absorbance ratios of 260/280 nm and 260/230 nm. The integrity of the purified RNA was confirmed with 2% agarose gel electrophoresis. The first-strand cDNA was synthesized from 1 μg of total RNA using a TIANGEN FastKing cDNA Kit (TIANGEN, Beijing, China) according to the manufacturer’s instructions. The 20 μL cDNA samples were diluted to a 1:5 ratio using 80 μL RNase-free water for all gene amplification reactions.

qRT-PCR reactions were performed using SYBR Green detection chemistry by running on 96 well-plates with a CFX connect Real-Time PCR System (BIO-RAD, Shanghai, China) with three replicates for each sample. Primers for each target were designed using Primer3 web service (https://primer3.ut.ee/), and the sequences of all primers are listed in Supplemental Table S4 . qRT-PCR analysis data were averaged by three biological replicates and calculated using the 2^-△△Ct method (Livak and Schmittgen, 2001). All expression levels were normalized with internal control, AaGAPDH (Li S. et al., 2021).

The upstream sequences (5,000 bp) of AaLAR (GenBank accession no. OR241534), AaANR (GenBank accession no. OR241535), AaHemB (GenBank accession no. OR241536), and AaPor (GenBank accession no. OR241537) genes were obtained from the Anthurium Genome Database of A. andraenum cv. Xiaojiao-Texana. Potential binding motifs for TFs of MYB family on AaLAR, AaANR, AaHemB, and AaPor promoter regions were predicted by the online website JASPAR²⁰²² (http://jaspar.genereg.net/) (Fornes et al., 2019). The TF binding sites and motifs were further screened by a relative profile score threshold of 90%. The binding sites on the positive strand were presented and the common MYB motifs in the anthocyanin or chlorophyll biosynthesis pathway were marked. As the motifs in the JASPAR database were only included in the model species, Arabidopsis thaliana, Oryza sativa, and Zea mays, the MYB proteins in A. andraeanum were aligned together by Clustalw for conducting phylogenetic analysis using the neighbor-joining method in MEGA version 7.0 to determine the genetic relationship ( Supplementary Figures S2 , S3 ).

Data for pigment contents and gene relative expression levels were statistically analyzed using SPSS (SPSS statistics 19, IBM Corp., Armonk, NY, USA). If significance occurred among cultivars, means were separated using Fisher’s LSD at P< 0.05 level. Data are presented as mean ± S.E. with either 3 or 5 replicates.

Spathes of eight anthurium cultivars had eight different colors, including ‘Rosa’ (Ros) with red spathe; ‘Wu Yang’ (WY), orange spathe; ‘Sante’ (San), coral spathe; ‘Guang Hua Zi Yun’ (GHZ), purple spathe; ‘Pink Champion’ (Pin), pink spathe; ‘Sonate’ (Son), light pink; ‘Midori’ (Mid), green spathe; and ‘Acropolis’ (Acr) ( Figure 1A ). Anthocyanin contents in spathes of Ros and GHZ were the highest, 14.0 and 15.0 ΔA/g·FW, respectively, which were substantially greater than the others ( Figure 1B ). Anthocyanin content in Acr was the lowest, only 0.5 ΔA/g·FW. The other cultivars had anthocyanin ranging from 1.0 to 3.8 ΔA/g·FW, suggesting that anthocyanin contents were associated with redness of spathes. The main component of anthocyanins in GHZ was cyanidin, up to 0.25 mg·g^-1·FW, rather than pelargonidin, which was only 0.025 mg·g^-1·FW ( Supplementary Figure S4A ). Proanthocyanidin (PA) contents of Ros and San were highest, 11.4 and 14.5 mg·g^-1·FW, respectively, following by Son and Mid. PA content in Acr was about 5.27 mg·g^-1·FW ( Supplementary Figure S4B ).

The Chl and carotenoid levels were the highest in Mid at 0.16 ΔA/g·FW and 0.08 ΔA/g·FW, respectively ( Figures 1C, D ). Notably, the contents of both Chl and carotenoid in the other seven cultivars were not significantly different, ranging from 0.02 to 0.04 ΔA/g·FW for Chl, and 0.002 to 0.02 ΔA/g·FW for carotenoid, respectively. These results showed that the pigment composition and content in green and white spathes were different from those of the other colored spathes.

Raw data or reads generated from the transcriptome sequencing of 24 samples (8 cultivars, each with three biological replicates) ranged from 40,964,478 to 63,961,814. After stringent quality control checks and filtering, a total of 1,947,453,054 paired-end clean reads were obtained with Q30 values higher than 92% in each sample library. The data volume per library was greater than 6 Gb. The clean reads were spliced and assembled using Trinity. N50 length of the Trinity assembly ranged from 1610 to 1799 bp ( Supplementary Table S2 ). Trinity sequence of all samples were merged. According to the sequence similarity principle, the merged Trinity assembly sequences were clustered for a set of unigene sequences, which resulted in a total of 62,013 unigenes. The length of assembled unigenes varied from 220 to 21,222 bp. The average length of the unigenes was 2,038.45 bp, and the N50 was 2,519 bp ( Supplementary Figure S5A , Supplementary Table S3 ). These data showed that the throughput and sequencing quality were high enough for further analysis. The complete BUSCOs of the transcriptome assembly was 1387, accounting for 85.9% ( Supplementary Figure S6 ). In this study, the tissue used for RNA-Seq was spathes only. Some genes were hardly detectable in the spathe, suggesting that some genes were not expressed in this special organ. Thus, it was reasonable for the incompleteness of this transcriptome assembly.

A total of 48,582 coding sequences (CDS) were predicted, and their lengths varied from 261 to 15,300 bp with the majority ranging from 972 to 1,206 bp ( Supplementary Figure S5B ). These results implied that the reads derived from transcriptome sequencing had high quality and reliability. Further, 33,392 unigenes were assigned to 23 categories in the KOG database ( Supplementary Figure S5C ). Among them, 6,806 unigenes, “Posttranslational modification, protein turnover, chaperones” represented the largest group with 3,453 unigenes, followed by “Signal transduction mechanisms” (2,858 unigenes)”, and Transcription (1,965 unigenes). In the NR database, the identified unigenes were used to detect their homologous genes in other species by blast, and the results showed that the ratio of homologous genes in Elaeis guineensis, Phoenix dactylifera, Nelumbo nucifera, Musa acuminata subsp. Malaccensis, and Ananas comosus were 25.60%, 21.40%, 8.53%, 5.92%, and 4.59%, respectively ( Supplementary Figure S5D ).

Transcriptome analysis identified differentially expressed genes (DEGs) in different colored spathes based on their FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) levels. To identify the key DEGs involved in different colorations, seven pair-wise comparisons were conducted which included Ros vs. Acr; WY vs. Acr; San vs. Acr; Mid vs. Acr; GHZ vs. Acr; Pin vs. Acr; and Son vs. Acr ( Supplementary Figure S7A ). Among these comparisons, the greatest abundance of DEGs (29,588) was found in Pin and Acr libraries, of which 13,891 and 15,769 genes were up-regulated and down-regulated, respectively. Conversely, the lowest abundance of DEGs (14,015) were recorded in GHZ and Acr libraries ( Supplementary Figure S7A ). Gene functional annotation was carried out to understand the function of DEGs. According to the GO annotation, a total of 51,469 DEGs were spread across 1,123 terms consisting of three domains: molecular functions accounting for 57.52%, cellular components 12.54%, and biological processes 29.94%. The most frequent terms under the molecular function were protein binding, ATP binding, protein kinase activity, and nucleic acid binding, DNA binding, and zinc ion binding ( Supplementary Figure S7B ). To further examine the metabolic process of the DEGs, the KEGG pathway for DEGs in seven comparison pairs were conducted. The top 20 enriched pathways were listed in all comparable groups ( Supplementary Figure S7C ). From the results of KEGG enrichment analysis, the pathways for flavonoid biosynthesis, and porphyrin and chlorophyll metabolism were enriched in the comparisons of Ros vs. Acr; GHZ vs. Acr; WY vs. Acr; and Mid vs. Acr ( Supplementary Figures S7C , S8 ).

The DEGs associated with the flavonoid biosynthesis were analyzed through the examination of KEGG enrichment and gene annotations, from which 56 DEGs in the flavonoid pathway were identified, including nine MYB genes and seven bHLH genes regulating flavonoid biosynthesis, and one C4H, one 4-coumarate-CoA ligase-like (4CL), five CHS, four CHI, one F3H, two F3’H, one flavonol synthase (FLS), four DFR, six LAR, one ANS, three anthocyanidin 3-O-glucosyltransferase (3GT), two anthocyanidin 5,3-O-glucosyltransferase (5,3GT), four ANR, two caffeoyl-CoA O-methyltransferase (CCoAOMT), and three glutathione S-transferase F12 (also named TT19) in seven unrepetitive comparison pairs of Ros, WY, San, GHZ, Pin, Son, Mid, and Acr. According to the expression level of DEGs shown in heatmap ( Figure 2 ), MYBs were highly expressed in Ros, WY, San, and GHZ compared to their expression in the other cultivars, of which MYB2 Unigene32682 exhibited notable changes with the cultivars or the content of anthocyanin in spathes ( Figure 2 , Supplementary Figure S9 ). The results implied that MYB genes could play an important role in regulation of anthocyanin biosynthesis in anthuriums. For CHI, its expression levels in Pin, Son, Mid, and Acr were relatively lower than in other cultivars ( Figure 2 ), indicating that more naringenin chalcones was transformed into naringenin in Ros, WY, San, and GHZ. The expression level of F3’H in Son was significantly higher than other spathes, suggesting that more dihydroquercetin (DHQ) was synthesized in Son than the other spathes. F3H gene expression in San was the highest among eight anthurium cultivars, indicating that naringenin in San was catalyzed into more dihydrokaempferol (DHK). Higher levels of DFR expression occurred in Ros, WY, and GHZ, converting more dihydroquercetin/dihydrokaempferol to leuco-cyanidin/pelargonidin. The expression of ANS was increased in WY, San, and Son, with cyanidin and pelargonidin accumulation. The gene expression of 5,3GT was highly detected in WY and GHZ, accompanied with 3,5-O-glucoside product. Downregulation of a decorating enzyme CCoAOMT gene seemed to activate the anthocyanin biosynthesis in red- and purple-colored spathes. Anthocyanin transport carrier TT19 was induced in Ros, WY, San, GHZ, Pin, and Son, supporting anthocyanin accumulation and especially cyanidin accumulation. The increased expression of the LAR gene, such as Unigene59237 ( Figure 2 ), could result in converting leucocyanidin and leucopelargonidin to catechin and afzelechin, respectively. The highly expressed ANR gene, especially Unigene7781, could make cyanine and pelargonidin flow into epicatechin and epiafzelechin, respectively. The results showed that the expression of LAR and ANR were significantly different in green and white spathe cultivars compared with the six other cultivars. The increased expression of LAR and ANR could reduce the catalytic substrate to synthesize anthocyanins, resulting in accumulation of proanthocyanidins in spathes.

Except for differential expression of flavonoid biosynthesis genes, we also analyzed if MYB-TFs were implicated in the regulation of some of the genes. A MYB-TF (Unigene32682, MYB2), a LAR gene (Unigene59237), and a ANR gene (Unigene7781) were chosen from DEGs to analyze their expressions in the eight cultivars ( Figure 2 ). The transcript level of MYB2 was significantly lower in the green or white spathed cultivar compared with the color intensified ones. Conversely, in comparison to the color intensified spathes, green and white spathes had significantly high levels of LAR and ANR transcripts ( Figure 3A , Supplementary Figure S10 ).

Promoter analysis of LAR and ANR genes showed that the promoter of LAR contained two motifs bound by MYB2, and the promoter of ANR also had one motif bound by MYB2 ( Figure 3B , Supplementary Figure S2 ), implying that the expression levels of LAR and ANR could be likely regulated by the TF-MYB2 ( Figure 3B ). In addition, expression analysis revealed that the expression pattern of MYB2 was opposite to that of LAR and ANR, suggesting MYB2 might negatively regulate the expression of LAR and ANR ( Figure 3A ).

bHLH family members were also analyzed. Expression levels of bHLH113 and bHLH49 increased in purple spathe GHZ, while bHLH62 and bHLH57 were highly detected in red spathe Ros ( Figure 2 ), indicating bHLH might be involved in the positive regulation of anthocyanin biosynthesis.

Comparative analysis of DEGs showed that the expression levels of a set of MYB-TFs, including MYB124 (Unigene12793) were significantly higher in Mid (green spathe) than in the others, which could suggest that those MYB-TFs may also participate in the regulation of Chl synthesis ( Figure 4A ). Thus, the DEGs in the chlorophyll biosynthesis pathway were further analyzed, and data indicated that the expressions levels of HemB (Unigene51685) and Por (Unigene38255) were significantly greater in Mid than in the other cultivars ( Figure 4B , Supplementary Figure S10 ). To explore the expression pattern of the DEGs, qRT-PCR was conducted with three genes, MYB124, HemB, and Por, involved in the biosynthesis of Chl. The results showed that the relative expression levels of the three genes were similar to the FPKM values ( Figure 4C ), suggesting that MYB124 was positively correlated with the expression of HemB and Por in green spathe. These data implied that MYB124 could promote Chl synthesis by regulating HemB and Por genes.

To further investigate the regulatory relationships among MYB124, HemB, and Por, we analyzed the promoter sequences of HemB and Por, respectively. The analysis showed that upstream of HemB and Por contained MYB124-bounding motifs, respectively ( Figure 4D , Supplementary Figure S3 ). The results supported the above analysis that TF-MYB124 enhanced the expression of the downstream genes, HemB and Por and promoted Chl accumulations, thus the green spathe in cultivar Mid.

Since Mid with green spathe was loaded with carotenoids ( Figure 1D ), we took a closer look at gene expression in carotenoid biosynthesis pathway. Five unigenes were identified, including one PDS (Unigene54090), two LUT5 (Unigene14810 and Unigene58790), one CrtR-b (Unigene54067), and one ZEP (Unigene8259). Upstream gene phytoene desaturase (PDS), which is associated with carotenoid content, was detected to be highly expressed in Mid. Increased expression of lutein deficient (LUT5) and β -carotene hydroxylase (CrtR-b) in Mid suggested more lutein accumulation ( Supplementary Figure S11 ). Meanwhile, the level of zeaxanthin epoxidase (ZEP) was significantly greater in Mid than other cultivars, indicating that ZEP may be involved in carotenoid synthesis through the epoxidation of β-carotene and β-cryptoxanthin ( Supplementary Figure S11 ). Thus, increased carotenoid content in Mid might be caused by these elevated expression genes in carotenoid biosynthesis pathway.

To further confirm key genes involved in the change of spathe colors, we used qRT-PCR to analyze the expression pattern of seven candidate genes in nine additional cultivars with red, green, and white spathes at S7 stage ( Figure 5 , Supplementary Figure S1 ). As expected, the expression levels of LAR and ANR in green or white spathes were significantly higher than in red spathes. In contrast, MYB2 gene exhibited the opposite expression pattern of LAR and ANR ( Figure 5 , Supplementary Figure S9 ). Therefore, these results confirmed that the expression of MYB2 was negatively correlated with the expression of LAR and ANR.

HemB and Por genes were highly expressed in green spathes but at very low levels in red and white spathes ( Figure 5 ). Such expression levels were closely associated with Chl contents in the respective spathes ( Figure 1C ). Additionally, MYB124 was highly correlated with HemB or Por expression and the total Chl content, suggesting that HemB and Por may be regulated by MYB124 in regulation of Chl biosynthesis.

Taken together, these results showed that the white colored spathe was probably due to the increased expression of LAR and ANR, resulting in the production of colorless PA. The green-colored spathe was likely attributed to the elevated expression of HemB and Por, leading to chlorophyll biosynthesis. While MYB2 could negatively regulate the expression of LAR and ANR, and MYB124 enhanced the expression of HemB and Por.