Integrative Genome-Wide Analysis of R2R3-MYB Transcription Factors in Oryza sativa subsp. japonica Reveals Candidate Regulators of Anthocyanin Biosynthesis

1 Introduction

Transcription factors (TFs) are regulators of gene expression, acting as molecular switches that modulate the transcriptional system in response to developmental and environmental signals [1]. Among the diverse families of TFs, the MYB (myeloblastosis) superfamily represents one of the largest and most functionally diverse classes in plants. Initially discovered in avian myeloblastosis virus, the MYB domain is defined by a conserved helix-turn-helix (HTH) DNA-binding motif, typically composed of up to three imperfect repeats (R1, R2, and R3) of ~52 amino acids each [2].

Plant MYB genes are classified into four groups based on repeat numbers, namely R1-MYB, R2R3-MYB, R1R2R3-MYB (3R-MYB), and deviant MYBs. Among these, R2R3-MYB TFs are the most abundant and plant-specific, regulating diverse processes such as metabolism, stress tolerance, hormone signaling, and development [3]. Importantly, R2R3-MYB TFs regulate phenylpropanoid and flavonoid biosynthesis, including anthocyanins, by activating structural genes such as Chalcone Synthase (CHS), Dihydroflavonol 4-Reductase (DFR), and Anthocyanidin Synthase (ANS) [4,5]. This regulation is mediated through the MYB–bHLH–WD40 (MBW) complex. Within this complex, the MYB component confers DNA-binding specificity by recognizing conserved cis-regulatory motifs in the promoters of target genes, thereby determining the transcriptional activation of anthocyanin biosynthetic enzymes. The formation and activity of the MBW complex are also influenced by upstream regulatory signals, as MYB TFs often bind promoter regions enriched with light- and stress-responsive cis-elements that mediate environmental regulation of pigment biosynthesis. Such regulatory interactions allow plants to modulate anthocyanin accumulation in response to external stimuli, including light intensity and abiotic stresses, which are known to induce flavonoid pathway genes. Moreover, anthocyanin production is frequently restricted to specific tissues or developmental stages, reflecting spatially regulated expression patterns of MYB TFs that determine where MBW complexes are assembled and activated. For example, functional characterization of R2R3-MYB regulators in Rosa rugosa demonstrated that RrMYB12 and RrMYB111 directly modulate the transcription of flavonoid biosynthetic genes, leading to coordinated accumulation of flavonols and anthocyanins in pigmented tissues [4]. Therefore, identifying candidate anthocyanin regulators requires integrating phylogenetic relationships with promoter cis-element architecture and tissue-specific expression profiles, enabling a multi-layered evaluation of MYB genes that are both evolutionarily related to known regulators and transcriptionally active under relevant regulatory contexts. Flavonoids and anthocyanins not only contribute to pigmentation but also provide photoprotection, pathogen defense, and antioxidative functions [4]. Given their central regulatory roles, R2R3-MYB TFs were selected as the primary focus of the present investigation, leading to the identification and characterization of MYB TFs of rice, which are potentially involved in anthocyanin biosynthesis.

In Arabidopsis thaliana, R2R3-MYB TFs such as AtMYB75 (PAP1), AtMYB90 (PAP2), AtMYB113, AtMYB15, and AtMYB114 have been well-characterized as key activators of anthocyanin biosynthesis [6]. AtMYB15 overexpression upregulated genes in the shikimate pathway, a precursor route for phenylpropanoids, possibly enhancing anthocyanin biosynthesis under wounding and stress conditions [7]. These TFs coordinate with bHLH proteins (TT8, GL3) and the WD-repeat protein TTG1 to direct hierarchical activation of biosynthetic genes. Similar regulatory modules have been identified in monocots such as Zea mays, with species-specific adaptations and gene duplications that reflect evolutionary divergence [8].

Rice (Oryza sativa), the staple food for over half the global population, contains a large number of MYB genes, many yet uncharacterized. Genome-wide analyses have identified a large MYB TF family in rice, with over 100 MYB-like sequences reported and approximately 70–90 belonging to the R2R3-MYB subgroup [9,10,11,12]. Comparative studies across plant species identified 197 MYB genes in rice and 198 in Arabidopsis, classifying them into multiple subfamilies based on conserved domains and motif composition and suggesting that gene duplication events contributed to both conserved and species-specific expansion patterns [9,10]. Subsequent analyses of 124 rice R2R3-MYB genes revealed diverse expression responses to abiotic stresses, including drought, cold, and salinity, indicating functional diversification within the family [11,12]. Despite these advances in characterizing rice MYB genes, comprehensive phylogenetic and functional analyses of MYBs associated with pigment biosynthesis remain limited. Nevertheless, specific regulators such as OsMYB3 have been shown to activate anthocyanin biosynthesis by upregulating OsDFR and OsANS in black rice pericarp, highlighting the regulatory potential of MYB TFs in pigmentation pathways [13]. The deep regulatory network of MYBs influencing pigment biosynthesis and their orthologous relationships with well-studied Arabidopsis MYBs need further investigation. Studying pigment production in rice and identifying novel MYB genes is vital for improving grain quality and nutritional value, especially in South Asia, where rice is a staple for over 70% of the population [14]. Pigmented rice varieties rich in anthocyanins and flavonoids offer enhanced antioxidant properties and reduced risk of chronic diseases [15].

Genome-wide analyses enable systematic identification of gene families and provide insights into their evolution, functions, and stress-related roles [9,10,11,12]. With a compact 389 Mb genome and high-quality annotation, O. sativa serves as a model for monocots and cereal crops. The japonica cultivar “Nipponbare” was the first rice genome to be sequenced [16], and it’s well-assembled. Japonica rice also holds considerable economic and regional importance, being widely cultivated in East Asia and temperate regions, where it contributes significantly to food security. Functionally, japonica varieties exhibit traits such as cold tolerance and unique grain quality, making them valuable for breeding and comparative genomics. Previous studies on Oryza sativa MYB TFs have primarily focused on genome-wide identification and basic classification of MYB family members, providing valuable catalogs of these regulators. However, many of these analyses did not systematically integrate evidence such as phylogenetic proximity to known anthocyanin-regulating MYBs, promoter cis-regulatory architecture, and expression patterns in pigmented tissues to prioritize candidates potentially involved in anthocyanin biosynthesis [17]. Consequently, the functional relevance of many MYB genes in pigment regulation remains unclear. To address this gap, the present study applies a comprehensive genome-wide analysis of MYB TFs in Oryza sativa subsp. japonica, with a particular focus on the R2R3-MYB subgroup. The analysis included structural classification, chromosomal distribution, gene duplication events, conserved motif composition, and phylogenetic relationships. Because the MYB superfamily in rice is large and functionally diverse, such a genome-wide investigation is necessary to systematically identify candidate regulators potentially associated with anthocyanin biosynthesis.

To infer potential functions, phylogenetic relationships were examined through comparative analysis with MYB TFs from Arabidopsis thaliana. Candidate pigment-related MYBs were further predicted by integrating multiple lines of evidence, including clustering with known anthocyanin regulators, promoter cis-element analysis, and gene expression profiling. Although several MYB regulators have previously been reported in rice, a unified analysis combining phylogenetic relationships with expression-based candidate prioritization and known anthocyanin activators has not been performed. The integrative framework used in this study therefore provides a systematic approach for identifying candidate MYB regulators involved in flavonoid and anthocyanin biosynthesis. Overall, this work identifies a set of preliminary candidate MYB genes for future functional validation and establishes a foundation for understanding transcriptional regulation of pigment biosynthesis in rice, with potential applications in crop improvement and nutritional enhancement.

2 Materials and Methods

2.1 Database Search and Sequence Retrieval

To identify MYB TFs, the conserved MYB DNA-binding domain (PF00249) was used as the primary query motif. The amino acid sequence of the MYB transcription factor AtMYB15 (Gene ID: At3g23250) from Arabidopsis thaliana was initially retrieved from the UniProt database (https://www.uniprot.org/) to confirm the presence of the MYB-like DNA-binding motif using the Motif Finder tool available at the KEGG GenomeNet platform (https://www.genome.jp/tools/motif/). The identified PF00249 domain was subsequently used to construct a hidden Markov model (HMM) profile for genome-wide searches. The presence of PF00249 in well-characterized anthocyanin regulators, such as PAP1, was also confirmed to ensure conservation among functionally validated MYBs. Genome-wide identification of MYB proteins was performed using the HMMER v3.3 package against protein datasets obtained from Phytozome v13 (https://phytozome-next.jgi.doe.gov/) for Oryza sativa subsp. japonica (cv. Nipponbare) and Arabidopsis thaliana and from Ensembl Plants (https://plants.ensembl.org/) for Oryza rufipogon. Searches were conducted using the PF00249 HMM profile with an E-value cutoff of 1e−5 to identify candidate MYB proteins. Redundant sequences and truncated entries were removed to retain non-redundant protein sets. Candidate sequences were then validated using multiple conserved domain databases, including Pfam and the NCBI Conserved Domain Database (CDD). Sequences lacking a complete MYB DNA-binding domain were excluded from further analysis. Using this pipeline, a total of 109, 104, and 108 MYB proteins containing the PF00249 domain were identified from Oryza sativa subsp. japonica, Arabidopsis thaliana, and Oryza rufipogon, respectively. While initial retrieval included all MYB-domain-containing proteins, downstream analyses focused on candidate R2R3-MYBs, which are the predominant regulators of anthocyanin biosynthesis.

2.2 Multiple Sequence Alignment and Phylogenetic Analysis

Multiple sequence alignment of MYB protein amino acid sequences was performed using MUSCLE v2.1. Poorly aligned regions and positions containing excessive gaps were removed using the partial deletion option implemented in MEGA X to improve alignment quality prior to phylogenetic reconstruction. Phylogenetic analysis was conducted using the Maximum Likelihood (ML) method. The Jones–Taylor–Thornton (JTT) amino acid substitution model was used for tree construction. Rate variation among sites was assumed to follow uniform rates across all positions. MYB TFs contain highly conserved R2R3 DNA-binding domains; the alignment primarily consisted of homologous regions with limited sequence length variation after the removal of poorly aligned sites. Under these conditions, assuming uniform substitution rates provides a stable framework for reconstructing higher-level clade relationships within large gene families. The ML heuristic search was performed using the Nearest Neighbor Interchange (NNI) method to explore tree topology. Branch support was evaluated using bootstrap analysis with 1000 replicates, and bootstrap values ≥80% were considered to indicate well-supported clades. A total of 321 MYB protein sequences containing the PF00249 domain were included in the initial phylogenetic analysis, comprising 109 sequences from Oryza sativa subsp. japonica, 104 sequences from Arabidopsis thaliana, and 108 sequences from Oryza rufipogon.

To facilitate functional inference, the phylogenetic dataset was expanded to include previously characterized rice MYB TFs and well-established anthocyanin-activating MYBs from other plant species. In total, 109 MYB proteins from Oryza sativa subsp. japonica were analyzed together with 20 functionally characterized rice MYB TFs and representative anthocyanin-related MYBs reported in other plant species. The inclusion of these reference proteins enabled the identification of rice MYB genes clustering with known pigment regulators. Candidate anthocyanin-related MYBs were inferred based on their phylogenetic proximity to experimentally validated anthocyanin activators and supported by bootstrap values ≥70%. Because the objective of the phylogenetic analysis in this study was to identify clustering patterns with known anthocyanin regulators rather than to estimate evolutionary rates or divergence times, the uniform-rate model was retained.

2.3 Physicochemical Properties and Subcellular Localization Prediction

Information on Oryza sativa MYB genes was retrieved from the Phytozome (https://phytozome-next.jgi.doe.gov/) and ProtParam databases (https://web.expasy.org/protparam/). Phytozome served as the primary source for genomic features, including chromosome number, gene location, and orientation, as well as the lengths of mRNA and corresponding peptides. In contrast, the ProtParam tool was used to predict various physicochemical properties of the MYB proteins, such as molecular weight, isoelectric point (pI), instability index, and Grand Average of Hydropathicity (GRAVY) values (Table S1). Additionally, subcellular localization prediction of the MYB proteins was assessed using the Wolf PSORT server (https://wolfpsort.hgc.jp/), which analyzes protein sequences to predict potential cellular compartments along with associated confidence scores. This analysis aimed to determine the intracellular distribution of MYB proteins, supporting functional characterization [10].

2.4 Conserved Motifs and Domain Analysis of OsMYB Genes

Domain annotation of MYB proteins from Oryza sativa subsp. japonica was performed using the NCBI Conserved Domain Database (CDD) (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) to confirm the presence of the MYB DNA-binding domain and associated conserved regions. The identified domains were further verified and visualized using TBtools. To identify conserved sequence motifs within the OsMYB protein sequences, motif discovery was performed using the MEME Suite (version 5.5.9) (https://web.expasy.org/protparam/). Protein sequences were analyzed in classic discovery mode with the site distribution set to Zero or One Occurrence Per Sequence (ZOOPS). The maximum number of motifs was limited to 10, and the motif width was restricted to a range of 6–50 amino acids. The significance threshold for motif detection was set to a maximum E-value of 0.05, and a first-order Markov background model generated from the input sequences was used to estimate amino acid frequencies. Motif discovery was carried out using the expectation–maximization algorithm implemented in MEME. The resulting motif sequences and motif logos were extracted to characterize conserved patterns among OsMYB proteins. Motif composition and distribution were further visualized by integrating MEME outputs with phylogenetic tree data using TBtools. Motif 1 and Motif 2 corresponded to the conserved MYB DNA-binding domain (R2/R3 repeats). The PF00249 MYB-like DNA-binding domain of OsMYB was manually scanned for the presence of a bHLH-interaction motif. The sequence of the R3 repeat was examined for the standard pattern [DE]Lx2 [RK]x3Lx6Lx3R, with conserved leucine (L) and arginine/lysine (R/K) residues identified as potential interaction sites.

2.5 Intron-Exon Structure Analysis of OsMYB Genes

The genomic DNA and coding sequences (CDS) of 109 putative OsMYB genes were obtained from the Phytozome database. These sequences were then analyzed using the Gene Structure Display Server 2.0 (GSDS; http://gsds.cbi.pku.edu.cn/) to investigate the exon–intron organization within each gene [10].

2.6 Cis-Regulatory Elements of OsMYB Genes

Promoter sequences of OsMYB genes from Oryza sativa subsp. japonica were retrieved from the Phytozome database. For each gene, a 1000 bp region upstream of the translation start site (ATG) was extracted and considered as the putative promoter region. These promoter sequences were analyzed using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html) to identify all potential cis-regulatory elements present within the sequences. Following the initial identification of cis-acting regulatory elements, the analysis specifically focused on motifs associated with flavonoid and anthocyanin biosynthesis. Particular attention was given to regulatory elements such as MYB-binding sites (MBS), light-responsive elements including ACE, Chs-CMA1a, Chs-CMA2a, and G-box, and hormone- or stress-responsive elements such as ABRE, which are known to influence flavonoid pathway gene expression [18]. Subsequently, the distribution of these anthocyanin-related cis-elements was examined in the promoter regions of the newly identified R2R3 OsMYB genes to assess their potential regulatory involvement in pigment biosynthesis. Cis-element distributions and positional patterns were visualized using TBtools. This analysis was conducted as a descriptive survey of regulatory motifs, and no statistical enrichment testing against genomic background sequences was performed.

2.7 The miRNA Analysis for OsMYB Genes

Potential miRNA interactions with OsMYB genes from Oryza sativa subsp. japonica were investigated to explore post-transcriptional regulatory mechanisms. Known rice miRNA sequences were retrieved from the PmiREN database, and coding sequences (CDS) of the OsMYB genes were used as query sequences for target prediction. miRNA target prediction was performed using the psRNATarget web server with default scoring parameters. Predicted miRNA–mRNA interactions were filtered based on expectation score and target accessibility criteria provided by the psRNATarget algorithm to reduce false-positive predictions. Following the initial prediction, analysis focused specifically on miRNAs targeting the newly identified R2R3 OsMYB TFs, which were prioritized as potential regulators of anthocyanin biosynthesis. Reported functional roles of these miRNAs were obtained from literature surveys, with emphasis on their known or potential involvement in pigment biosynthesis or integration with MBW complex pathways. The predicted miRNA–target interactions were compiled and summarized to support candidate prioritization of R2R3 OsMYBs potentially regulated by miRNAs in anthocyanin-related pathways (Table S2).

2.8 Gene Ontology of OsMYB Genes

Gene Ontology (GO) analysis was performed to investigate the functional roles of OsMYB genes from Oryza sativa subsp. japonica. GO annotations associated with molecular function and biological process categories were initially retrieved from the Phytozome database. The identified OsMYB gene set (109 proteins) was submitted to ShinyGO v0.741 (http://bioinformatics.sdstate.edu/go/) to conduct GO enrichment analysis. During enrichment analysis, the complete annotated gene set of Oryza sativa subsp. japonica available in ShinyGO was used as background reference. Statistical significance of GO term enrichment was evaluated using a hypergeometric test, and p-values were adjusted for multiple comparisons using the Benjamini–Hochberg false discovery rate (FDR) correction method. GO terms with an adjusted FDR value ≤ 0.05 were considered significantly enriched. To quantify enrichment significance, fold enrichment values were calculated and expressed as −log10(FDR) scores. These values were used to rank enriched GO categories and visualize the relative significance of each functional group. The enriched GO terms were summarized to identify major functional roles of OsMYB TFs related to anthocyanin biosynthesis.

2.9 Protein–Protein Interactions

To further investigate potential protein–protein interactions of candidate OsMYB TFs, interaction network analysis was performed using the STRING database (v0.74111). In addition to OsMYB proteins, known bHLH and WD40 components implicated in anthocyanin biosynthesis in Oryza sativa OsB2, OsRc, and OsTTG1 were included in the analysis. This allowed the visualization and prediction of interaction networks specifically involving MYB–bHLH–WD40 (MBW) complexes, which are central to the regulation of anthocyanin structural genes. Predicted interactions were examined to identify candidate OsMYBs that may associate with these MBW components, providing evidence for their potential regulatory role in anthocyanin biosynthesis.

2.10 Gene Duplication and Synteny Analysis

The analysis based on non-synonymous substitution rate (Ka) and synonymous substitution rate (Ks) values was used to estimate the divergence time of the MYB gene family in Oryza sativa. The Ka/Ks ratios for each pair of paralogous genes were calculated using the TBtools software. These ratios helped determine the molecular evolutionary rates of gene pairs. The divergence time (T) was calculated using the formula T = Ks/(2λ), where the neutral substitution rate (λ) was assumed to be 6.5 × 10−9 substitutions per site per year [19].

Gene duplication events within the O. sativa MYB gene family were identified using the MCScanX toolkit with default settings, implemented through TBtools, which allowed for the detection of both segmental and tandem duplications. To visualize the genomic relationships of paralogous MYB genes, micro-synteny analysis was conducted using the Micro Synteny View function in TBtools. Furthermore, dual syntenic comparisons were performed between O. sativa and A. thaliana and O. sativa and O. rufipogon to assess the conservation and divergence of MYB genes across species.

2.11 Expression Profiling of OsMYB Genes in Rice Pericarps

Transcriptome data (GSE67987) were retrieved from the NCBI Gene Expression Omnibus (GEO) database to examine the expression of OsMYB TFs associated with anthocyanin biosynthesis. Expression values of 109 OsMYB genes were obtained from pericarp tissues of black rice near-isogenic lines and mixed rice genotypes at 7 and 14 days after heading. Fold-change (FC) values were calculated to assess developmental regulation in black rice [FC (14 vs. 7 DAF, Black)] and genotype-specific differences [FC (Black vs. Mixed, 7 DAF) and FC (Black vs. Mixed, 14 DAF)]. Genes showing FC > 1 across all three conditions were considered as candidate upregulated regulators. Heatmaps based on fold-change values were generated using TBtools to visualize expression patterns across developmental stages and genotypes, where color intensity (white to dark red) represents increasing fold-change levels. Notably, these observations are based solely on fold-change thresholds without statistical testing, and the mixed rice genotype group serves as an approximate comparator rather than a clearly defined control; therefore, these results were interpreted as preliminary signals for candidate gene prioritization rather than definitive evidence of differential expression.

3 Results

3.1 Comprehensive Genome-Wide and Phylogenetic Characterization of OsMYB Genes

A phylogenetic analysis using protein sequences from O. sativa, O. rufipogon, and A. thaliana clustered the MYB proteins into three major clades (Fig. 1). In Clade I, several putative orthologous groups were identified. AtMYB70 clustered with OsMYB18 (bootstrap 1.0000), and AtMYB59 clustered with OsMYB36 (bootstrap 0.9820). The Arabidopsis proteins AtMYB53, AtMYB81, and AtMYB60 formed a subcluster (bootstrap 0.9920 and 0.9681), which then grouped with OsMYB94 (bootstrap 0.9401). AtMYB31 clustered with OsMYB98 (bootstrap 0.9142), and AtMYB83 grouped with OsMYB50 and OsMYB44, where the OsMYB44/OsMYB50 subcluster had bootstraps of 1.0000 and 0.9920, the combined OsMYB cluster 0.8563, and the addition of AtMYB83 0.8623. Clade II (bootstrap 0.9960) contained the AtMYB63/OsMYB61 pair. Clade III (bootstrap 0.8683) encompassed multiple putative orthologs, including AtMYB92/OsMYB75, AtMYB37/OsMYB7, AtMYB86/OsMYB26, AtMYB62/OsMYB12/OsMYB62, AtMYB50/AtMYB51/OsMYB1, AtMYB102 with OsMYB35/OrMYB46 (bootstrap 0.8523), AtMYB39 with OsMYB9/OrMYB43 (bootstrap 0.9960), and AtMYB57/AtMYB101 with OsMYB54/OsMYB87 (bootstrap 0.9960) and AtMYB55/OsMYB25 (bootstrap 0.9501). Bootstrap values across the tree ranged from 0.85 to 1.00, indicating strong support for these clades and subclades. Phylogenetic analysis further supported the potential role of several newly identified candidate MYB genes in anthocyanin-related regulatory networks through downstream analysis. Notably, candidates such as OsMYB44, OsMYB50, OsMYB61, OsMYB65, and OsMYB94 clustered in close proximity to Arabidopsis MYBs, including Arabidopsis thaliana MYB31, MYB83, MYB63, MYB60, and MYB53, suggesting possible functional relationships. It is important to note that orthology is suggested based on sequence similarity and clustering, and functional conservation should not be inferred without experimental validation. This analysis provides a comprehensive overview of the OsMYB gene family and their phylogenetic relationships, highlighting potential evolutionary patterns among rice and Arabidopsis. The subcellular localization data indicates distinct patterns of presence of 109 OsMYB proteins across various cellular compartments. The majority of OsMYB proteins were predominantly localized in the nucleus (Fig. 2).

Figure 1: The phylogenetic relationships among MYB genes from O. sativa, A. thaliana, and O. rufipogon were investigated using a maximum likelihood (ML) approach with 1000 bootstrap replicates. A total of 109 MYB gene sequences were analyzed using MEGA X v12 software. The resulting phylogenetic tree revealed three distinct clades, each visually represented by a unique color. MYB genes orthologous from A. thaliana and O. sativa were clearly distinguished by green and red markers, respectively.

Figure 2: Subcellular localization prediction of OsMYB proteins revealed their predominant presence in the nucleus, peroxisomes, and cytoplasm. The majority of proteins were represented by red-colored markers, indicating their primary localization within these cellular compartments.

3.2 Gene Structure Analysis

Gene structure analysis revealed substantial variation in exon–intron organization among the 109 OsMYB genes (Fig. 3). The observed exon/intron configurations included 1/0, 2/1, 2/2, 2/3, 3/2, 3/3, 4/0, 6/5, 7/6, and 10/9. Among these, the most common structure was three exons and two introns, followed by two exons and one intron, indicating a general conservation of gene architecture within the family alongside notable diversification.

Focusing on the newly identified candidate genes associated with anthocyanin regulation (OsMYB5, OsMYB6, OsMYB7, OsMYB8, OsMYB9, OsMYB15, OsMYB24, OsMYB25, OsMYB32, OsMYB41, OsMYB43, OsMYB44, OsMYB45, OsMYB47, OsMYB50, OsMYB59, OsMYB60, OsMYB61, OsMYB62, OsMYB65, OsMYB68, OsMYB69, OsMYB70, OsMYB71, OsMYB81, OsMYB86, OsMYB89, OsMYB94, OsMYB98, OsMYB101, OsMYB104, and OsMYB105), a clear pattern emerged. The majority of these candidates exhibited the 3 exon/2 intron (2/3) or 2 exon/2 intron (2/2) structures, which are characteristic features commonly observed in R2R3-MYB genes implicated in flavonoid and anthocyanin biosynthesis in plants. This structural conservation supports their potential functional relevance. In contrast, a few candidate genes displayed more divergent architectures. For example, OsMYB47 showed the most complex structure with 10 exons and 9 introns, while OsMYB71 (6/5) also exhibited relatively higher exon–intron numbers. Conversely, simpler structures such as 2/1 (e.g., OsMYB32, OsMYB44, OsMYB60, OsMYB70, OsMYB81, and OsMYB89) and 1/0 (OsMYB50) were also observed, indicating structural variability even within candidate regulators. Gene structures alone is insufficient to infer function; therefore, candidate gene identification in this study was based on an integrative approach combining phylogenetic relationships, promoter cis-element analysis, conserved motif detection, and expression profiling.

Figure 3: The exon-intron structures were depicted using shapes colored in Yellow and black to represent the exons contained within the genes.

3.3 Conserved Motif Analysis and Domain Prediction Analysis of OsMYBs

Motif analysis of the MYB gene family in Oryza sativa confirmed the presence of conserved domains and motifs characteristic of MYB TFs (Fig. 4). The identified MYB genes possessed several key domains, including MYB_DNA-binding, MYB_DNA-binding superfamily, MYB_DNA-bind_6, the P_C superfamily, and the SANT superfamily. Specifically, Motif 1 (associated with MYB_DNA-binding) and Motif 2 (MYB_DNA-binding superfamily) were the most frequently occurring domains across the MYB gene family suggesting a predominant representation of R2R3-type MYB TFs (Fig. 4). These domains, often located toward the N-terminal region, are key components of the R2 and R3 repeats responsible for DNA binding in anthocyanin biosynthesis. For instance, AtMYB75 from Arabidopsis harbor similarly conserved R2R3 motifs and are experimentally validated to activate anthocyanin biosynthesis [20], thus reinforcing the functional relevance of these motifs.

In addition to these, five motifs (Motifs 6, 7, 8, 9, and 10) were conserved across several MYB proteins, indicating shared functional or structural roles. Interestingly, OsMYB70 and OsMYB33 were the only genes that harbored the SANT superfamily domain, while OsMYB54 uniquely contained the P_C superfamily domain. The presence of additional domains such as the SANT and P_c superfamily in a few proteins further indicates domain expansion and potential specialization in regulatory roles. Domains unique to these, such as SANT or P_C, possibly point to particular regulatory functions for them. In other TFs, variation in motifs located in the C-terminal regions was found to affect the specificity of secondary metabolic gene targets, including flavonoid biosynthetic genes. Domains in OsMYBs that have expanded might indicate possible functional divergence for these to be verified.

Figure 4: The analysis of OsMYB genes involved the identification of conserved domains and motifs, alongside the construction of their phylogenetic relationships. Using MEME version 5.5.2, ten distinct motifs were identified and visualized through a color-coded bar graph. This motif distribution was aligned with a phylogenetic tree, allowing for a comparative view of evolutionary relatedness among the proteins.

3.4 Cis-Element Analysis of OsMYBs

Cis-element analysis of the promoter regions of the 109 OsMYB genes identified 105 regulatory elements belonging to several functional categories, including 27 light-responsive elements, 22 stress-responsive elements, 20 development- and metabolism-related elements, 16 hormone-responsive elements, 6 core promoter/general transcription elements, 1 signal transduction/calcium/redox element, and 13 uncharacterized elements (Table 1; Fig. 5). Among these, core promoter elements such as the CAAT-box and TATA-box were the most abundant, consistent with their fundamental role in transcription initiation. Moderate frequencies of regulatory motifs such as AT-TATA-box, ABRE, G-box, MYC, and MYB elements were also detected across multiple OsMYB promoters.

Particular attention was given to cis-elements associated with flavonoid and anthocyanin pathway regulation in the promoters of the prioritized candidate MYB genes. Several light-responsive motifs, including the MYB-recognition element (MRE), Chs-CMA1a, Chs-CMA2a, G-box, and Box 4, were identified in these promoters [21]. These elements are known to mediate transcriptional activation of flavonoid biosynthetic genes under light-regulated condition. In addition, development-related elements such as AC-I and AC-II, which are associated with tissue-specific expression and flavonoid biosynthesis, were detected in several candidate OsMYB promoters. The presence of these regulatory motifs suggests that the prioritized MYB genes may participate in the transcriptional regulation of anthocyanin biosynthesis through integration of environmental and developmental signals.

Figure 5: The investigation of cis-regulatory elements within the promoter regions of O. sativa MYB genes revealed their potential roles in a wide range of plant developmental and physiological processes. This analysis detailed the presence and frequency of various regulatory elements associated with growth, stress response, hormone signaling, and environmental cues. Furthermore, the spatial distribution of these elements across the promoter regions was assessed, offering insights into their organization.

3.5 Protein–Protein Interaction

To explore potential regulatory roles of candidate OsMYBs in anthocyanin biosynthesis, protein–protein interaction (PPI) networks were analyzed using STRING, including the key MBW components in O. sativa: OsB2, OsRc, and OsTTG1 (Fig. 6). The resulting network comprised 95 nodes and 682 edges, with an average node degree of 14.4 and an average local clustering coefficient of 0.639, indicating substantial connectivity among the proteins. The expected number of edges based on random chance was 46, and the PPI enrichment p-value was <1.0e−16, confirming that the observed interactions are significantly more interconnected than expected by chance.

Within the network, several candidate OsMYBs were predicted to interact directly with MBW components. Specifically, OsTTG1 interacted with OsMYB25, OsMYB44, OsMYB65, OsMYB104, and OsMYB105, while OsRc was predicted to interact with OsMYB15 and OsMYB44. These predicted interactions suggest that the corresponding OsMYBs may potentially function as components of the MYB–bHLH–WD40 (MBW) complex, which is known to regulate anthocyanin biosynthesis in plants. The network analysis also indicated clusters of OsMYBs with multiple putative interactions, highlighting them as promising candidates for involvement in tissue-specific pigment regulation. However, these interactions are based on computational predictions and should be interpreted with caution. Notably, these OsMYBs were independently identified as candidate regulators through integrative phylogenetic and expression analyses, providing additional support for their possible roles in pigment biosynthesis.

Figure 6: Protein–protein interaction network of OsMYB TFs involved in anthocyanin biosynthesis in Oryza sativa. The network illustrates interactions between OsMYBs and the known bHLH and WD40 components, including OsB2, OsRc, and OsTTG1, highlighting connectivity and potential regulatory relationships in the anthocyanin pathway.

3.6 Gene Duplication and Synteny Analysis of OsMYB

The expected chromosomal positions indicated that the OsMYB genes were widely distributed across all 12 chromosomes of O. sativa, with additional loci present on unassembled (ChrUn) and synthetic regions (ChrSy), as shown in Fig. 7. A notable concentration of OsMYB genes was observed on chromosome 1 and chromosome 4. Several candidate anthocyanin-associated OsMYB genes showed chromosomal clustering patterns. For instance, OsMYB5, OsMYB6, OsMYB7, OsMYB8, and OsMYB9 were all located on chromosome 8, while another group including OsMYB60, OsMYB62, OsMYB68, OsMYB70, and OsMYB71 was located on chromosome 2. Such local clustering of candidate MYB genes may reflect historical tandem or segmental duplication events that contributed to the expansion of regulatory MYB TFs in rice.

Figure 7: The chromosomal distribution of the OsMYB gene family in O. sativa. The figure’s red bars represent O. sativa chromosomes, while the black lines indicate the functional relationships or co-regulation between OsMYB genes. This analysis provides insights into the spatial arrangement and potential interactions among OsMYB genes within the O. sativa genome.

The circular heatmap presents the estimated divergence times of paralogous O. sativa MYB gene pairs, calculated using synonymous substitution rates (Ks) and the molecular clock formula T = Ks/(2λ), where λ is set at 6.5 × 10−9 substitutions per site per year, a rate commonly used for monocot plants [19] (Fig. 8). The color gradient from white to deep red represents increasing divergence time in millions of years ago (MYA), with darker shades indicating more ancient duplication events. Several gene pairs, such as OsMYB34/OsMYB39, OsMYB2/OsMYB94, and OsMYB15/OsMYB39, exhibit high Ks values and are shaded in deep red, corresponding to estimated divergence times of approximately 400–500 MYA under the assumed substitution rate, likely during ancient whole-genome duplication events. Notably, some of these duplicated genes correspond to candidate anthocyanin regulators identified in this study. For example, both OsMYB15 and OsMYB94, which displayed relatively high Ks values in duplication analysis, were also prioritized as potential anthocyanin regulatory candidates based on integrative phylogenetic, promoter, and expression analyses. In contrast, many other gene pairs show much lower Ks values, indicating more recent duplications possibly resulting from segmental or tandem duplication events within the rice genome. A Ka/Ks ratio 1.0 for paralogous OsMYB gene pairs originated over 400–500 MYA indicates that these genes are likely evolving neutrally, with no strong evolutionary pressure to maintain or change their protein-coding sequences. Paralogous OsMYB gene pairs with a Ka/Ks ratio of 0 are under strong purifying selection, indicating evolutionary pressure to maintain their protein sequence and function. This observation suggests that certain duplication events within the MYB family may have contributed to the diversification of TFs potentially involved in pigment-related pathways, although functional validation would be required to confirm these roles [70].

Figure 8: The time of gene duplication for various paralogous OsMYB gene pairs in O. sativa based on synonymous (Ks) and nonsynonymous (Ka) substitution rates. The analyses were performed using the TBtools software. The Ka/Ks ratio, representing the rate of nonsynonymous (Ka) to synonymous (Ks) substitutions, was calculated to infer the selection pressure acting on these gene pairs.

Several OsMYB genes in rice exhibit strong syntenic conservation with genes in Arabidopsis, indicating evolutionary retention of key regulatory functions (Fig. 9a,b). Both OsMYB5 and OsMYB7 were syntenically linked with Chr5 genes, likely arising from segmental duplications and retaining conserved domains seen in Arabidopsis MYBs [71]. OsMYB31 and OsMYB50, both on Chr1, also showed syntenic conservation, hinting at their descent from ancestral MYB genes. OsMYB50 identified in this study as a candidate regulator of anthocyanin biosynthesis shares conserved genomic context, supporting its potential functional relevance. In Chr3, OsMYB54 and OsMYB60 exhibited multiple syntenic links, suggesting stable evolutionary roles. OsMYB60, along with OsMYB62, OsMYB68, and OsMYB71 on chromosome 2, were among the newly proposed candidate anthocyanin regulators and showed conserved syntenic relationships with regions on chromosomes 3 and 5. These patterns suggest that some of the candidate genes may have retained ancestral regulatory functions following duplication events. Gene duplication and synteny analysis was mainly performed in order to place the OsMYB family into an evolutionary context.

Figure 9: Dual synteny plots were generated to illustrate the genomic collinearity of putative OsMYB genes between O. sativa and two other plant species. Red lines indicate orthologous gene pairs of OsMYB genes across species. Chromosomes of O. sativa are shown in orange for clear distinction. (a) Synteny map between O. sativa and A. thaliana, (b) Synteny map between O. sativa and O. rufipogon.

3.7 Osa-miRNA–OsMYB Regulatory Interactions in O. sativa

Predicted Osa-miRNA interactions were analyzed to identify potential post-transcriptional regulators of the newly prioritized R2R3 OsMYB genes implicated in anthocyanin biosynthesis (Table S2). Several miRNA families were found to target candidate anthocyanin-related OsMYBs, suggesting potential roles in modulating their expression. For example, Osa-miR159a–e was predicted to target OsMYBs associated with reproductive tissue expression, whereas Osa-miR164a, b, d, and e were linked to OsMYBs expressed in leaves and developing grains [72]. Similarly, Osa-miR166a, b, d, f, g, and h were predicted to interact with OsMYBs showing tissue-specific expression patterns consistent with anthocyanin accumulation. Several additional miRNAs, including Osa-miR319a/b, Osa-miR396a–c, Osa-miR408, Osa-miR528, Osa-miR5505, and Osa-miR5809b, were identified as potential regulators of candidate OsMYBs expressed under environmental and developmental cues known to influence pigment accumulation [73,74,75,76]. Development-related miRNAs such as Osa-miR171a/h, Osa-miR172c, Osa-miR444a/d, and Osa-miR531b were also predicted to target OsMYBs with spatial or temporal expression patterns relevant to anthocyanin deposition in reproductive tissues and grains.

Overall, these predicted miRNA–candidate OsMYB interactions highlight specific miRNA families that may contribute to post-transcriptional regulation of anthocyanin biosynthesis, either by fine-tuning expression in pigmented tissues or integrating developmental and environmental signals. These interactions are based on computational predictions and literature-based functional annotations, and experimental validation is necessary to confirm their direct regulatory roles in anthocyanin accumulation.

3.8 Gene Ontology

GO enrichment analysis of Oryza sativa MYB proteins identified several biological processes associated with transcriptional regulation and metabolic control that are relevant to secondary metabolite biosynthesis. Significant enrichment was observed for regulation of transcription, DNA-templated (3 genes, fold enrichment = 21.7, FDR = 2.0 × 10−3), regulation of RNA biosynthetic process (3 genes, fold enrichment = 21.7, FDR = 2.0 × 10−3), and regulation of RNA metabolic process (3 genes, fold enrichment = 21.1, FDR = 2.0 × 10−3), indicating that a subset of OsMYB proteins may function in transcriptional control of downstream metabolic pathways (Fig. 10). Additional enriched categories included regulation of biosynthetic processes (3 genes, fold enrichment = 19.7, FDR = 2.0 × 10−3) and regulation of metabolic process (3 genes, fold enrichment = 13.9, FDR = 3.3 × 10−3), which are consistent with regulatory roles in plant secondary metabolism. Since MYB TFs are widely recognized regulators of genes involved in the Anthocyanin biosynthesis pathway and the related Flavonoid biosynthesis pathway, these enriched processes suggest that some OsMYB proteins may participate in transcriptional networks controlling pigment biosynthesis [77]. In addition, developmental processes such as anther development (1 gene, fold enrichment = 216.2, FDR = 1.2 × 10−2) and flower development (1 gene, fold enrichment = 34.9, FDR = 4.2 × 10−2) were also enriched, which may be relevant because anthocyanin pigments commonly accumulate in reproductive tissues. Overall, the enrichment pattern supports a potential regulatory association between OsMYB TFs and pathways involved in secondary metabolism and pigment-related biological processes.

Figure 10: The bar plot illustrates significantly enriched Gene Ontology (GO) categories of 109 OsMYB proteins. This highlights key molecular functions and biological pathways that MYB genes are likely involved in to understand the functional roles of these transcription factors in rice. The numerical values of each bar represent the fold enrichment scores expressed as −log₁₀ of the false discovery rate (FDR), indicating the statistical significance of enrichment for each GO term. Higher log₁₀(FDR) values indicate stronger statistical significance, meaning that the associated GO category is very unlikely to have appeared by random chance.

3.9 Functional Interpretation of OsMYB Genes Based on Arabidopsis Orthologues

The functional interpretation of the OsMYB gene family based on Gene Ontology annotations and A. thaliana orthologues reveals their diverse regulatory roles in rice (Table S3). Functional interpretation of the candidate OsMYB genes using orthologues from Arabidopsis thaliana suggests that several of these TFs may be associated with the regulation of flavonoid metabolism, including pathways linked to anthocyanin accumulation. Among the identified candidates, OsMYB44 shows the strongest association with anthocyanin regulation. Its Arabidopsis orthologue, AtMYB112, has been reported to positively regulate anthocyanin biosynthesis and enhance tolerance to abiotic stresses such as salinity and drought by activating antioxidant and stress-responsive genes [78]. This functional similarity indicates that OsMYB44 may contribute to anthocyanin accumulation in rice as part of a broader protective response to environmental stress.

OsMYB45 also represents a potential regulator of flavonoid-derived metabolites. Its orthologue AtMYB60-related TFs participate in jasmonate-mediated developmental processes and promote flavonol biosynthesis through activation of the FLS1 gene [79]. Because flavonols and anthocyanins share common upstream biosynthetic steps within the flavonoid pathway, this relationship suggests that OsMYB45 may indirectly influence anthocyanin biosynthesis by regulating flavonoid metabolic flux. In addition, several candidate genes including OsMYB15, OsMYB61, and OsMYB101 are orthologous to AtMYB97-like regulators that control secondary metabolite pathways in Arabidopsis [80]. Although these TFs are primarily associated with glucosinolate metabolism and stress-responsive signaling, MYB proteins involved in secondary metabolism frequently exhibit functional diversification across plant species. Therefore, their identification through expression profiling, phylogenetic relationships, and cis-regulatory element analysis suggests that these OsMYBs may participate in broader phenylpropanoid or flavonoid metabolic networks that include anthocyanin biosynthesis.

Taken together, orthologue-based functional inference indicates that OsMYB44 and OsMYB45 are the most plausible anthocyanin-associated regulators among the candidate genes, while OsMYB15, OsMYB61, and OsMYB101 may play supportive or indirect roles in secondary metabolite regulation. These findings are consistent with the phylogenetic placement and expression patterns of these genes, supporting their potential involvement in the transcriptional networks governing anthocyanin accumulation in rice.

3.10 Identification of Novel R2R3 OsMYBs as Putative Anthocyanin Regulators

A phylogenetic analysis was carried out of 109 OsMYBs along with the 20 rice MYB TFs documented to date [81] and well-characterized anthocyanin-activating MYBs from other species (Fig. 11). Based on phylogenetic clustering, several OsMYBs were identified as candidates for involvement in pigment biosynthesis due to their orthology with functionally validated anthocyanin regulators from other species. OsMYB36 clustered with DcMYB6, VlMYBA2, CmMYB6, and PhAN2, all known to activate anthocyanin pathways (bootstrap = 0.7924). Similarly, OsMYB54 showed strong orthology to SbY1 (Sorghum bicolor) (bootstrap = 0.9501), and OsMYB31 showed strong orthology to BvMYB12 (Beta vulgaris) (bootstrap = 0.8543). The domain analysis also confirms these as R2R3 OsMYB TFs. Conversely, OsMYBs that were most closely related to the three previously characterized rice pigment regulators were excluded to prioritize novel candidates.

Figure 11: Circular Phylogenetic Tree of 109 O. sativa MYB Proteins Integrated with Known Pigment-Related MYB Regulators from Rice and Other Plant Species. Red arrows indicate literature-reported O. sativa MYB genes involved in pigmentation; green arrows highlight well-characterized anthocyanin-activating MYB regulators from other species. Literature-reported OsMYBs are labeled with their original gene IDs, while similarly identified gene IDs in our dataset are renamed in the ‘OsMYB.1’ format to avoid ID overlapping.

Fold-change analysis of transcriptome data identified several OsMYB TFs showing consistent upregulation during pericarp development and higher expression in black rice compared with mixed rice genotypes (Fig. 12). Using the selection criterion of FC > 1 for developmental induction [FC (14 vs. 7 DAF, Black)] and genotype comparisons [FC (Black vs. Mixed, 7 DAF) and FC (Black vs. Mixed, 14 DAF)], a total of 32 OsMYB genes were identified as candidate upregulated regulators. Among these, several genes displayed strong expression differences, including OsMYB8 (FC 9.428, 1.404, and 2.627), OsMYB24 (FC 1.929, 2.296, and 4.345), OsMYB45 (FC 3.043, 3.404, and 1.747), OsMYB60 (FC 3.203, 1.931, and 1.044), OsMYB61 (FC 1.328, 4.148, and 4.821), OsMYB65 (FC 1.324, 27.526, and 16.021), and OsMYB81 (FC 2.440, 7.766, and 2.907), corresponding to FC (14 vs. 7 DAF, Black), FC (Black vs. Mixed, 7 DAF), and FC (Black vs. Mixed, 14 DAF), respectively, indicating strong developmental and genotype-associated upregulation. Notably, OsMYB65 exhibited the highest genotype-specific expression, with FC values of 27.53 and 16.02 in Black vs. Mixed rice at 7 and 14 DAF, respectively.

Figure 12: Expression Profiling of OsMYB Genes Associated with Anthocyanin Biosynthesis in Developing Rice Pericarps (GSE67987). Heatmaps illustrate fold-change (FC) patterns of 109 OsMYB TF genes in pericarp tissues of black rice near-isogenic lines and mixed rice genotypes 7 and 14 days after heading. Fold changes were calculated for developmental regulation in black rice [FC (14 vs. 7 DAF, Black)] and genotype comparisons [FC (Black vs. Mixed, 7 DAF) and FC (Black vs. Mixed, 14 DAF)] using expression data retrieved from the NCBI GEO dataset GSE67987. Color intensity from white to dark red represents increasing fold-change values across developmental stages and genotypic comparisons.

The predicted bHLH-interaction motif in these OsMYBs was located within the R3 repeat of the PF00249 MYB-like DNA-binding domain, which was present in all of these OsMYBs. It spans roughly the residues EDEQLR...CRLL...HRR, where conserved leucines (L) and arginines (R) were key for interaction. This short motif likely mediates binding to bHLH TFs in rice. Promoter analysis of these potential OsMYB anthocyanin candidates revealed several cis-regulatory elements associated with flavonoid and anthocyanin biosynthesis (Table 2). Notably, the MYB recognition element (MRE), a key motif involved in MYB-mediated regulation of flavonoid pathway genes, was identified in the promoters of OsMYB44 and OsMYB98, suggesting their potential role in regulating structural genes of the anthocyanin biosynthetic pathway. In addition, chalcone synthase-associated motifs, including Chs-CMA1a and Chs-CMA2a, which are known to control light-induced expression of flavonoid biosynthetic genes, were detected in OsMYB59, OsMYB69, OsMYB70, OsMYB101, and OsMYB62, indicating possible regulation of early anthocyanin pathway enzymes such as chalcone synthase. The presence of AC-I elements, associated with phenylpropanoid metabolism and flavonoid biosynthesis, was observed in OsMYB69 and OsMYB70, further supporting their potential involvement in pigment-related secondary metabolism. Additionally, MYB binding motifs and MYB-like sequences were widely distributed across promoters of genes such as OsMYB5, OsMYB6, OsMYB25, OsMYB43, and OsMYB81, suggesting potential autoregulatory or cross-regulatory interactions within the MYB TF network. Together, the occurrence of these anthocyanin-related cis-elements in specific OsMYB promoters indicates that these genes may participate in the transcriptional regulation of flavonoid and anthocyanin biosynthesis in rice, serving as preliminary candidate regulators for further experimental validation.

4 Discussion

This study aimed to identify candidate MYB TFs potentially involved in anthocyanin regulation in rice through an integrative genome-wide analysis. A total of 109 OsMYB genes were identified in Oryza sativa using domain-based screening. Previous studies have reported larger MYB families in rice, reflecting differences in search strategies and classification criteria [9,10,11,12]. By focusing on MYB proteins containing the conserved PF00249 domain associated with anthocyanin regulators such as AtMYB15 and PAP1 (AtMYB75), our analysis targeted MYB members with potential regulatory roles in pigment biosynthesis. This targeted approach is particularly relevant for exploring untapped pigmented rice germplasm in South Asia, including Sri Lankan traditional varieties that exhibit diverse grain pigmentation [82,83,84]. Identifying MYB regulators from such unexploited germplasm would enable targeted breeding of nutritionally superior cultivars, aligning with growing consumer demand for functional foods and dietary diversification [85].

Phylogenetic relationships and conserved motif analysis indicate that several OsMYBs cluster with experimentally validated anthocyanin regulators from other plant species [81]. Anthocyanin biosynthesis in plants is typically controlled by the MYB–bHLH–WD40 (MBW) regulatory complex, which activates transcription of late flavonoid biosynthetic genes [4]. Several newly identified candidate OsMYBs showed predicted interactions with key MBW components, including the bHLH proteins OsB2 and OsRc, and the WD40 protein OsTTG1. This suggests that they may function as regulatory hubs with the capacity to modulate MBW complex assembly or stability; such network centrality, when considered alongside their conserved domain architecture and expression patterns, raises the possibility that these factors contribute to the fine-tuning of spatiotemporal anthocyanin biosynthesis rather than acting solely as primary transcriptional activators. The presence of conserved R2R3 DNA-binding motifs and similarity to characterized MYBs suggests that several identified OsMYBs may participate in similar regulatory modules controlling pigment accumulation in rice grains. Gene structure and domain organization further support the regulatory potential of these TFs. R2R3-MYB genes involved in anthocyanin biosynthesis across plant species commonly possess relatively simple exon–intron architectures and conserved DNA-binding domains that enable interaction with promoters of structural genes in the flavonoid pathway [86]. The conservation of these features among several OsMYBs identified in this study suggests evolutionary preservation of transcriptional mechanisms controlling flavonoid metabolism.

Promoter analysis suggests that many OsMYB genes may integrate environmental and developmental signals relevant to anthocyanin accumulation [87]. Light-responsive elements were particularly abundant, consistent with the well-established role of light signaling in activating anthocyanin biosynthesis [21,30]. The presence of additional stress- and hormone-related cis-elements suggests that pigment accumulation may be coordinated with broader physiological responses, as anthocyanins frequently function as protective antioxidants during environmental stress [35]. However, it should be noted that cis-element predictions are based solely on sequence motifs and do not necessarily indicate functional regulatory activity. Experimental validation such as promoter-reporter essays would therefore be required to confirm the regulatory relevance of these elements [88].

Expression profiling provided an additional layer of evidence for prioritizing candidate regulators. Differential expressions during grain pigmentation enabled the identification of 32 OsMYB genes potentially associated with anthocyanin accumulation. Among these, genes such as OsMYB65 displayed strong expression patterns in pigmented grains, suggesting potential roles in activating flavonoid biosynthetic pathways. The prioritization of candidates based on combined phylogenetic relationships, regulatory motifs, and expression patterns strengthens the likelihood that these TFs participate in pigment regulation. Physicochemical characterization showed that most of these candidate proteins possess negative GRAVY values, indicating hydrophilic properties typical of TFs functioning in the nucleus [89] (Table S1). Their predicted isoelectric points (pI) varied widely, suggesting structural diversity that may influence DNA-binding affinity and interactions with regulatory complexes controlling secondary metabolism.

Comparative analysis with functionally characterized MYB genes highlights the limited number of experimentally confirmed pigment regulators currently known in rice. To date, only a few MYBs including OsC1, OsMYB3/Kala3, and OsPL/OsMYB55 have been demonstrated to regulate grain pigmentation [81]. The identification of additional candidate regulators in this study therefore expands the potential regulatory network underlying anthocyanin biosynthesis in rice and provides targets for future functional characterization.

Predicted interactions with several Osa-miRNA families further suggest an additional regulatory layer controlling MYB activity. miRNA-mediated regulation has been shown in other plant systems to fine-tune TFs involved in secondary metabolism, including flavonoid biosynthesis [90,91,92,93]. These predictions indicate that post-transcriptional control mechanisms may contribute to the spatial or environmental regulation of pigment accumulation, although experimental validation will be necessary to confirm these interactions.

Despite the integrative nature of the analysis, several methodological limitations should be considered. First, promoter cis-element analysis is descriptive and does not demonstrate functional TF binding. Second, the expression dataset used for candidate prioritization represents specific developmental stages and environmental conditions, which may not capture all regulatory contexts in which anthocyanin-related MYBs are active. Some candidate genes may therefore exhibit condition-specific or tissue-specific expression patterns that were not detected in the available data. In addition, the fold-change analysis provides an initial prioritization of candidate OsMYB genes that may regulate anthocyanin accumulation in black rice pericarp. However, due to the absence of replicates and statistical testing, these results cannot be interpreted as definitive evidence of differential expression. Future studies with biological replicates and experimental validation are needed to confirm the roles of these candidates in flavonoid and anthocyanin biosynthesis. Finally, phylogenetic similarity alone cannot establish functional equivalence between MYB orthologs across species.

Future studies should therefore focus on functional validation of prioritized candidates. Approaches such as transgenic overexpression, gene knockout using CRISPR/Cas systems, promoter activation assays, and co-expression analysis with structural genes of the flavonoid pathway (e.g., DFR, ANS, and UFGT) would help determine whether these MYBs directly regulate anthocyanin biosynthesis [94]. Validation of predicted miRNA–MYB interactions through degradome sequencing or reporter assays could further clarify post-transcriptional regulatory mechanisms [95].

Overall, this integrative genome-wide analysis identifies a prioritized set of MYB TFs potentially involved in anthocyanin regulation in rice. These findings provide a foundation for future functional studies and may support the development of pigment-rich rice varieties with enhanced nutritional and antioxidant properties.

5 Conclusions

This study provides a genome-wide characterization of 109 OsMYB TFs in rice, identifying candidate genes potentially involved in anthocyanin biosynthesis and pericarp pigmentation. Phylogenetic analysis, conserved R2R3-MYB domains, and expression profiling highlight genes that may act as key regulators, providing a framework for future functional studies to validate their roles. The analysis of miRNA involved cis-regulatory elements and comparative orthology with Arabidopsis and other rice species further supports their potential involvement in pigment accumulation and stress-related pathways. These candidates OsMYBs represent promising targets for future genetic improvement and biofortification efforts, although their functional activity and regulatory impact on anthocyanin biosynthesis require experimental confirmation. Overall, this work lays a foundation for targeted studies to elucidate regulatory networks controlling pigmentation and may ultimately contribute to the development of nutritionally enriched and stress-resilient rice varieties.