BLASTN and BLASTP alignments (coverage = 60%, identity = 60%, e-value = 1.00e-20) were performed on 31 species of Brassicaceae, using the MYBL2 (AT1G71030) nucleic acid sequence and the protein sequence of Arabidopsis thaliana as references to obtain homologous sequences. The protein sequence of the homologous gene was extracted and submitted to BatchCDD (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) for conservative domain analysis, and the results were compared with those from BLAST to determine the gene family members. The IQ-TREE software was employed to construct the phylogenetic tree with maximum likelihood, while iTOL (https://itol.embl.de) was utilized to enhance the visual representation of the phylogenetic tree. In six Brassica species, SynOrthos software was used to analyze the collinear relationships between A. thaliana and Brassica species (Cheng et al., 2012), and TBtools-II software (v2.118) facilitated the visual analysis (Chen et al., 2023). The coding sequences (CDS) and protein sequences of MYBL2 family members were retrieved from 30 Brassicaceae genomes, and the Ka/Ks values for each MYBL2 family member were computed using TBtools-II software (v2.118) (Chen et al., 2023). The R software package pheatmap was used to create the heat map for MYBL2 members with Ka/Ks values.

MEME (https://meme-suite.org/meme/) was employed to analyze the conserved characteristics of the MYBL2 protein sequence, while BatchCDD (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) was utilized to assess the conserved domains of the MYBL2 protein. The 2 kb sequence upstream of the MYBL2 start codon was extracted and submitted to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for promoter cis-acting element prediction analysis, with visualization performed using TBtools-II (v2.118) software (Chen et al., 2023). The protein sequence of the Brassica MYBL2 homologous gene was submitted to CLUSTALW (https://www.genome.jp/tools-bin/clustalw) for sequence alignment, and the resulting protein sequence alignments were visualized using ESPript3.0 (https://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi) and WebLogo3 (https://weblogo.threeplusone.com), respectively.

To investigate the expression patterns of the MYBL2 gene across various species of Brassica and different tissues of B. napus, we collected transcriptome data from the leaves (green and purple) of B. rapa, B. juncea, and B. oleracea. Additionally, the transcriptome data encompassed five distinct tissues of B. napus, which included siliques (green and purple), stems (green and purple), leaves (green and purple), flower colors (white and purple), seed coats (yellow and black), as well as six different flower colors (white, beige, yellow, purple, apricot pink, and orange). Furthermore, we obtained transcriptome data from four different developmental stages (S1, S2, S3, S4) of purple flower B. napus. Each group of RNA-seq data above contained three biological replicates. The raw RNA-seq data obtained were mapped to the pak choi reference genome (Thakur, 2024) using HISAT2 software (v2.1.0), ensuring that each read matched only one region. The htseq-count function of the HTSeq software (v0.11.2) package was employed to count the number of reads aligned to each gene (Anders et al., 2015). Subsequently, the DEseq2 package (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) was utilized to analyze differentially expressed genes (DEGs). The expression level of each gene was calculated using StringTie software (v2.1.1) (Thakur, 2024). Additionally, TBtools-II software (v2.118) was used to create a heat map of gene expression related to the anthocyanin synthesis pathway (Chen et al., 2023).

To further investigate the hub genes involved in regulating anthocyanin synthesis in the B. napus different tissues, we employed WGCNA to construct the interaction relationships among ABGs across leaves (green and purple), stems (green and purple) and flowers (white and purple) of B. napus. Correlation analysis was conducted using the WGCNA software package (v4.3.1) (Langfelder and Horvath, 2008) to assess the relationships between each co-expression module, phenotypic traits, and anthocyanin content, thereby identifying key hub genes that influence anthocyanin synthesis. Subsequently, the gene regulatory network data derived from the WGCNA analysis were imported into Cytoscape (v3.80) for visual representation (Shannon et al., 2003). For detailed methodology, please refer to the article by Liu et al (Liu et al., 2024).

The leaves (green and purple), stems (green and purple) and flowers (white and purple) were collected and then used to extract total RNA, from which cDNA was synthesized through reverse transcription for quantitative real-time PCR (qRT-PCR) analysis. The specific methodology is as follows: First, the Eastep^® Super Total RNA Extraction Kit from Promega Biotech Co., Ltd (Beijing, China) was utilized for total RNA extraction. Subsequently, the HiScript III 1st Strand cDNA Synthesis Kit (+gDNA wiper) from Vazyme Biotech Co., Ltd (Nanjing, China) was employed for cDNA synthesis via reverse transcription. The qRT-PCR was conducted using the SYBR qPCR Master Mix from Vazyme Biotech Co., Ltd (Nanjing, China). The reaction mixture comprised a total volume of 20.0 µL, consisting of 10.0 µL SYBR qPCR Master Mix, 0.5 µL each of forward and reverse primers, 4.0 µL of cDNA template, and 5.0 µL of ddH₂O. The amplification procedure included an initial denaturation step at 95 °C for 30 seconds, followed by 40 cycles of denaturation at 95 °C for 10 seconds and annealing/extension at 60 °C for 30 seconds. This was followed by a final denaturation at 95 °C for 15 seconds, annealing at 60 °C for 60 seconds, and a melting curve analysis at 95 °C for 15 seconds. The relative expression levels were calculated using the 2 ^-ΔΔCT method, with the fold change determined based on these relative expression levels. The average values were visualized using Prism 9 (v9.5.1) software. Primer information for the qRT-PCR analysis is provided in Supplementary Table S1 . For detailed methods, please refer to Chen et al (Chen, 2020).

To systematically identify MYBL2 genes in Brassicaceae plants, we utilized the MYBL2 protein sequence of A. thaliana (AT1G71030.1) as a reference. This approach enabled the identification of MYBL2 genes across 31 Brassicaceae plants with published reference genomes. The results revealed the presence of 48 MYBL2 genes, excluding Aethionema arabicum ( Table 1 ; Supplementary Table S2 ). Among these, B. napus exhibited the highest number of homologous MYBL2 copies, totaling four. Additionally, three copies were found in B. carinata, B. juncea, Raphanus sativus, and Lepidium meyenii. Two homologous MYBL2 copies were identified in B. nigra, B. oleracea, B. rapa, Camelina sativa, Leavenworthia alabamica, Orychophragmus violaceus, and Sinapis arvensis, while only one MYBL2 homologous gene was detected in the remaining species ( Table 1 ).

To investigate the evolutionary relationships among MYBL2 homologous genes across 31 species of Brassicaceae, we constructed a phylogenetic tree using protein sequences. The phylogenetic analysis revealed that 48 MYBL2 homologous genes can be classified into six subfamilies (Group I to Group VI) ( Figure 1 ). Notably, the rna-MERR-LOCUS39516 from Microthlaspi erraticum forms a separate group (Group I), while Group II comprises two genes: MIN02G3220.t1 from Matthiola incana and KFK41528.1 from Arabis alpina. Group III includes 16 members of the MYBL2 gene family, among which AT1G71030.1 is represented. Additionally, Thlar.0031s0691.1 from Thlaspi arvense and Thhalv10019036m from Eutrema salsugineum are classified within Group IV. Groups V and VI contain 12 and 15 genes, respectively, with 16 MYBL2 homologous genes from six species of Brassica distributed across these two groups. Ka/Ks analysis found that only 20 gene pairs had Ka/Ks values > 1, indicating that they were subject to positive selection, while the Ka/Ks values of other gene pairs were < 1, indicating that these genomes were mainly subject to purifying selection during evolution ( Supplementary Figure S1 ; Supplementary Table S3 ).

In six Brassica species of U’s triangle, a total of 16 MYBL2 homologous genes were identified ( Table 1 ). To clarify the evolutionary relationships among 16 MYBL2 homologous genes in six Brassica species, we conducted phylogenetic and collinearity analyses of these genes ( Figure 2 ; Supplementary Table S4 ). The results indicated that the 16 MYBL2 homologous genes were categorized into two primary branches, which could further be subdivided into three smaller branches corresponding to the three subgenomes: A, B, and C ( Figure 2A ; Supplementary Table S5 ). Notably, within these three subgenomes, the MYBL2 homologous genes located on chromosomes A07, C06, and B03 clustered within the same evolutionary branch. In contrast, the MYBL2 homologous genes found on chromosomes A02, C02, C03, and B05 were allocated to a different evolutionary branch, suggesting that the homologous genes in these two branches may have diverged early in their evolutionary history. Furthermore, we examined the co-linearity of these 16 MYBL2 homologous genes across the Brassica A, B, and C subgenomes. The findings revealed that in the A subgenome, the MYBL2 homologous genes were predominantly located on the A02 and A07 chromosomes of B. napus and B. rapa, while B. juncea contained only one homologous gene on the A02 chromosome. In the C subgenome, MYBL2 homologous genes were mainly found on chromosomes C02 and C06 of B. napus and B. oleracea, with only one homologous gene from B. carinata located on chromosome C03. In the B subgenome, MYBL2 homologous genes were primarily distributed across the B03 and B05 chromosomes of B. nigra, B. juncea, and B. carinata ( Figure 2B ; Supplementary Table S6 ).

To comprehensively study the evolutionary patterns of 16 MYBL2 homologous genes in Brassica, we analyzed the promoter cis-acting elements, gene structure, conserved motifs, and conserved domains of these genes ( Figure 3 ). Initially, we predicted the cis -acting elements within the 2 Kbp promoter sequence upstream of the 16 MYBL2 homologous genes, identifying a total of 39 distinct cis-acting elements ( Supplementary Table S7 ). These elements were classified into categories including anaerobic induction, auxin responsiveness, cell cycle regulation, circadian control, defense and stress responsiveness, dehydration, low-temperature and salt stresses, low-temperature responsiveness, endosperm expression, gibberellin responsiveness, abscisic acid responsiveness, light responsiveness, MeJA responsiveness, meristem expression, drought inducibility, MYBHv1 binding site, salicylic acid responsiveness, wound responsiveness, and zein metabolism regulation, totaling 18 types of cis-regulatory elements. Notably, cis -acting elements associated with light response constituted half of the total. The distribution of components within the same branch exhibited similarities ( Figure 3A (a)). The gene structure of the 16 MYBL2 homologous genes was highly conserved; with the exception of the deletion of the first exon in the gene Bol.MYBL2.C06, the remaining 15 homologous genes comprised three exons, and the total gene length, as well as the length and location of each exon, were highly similar ( Figure 3A (b)). MEME online software was utilized to analyze 16 homologous MYBL2 proteins motifs, resulting in the identification of a total of 10 distinct motifs. The analysis revealed that all homologous genes contained motifs 1, 3, and 4; however, the gene BjuMYBL2.A02 uniquely lost motif 2 and exhibited a partial deletion at the N-terminus ( Figure 3A (c); Supplementary Table S8 ).

Interestingly, both the cis-acting elements of the promoter and the gene structure, including conserved motifs and domains, reveal that the gene BjuMYBL2.A02 is significantly different from other members of the Brassica family. To further investigate this, we compared the protein sequences of these 16 homologous genes. The results indicated that the gene BjuMYBL2.A02 lacks a complete R3-MYB structure ( Figure 3B ), while the remaining 15 homologous genes possess a complete R3-MYB structure. Additionally, motif analysis of the 16 MYBL2 protein sequences in Brassica demonstrated that they exhibit a typical R3 conserved motif, EAR motif (LPDLNI(S/G)L) and a terminal TLLLF motif ( Figure 3C ). These findings suggest that, with the exception of BjuMYBL2.A02, MYBL2 is highly conserved across Brassica species.

To investigate the role of MYBL2 homologous genes in anthocyanin biosynthesis within Brassica crops, we utilized existing RNA-seq data to analyze the expression patterns of these genes in the leaves (purple and green) of B. rapa, B. oleracea, and B. juncea, as well as in various tissues of B. napus, including different flower colors and purple petals at various developmental stages ( Figure 4 ; Supplementary Table S9 ). The results indicated that BraMYBL2.A02 was significantly up-regulated in the purple leaves of B. rapa ( Figure 4A ), BolMYBL2.C06 was significantly up-regulated in the purple leaves of B. oleracea ( Figure 4B ), and BjuMYBL2.B03 was significantly up-regulated in the purple leaves of B. juncea. Additionally, the expression levels of BjuMYBL2.A02 and BjuMYBL2.B05 in purple leaves were also found to be higher than those in green leaves ( Figure 4C ).

In B. napus, we first examined the expression patterns of four MYBL2 homologous genes across different tissues: siliques (green and purple), leaves (green and purple), stems (green and purple), petals (white and purple), and seed coats (yellow and black). Our findings revealed that the expression levels of these four MYBL2 homologous genes were higher in purple leaves compared to green leaves, with BnaMYBL2.A07 showing significant up-regulation in tissues exhibiting anthocyanin accumulation ( Figure 4D ). Moreover, BnaMYBL2.A07 and BnaMYBL2.C06 were highly expressed in petals of various colors, particularly in those with anthocyanin accumulation ( Figure 4E ; Supplementary Figure S2 ). Furthermore, we analyzed the expression patterns of the four MYBL2 homologous genes during the four developmental stages of purple petals. We found that BnaMYBL2.C06 was significantly highly expressed in stages S1 to S3, while BnaMYBL2.A07 exhibited high expression across all stages (S1 to S4), with the highest expression level observed in stage S2 ( Figure 4F ).

To further investigate the role of MYBL2 in the regulation of anthocyanin synthesis across various tissues of B. napus and to identify hub genes, we selected leaves (green and purple), stems (green and purple), and petals (white and purple) from B. napus. Six sample groups, each consisting of three biological replicates, were subjected to RNA sequencing, and the resulting phenotypic data were analyzed using WGCNA ( Figure 5 ).The results indicated that these genes were primarily clustered into three modules: MEblue, MEturquoise, and MEgrey ( Figure 5A ), with the MEblue module exhibiting the strongest gene interactions ( Figure 5B ). Subsequently, we conducted a correlation analysis among the three modules across the six sample groups, revealing that the MEturquoise module displayed the strongest correlation in purple petals ( Figure 5C ). Furthermore, we examined the regulatory relationships of ABGs within this module, identifying 62 enriched genes and complex regulatory interactions among them ( Figure 5D ). Notably, among the four MYBL2 homologous genes in B. napus, only BnaMYBL2.C06 was present in this module, and it exhibited regulatory relationships with 40 ABGs, further supporting the hypothesis that BnaMYBL2.C06 may function as a hub gene regulating anthocyanin synthesis.

To further validate the results of the RNA-seq analysis, we conducted qRT-PCR validation using leaves (green and purple), stems (green and purple), and petals (white and purple) of B. napus. We selected eight genes related to the anthocyanin biosynthetic pathway (BnaCHS.A02, BnaCHI.C04, BnaF3H.A09, BnaDFR.C09, BnaANS.C01, BnaUGT.A08, BnaTT8.C09, BnaPAP1.A07.b) and four MYBL2 homologous genes in B. napus for qRT-PCR analysis ( Figure 6 ). The results indicated that eight ABGs were significantly differentially expressed in purple leaves and petals, particularly the key genes BnaDFR.C09 and BnaANS.C01, which catalyze the conversion of flavonoids to anthocyanins, showing significant up-regulation in all three purple tissues. Notably, four homologous copies of MYBL2 exhibited significant differential expression in purple petals and purple leaves, but not in the purple and green stems. This expression pattern was analogous to that of the genes BnaCHS.A02, BnaCHI.C04, BnaF3H.A09, BnaANS.C01, and BnaUGT.A08. Importantly, the expression pattern of the four MYBL2 homologs aligns with that of the known positive regulator of anthocyanin biosynthesis, BnaPAP1.A07.b.

MYBL2 is a typical R3-MYB transcription regulator, belonging to the S4-derived R3-MYB subclass, characterized by a partial R2 domain, an EAR-repression motif, and a TLLLFR inhibitory domain (Khusnutdinov et al., 2021; LaFountain and Yuan, 2021). In this study, we systematically identified MYBL2 homologous genes in 31 Brsssicaceae species with published reference genome information, based on the AtMYBL2 sequence of A. thaliana, resulting in a total of 48 identified genes. Notably, no MYBL2 homologous gene was found in A. arabicum, which may be attributed to the loss of the MYBL2 gene during evolution or to incomplete genome assembly. Interestingly, the copy number of MYBL2 homologous genes within the same species closely aligns with the phylogeny and evolution of the Brassicaceae family (Murat et al., 2015; Oh and Dassanayake, 2019). For instance, there is only one MYBL2 homologous gene in A. lyrata, A. haleri, and other plants that have not undergone genome-wide doubling events, whereas there are two MYBL2 homologous genes in B. rapa, B. oleracea, and other species that have experienced such doubling events. Allotetraploid species, such as B. juncea and B. carinata, possess three MYBL2 homologous genes, while B. napus has four. The Ka/Ks analysis results of 48 MYBL2 homologous genes showed that only 10 genes were positively selected during the evolution process, these genes were mainly concentrated in species with more human selection, such as Brassica, while the other MYBL2 homologous genes were mainly purifying selected during the evolution process. Furthermore, phylogenetic and collinearity analyses of 16 MYBL2 homologous genes in Brassica crops provide further evidence supporting this observation, indicating that the distribution of these 16 MYBL2 homologous genes across the A, B, and C subgenomes originated from a common ancestor but differentiated at an earlier stage. The distribution of these genes on different chromosomes of the A, B, and C subgenomes reflects the orthologous relationship of these two genes. At the same time, conserved domain analysis found that the BjuMYBL2.A02 gene exhibits a deletion of the R3 domain, which indicates that the gene structure of MYBL2 has also changed during species evolution. However, the genomes of the allotetraploids B. juncea, B. carinata, and B. napus are currently at the telomere-to-telomere level due to the presence of high repetitive sequences, necessitating further investigation as subsequent genomic information becomes available.

The regulation of MYBL2 on anthocyanin biosynthesis in plantswas first reported in A. thaliana. It was discovered that MYBL2 work as a repressor that competes with PAP1/2 for binding to GL3/TT8, thereby inhibiting the formation of the PAP1/2-GL3/TT8-TTG1 transcriptional regulatory complex, which results in the negative regulation of anthocyanin biosynthesis in A. thaliana (Dubos et al., 2008; Matsui et al., 2008). Knockout of AtMYBL2 in A. thaliana, the expression levels of the key gene DFR and the regulatory gene TT8 were found to be up-regulated, resulting in a significant increase in anthocyanin content (Dubos et al., 2008). In addition, HY5 can inhibit MYBL2 activity by activating the expression of MYBD, which in turn increases anthocyanin accumulation in A. thaliana (Nguyen et al., 2015). And it has been observed that MYBL2 can interact with AP2 to form the AP2-MYBL2-TT8/EGL3 complex, which disrupts the formation of the MBW complex. This disruption inhibits the expression of ANR, TT12, TT19, and AHA, thereby hindering the biosynthesis of proanthocyanidins in A. thaliana seeds (Jiang, ). In Brassica crops, BrMYBL2.1 is considered to negatively regulate anthocyanin biosynthesis in Zicaitai (B. rapa) and Chinese cabbage (Zhang et al., 2020; Kim et al., 2022), while the anthocyanin content in B. oleracea increased following MYBL2 mutation (Khusnutdinov et al., 2022). These findings suggest that the function of MYBL2 in B. rapa and B. oleracea is consistent with that in A. thaliana, where it negatively regulates anthocyanin biosynthesis. However, in this study, we observed that the expression patterns of four MYBL2 homologous genes in the leaves, stems, and petals of B. napus aligned with the trend of ABGs, with expression levels significantly up-regulated in tissues exhibiting high anthocyanin content. This indicates a close relationship between MYBL2 and anthocyanin biosynthesis in B. napus. Moreover, we found that the expression pattern of MYBL2 was consistent with that of BnaPAP2.A7b, and their competing bHLH transcription factor TT8, along with the regulated structural gene DFR, also exhibited similar patterns. These results further suggest that as plant anthocyanins are synthesized and accumulated, the expression of MYBL2 is up-regulated, competing with BnaPAP2.A7b for TT8 to inhibit further anthocyanin synthesis. Recent studies have revealed that MYBL2 interacts with PAP1/2 to diminish its transcriptional activation activity, consequently reducing the expression of key genes DFR and TT19 involved in anthocyanin biosynthesis (Xing et al., 2024). This finding is particularly insightful for our subsequent research, as it elucidates the similar expression trends of MYBL2 and PAP1/2 observed in this study. We speculate that the substantial synthesis of anthocyanins triggers a feedback regulation mechanism in B.napus. This feedback may disrupt the formation of the MBW complex by enhancing MYBL2 expression and its interaction with PAP1/2, thereby inhibiting further synthesis and accumulation of anthocyanins. The underlying mechanisms warrant further investigation.