The draft coconut genome served as the reference for obtaining the coconut MYB gene sequences (Xiao et al., 2017). A Hidden Markov Model (HMM) profile representing the MYB DNA-binding domain with accession (PF00249) was acquired from the Pfam protein family database (http://pfam.xfam.org/) (Finn et al., 2016). Subsequently, this profile was employed as a query (P < 0.001) for the identification of all potential CnMYB genes. In parallel, AtMYB gene sequences were utilized as query sequences to detect MYB genes within the coconut genome. The identified MYB genes containing conserved domains underwent further analysis, while those lacking the PF00249 conserved domain were excluded from the dataset. Amino acid sequences of Arabidopsis MYB proteins (AtMYBs) were obtained from the Arabidopsis Information Resource (TAIR) database (https://www.arabidopsis.org/). The BLAST searches were conducted against the coconut genome databases using Arabidopsis MYB protein sequences as queries to identify CnMYB gene families (Katiyar et al., 2012). In total, 179 CnMYB genes were identified from the coconut genome. Additional information regarding the number of amino acids, molecular weight (MW) and isoelectric point (pI) of each MYB protein of interest was gathered using the ExPASy proteomic website (https://web.expasy.org/compute_pi/). The CELLO tool was used to predict the intracellular distribution of all CnMYB genes.

A neighbor-joining (NJ) phylogenetic analysis was conducted by MEGA-X based on the alignment. Bootstrap analysis with 1000 replicates was performed to calculate the reliability of the NJ tree (Kumar et al., 2018). The necessary data, including mRNA sequences, CDS and gene annotation summaries of coconut MYB genes, were retrieved from Gigascience Database 2017 (http://dx.doi.org/10.5524/100347) ( Supplementary Table 1 ). To confirm the structures of the coconut MYB genes, the mRNA sequences were aligned with the complete coconut genome sequence. The Gene Structure Display Server was then employed to analyze and determine the intron-exon organization of the coconut MYB genes. For the analysis of conserved sequence motifs within MYB proteins, the Motif Elicitation (MEME) tool (http://meme-suite.org/tools/meme) (Bailey et al., 2009) was utilized, and the results are presented in Supplementary Table 2 . Additionally, the TBtools software (version 1.045) developed by Chen et al. (2020) was used to visualize the chromosomal localization of MYB genes in the coconut genome. This analysis utilized the annotated genomic data from the coconut genome database to accurately map the MYB genes onto the corresponding chromosomes.

Samples from four distinct species of coconuts, namely Red Dwarf (RD), Yellow Dwarf (YD), Brown Dwarf (BD), and Green Dwarf (GD), were carefully selected for analysis. The specific parts examined included the leaf stalk epidermis (LSE), epicarp (EP), sepal (SE), and flower spike branch (FSB) of 7-month-old fruits ( Figure 1 ). Nine coconut fruits from the same fruit bunch were taken. Three coconut fruits were taken for collection of peel, sepals, floral branches and petioles separately and mixed them to form one biological replicate. In this way two different group were prepared to make another two biological replicate. These research samples were generously provided by the Coconut Research Institute (CRI), Chinese Academy of Tropical Agricultural Sciences (CATAS), located in Wenchang, Hainan, China. Following collection, all samples were immediately cryopreserved using liquid nitrogen and subsequently transferred to a freezer set at -80°C for future use.

The RNA extraction was carried out by the MRIP method (Method for RNA isolation from Palm) (Xiao et al., 2012) and has been improved. The protocol of the RNA extraction method according to the Iqbal et al. (2019). The quality of the extracted RNA (i.e., degradation and contamination) was assessed using 1% agarose gels. Additionally, the integrity of the RNA was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States), while the concentration was measured using a Nanodrop Spectrophotometer (IMPLEN, CA, United States). Sequencing of the samples was performed on the BGISEQ-MGI2000 platform instrument at BGI Genomics (Shenzhen 518083, China) with three biological replicates for each sample. To ensure high-quality data, the raw sequencing reads were processed to obtain clean reads. This involved excluding reads with adapters, reads with unknown bases more than 5% and low-quality base ratios more than 20% were filtered using SOAPnuke (version 1.4.0) (Chen et al., 2018). The resulting clean reads were stored in FASTQ format. Subsequently, the data were aligned to the reference genome using HISAT (v2.1.0) (Kim et al., 2015), and then matched with the assembled unique genes using Bowtie2 (v2.2.5) (Langmead and Salzberg, 2012). RNA-seq by expectation maximization (RSEM) (version 1.2.8) was utilized to calculate the expression levels of the genes (Li and Dewey, 2011). For functional annotation, the assembled unigenes were annotated using databases such as KEGG and GO, and transcription factors were predicted as well. Differential gene analysis within groups was conducted using DESeq with the conditions of Fold Change ≥ 2 and adjusted P-value ≤ 0.001 (Wang et al., 2010).

To conduct a comprehensive analysis of coconut MYB genes, transcriptoe datasets from various coconut tissues, including LSE, EP, SE and FSB, were utilized. The expression levels of genes were quantified using Reads Per Kilobase Million (RPKM) values, which were further Log2 transformed to facilitate comparative analysis. To visualize the expression patterns of 12 R2R3-MYB genes, a heatmap was generated using TBtools.

Forward and reverse primers for qPCR analysis were designed by National Center for Biotechnology Information (NCBI) Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast) with melting temperatures of 55-60°C, primer length 19-25bp, GC content 50-60% and amplicon size of 80-200bp. To prevent amplification of non-target gDNA, the primers were designed to span intronic regions. The properties of each primer were evaluated using the PCR Primer Stats software. The primers used for qPCR analysis can be found in Supplementary Table 4 . The reference gene CnACTIN was used as an internal control (Xia et al., 2014).The Quick and Reliable RNA Extraction Method (QRREM) was employed to isolate total RNA from the epicarp of coconut fruit, following the protocol of Iqbal et al. (2019). Subsequently, the isolated RNA underwent quality and quantity assessment through agarose gel electrophoresis and Nanodrop spectrophotometer analysis. For reverse transcription, 1 μg of RNA was used with the MightyScript first-strand cDNA synthesis kit following the manufacturer’s instructions. The qPCR reactions were conducted using the 2 × SYBR Green qPCR ProMix in 96-well optical plates on a Mastercycler ep realplex⁴ machine. The qPCR reactions were carried out with a total reaction volume of 10 μL, consisting of amplification at 95°C for 5 s, 55°C for 15 s, and 68°C for 20 s. The melting stage involved heating from 60°C to 95°C for 20 min. Each experiment was conducted with biological and technical triplicates. The fold change in expression level for each sample was calculated by normalizing the CT value relative to a reference gene, using the 2^-ΔΔCt method (Livak and Schmittgen, 2001).

The experiments were conducted in triplicate (n=3) to ensure the reliability of the findings. Mean values, accompanied by standard errors of the mean, were used to present the data. SAS software (SAS Inc., Cary, NC, USA) was employed to perform ANOVA (analysis of variance) and DMRT (Duncan’s Test) with a significance level set at p < 0.05 and p < 0.01 to assess the significance among various treatments. Correlation coefficients were calculated based on the mean values.

A comprehensive analysis of MYB genes in the coconut genome revealed a total of 179 MYB genes, which wereeffectively profiled from the coconut genome. According to the numbering in the coconut genome, it is arranged from small to large and positioned as CnMYB01 to CnMYB179. Detailed information about each of these profiled MYB genes as shown in Supplementary Tables 1 and 2 . The length of the peptide chains in the CnMYB proteins varied from 61 to 1150 amino acids, as shown in Supplementary Table 1 . Additionally, the predicted proteins exhibited a molecular weight range of 7.8 to 128.1 and isoelectric points ranging from 4.09 to 10.11, as indicated in Supplementary Table 1 . Regarding the intracellular distribution of the CnMYB putative proteins the vast majority was predicted to be localized in the nucleus. Based on the analysis of the N-terminal aa region, 154 CnMYBs were predicted to be localized in the nucleus whereas CnMYB02 would be localized in the cytoplasm and CnMYB03 and CnMYB132 both in the cytoplasm and nucleus ( Supplementary Table 1 ).

Through a search using DNA sequence annotations, we observed that 179 identified coconut MYB genes were distributed across 16 chromosomes. Among these genes, 148 genes exhibited uneven distribution patterns across the chromosomes ( Figure 2 ). More specifically, our analysis revealed that 11 CnMYB genes were located on chromosome 1, 15 genes on chromosome 2, 11 genes on chromosome 3, 12 genes on chromosome 4, 10 genes on chromosome 5, 12 genes on chromosome 6, 11 genes on chromosome 7, 10 genes on chromosome 8, 12 on chromosome 9 held, 8 genes on chromosome 10 held, 4 genes on chromosome 11, 9 genes on chromosome 12, 13 genes on chromosome 13, 11 genes on chromosome 14, 2 genes on chromosome 15, and 3 genes on chromosome 16 ( Figure 2 ). Additionally, 29 CnMYB genes were assigned to chromosomes that could not be determined. Notably, chromosome 2 contained the highest number of CnMYB genes (15), while chromosome 4, 6 and 9 each contain 12 CnMYB genes. In contrast, chromosome 15 had the lowest number of MYB genes (2) ( Figure 2 ).

The exon-intron organization of the 179 CnMYB genes was examined using the Gene Structure Display Server program, as illustrated in Figure 3 . The analysis revealed that the majority MYB genes exhibited varying numbers of introns, ranging from 0 to 15 ( Figure 3 ). Among the 179 CnMYB genes, CnMYB166 had the highest number of introns (15), followed by CnMYB80 and CnMYB117 with 11 introns each. Additionally, CnMYB81 and CnMYB151 contained 10 introns, while CnMYB130 had 9 introns. Furthermore, CnMYB12 had 7 introns, and CnMYB145 had 6 introns. Notably, 28 genes, including CnMYB1, CnMYB3, CnMYB4, CnMYB15, CnMYB25, CnMYB31, CnMYB32, CnMYB34, CnMYB44, CnMYB46, CnMYB63, CnMYB65, CnMYB73, CnMYB74, CnMYB85, CnMYB86, CnMYB87, CnMYB93, CnMYB97, CnMYB100, CnMYB101, CnMYB125, CnMYB139, CnMYB143, CnMYB158, CnMYB173, CnMYB174, and CnMYB177, did not contain any introns. Furthermore, three genes (CnMYB80, CnMYB117 and CnMYB166) had the maximum number of exons (≥12). Most of the MYB genes exhibited a small number of introns, typically ranging from 0 to 3, suggesting a conserved pattern of intron distribution within the MYB genes. The disparity observed in the CnMYB gene composition suggests noteworthy deviation within the coconut genome. Additionally, we performed protein motif profiling to assess the variation in the CnMYB gene family of coconut. Our findings revealed the presence of 10 preserved motifs across all 179 identified MYB genes, as depicted in Figures 4 , 5 and Supplementary Table 3 .

To construct a phylogeny of MYB proteins, a maximum likelihood (ML) method was employed, using 179 MYB proteins sequences from coconut and 154 MYB proteins sequences from Arabidopsis thaliana ( Figure 6 ). The resulting dendrogram classified the MYB genes into 34 distinct subgroups, denoted as S1-S26, C1-C6, 3R-MYB, and 4R-MYB ( Figure 6 ). These subgroups represent four main types of MYB protein families: 4R-MYB, 3R-MYB, R2R3-MYB, and the coconut-specific subgroups ( Figure 6 ). Out of the 34 subgroups, 29 included proteins from both coconut and Arabidopsis, while the remaining five were specific to either coconut (S9, S12 and S15) or Arabidopsis (C1 and C2). Similar species-specific subgroups of MYBs have been found in other plant species, such as Solanum tuberosum (Li Y. et al., 2020), Casuarina equisetifolia (Wang et al., 2021) and Petunia (Chen et al., 2021). Within the 3R-MYB subfamily, we identified four CnMYB genes (CnMYB18, CnMYB81, CnMYB130, CnMYB139) and ten AtMYB genes. Similarly, the 4R-MYB subfamily consisted of four CnMYB genes (CnMYB70, CnMYB71, CnMYB151, CnMYB166) and four AtMYB genes. The C1 to C6 subfamily included forty-seven CnMYB genes and thirteen AtMYB genes. The remaining 124 CnMYB genes and 127 AtMYB genes were associated with the R2R3-MYB family. The R2R3-MYB family was further divided into subgroups S1-S26, with varying numbers of CnMYB and AtMYB genes. Notably, subgroups S3 and S11 comprised only one CnMYB member each, making them the smallest groups, while subgroup S14 comprised 15 members, making it the largest group. Interestingly, no CnMYB genes were found in subgroups S9, S12, or S15, suggesting potential gene loss during coconut genome evolution or acquisition in Arabidopsis evolution. The presence of more CnMYB genes than AtMYB genes in certain subgroups indicates functional distinction of MYB genes among various plant species. These phylogenetic findings suggest that MYB genes clustered within the same set may share preserved functions, which should be further investigated through experimental approaches.

We conducted an analysis of the expression profiles of 179 CnMYB genes in four distinct tissues of coconut, including LSE, EP, SE and FSB using RNA-seq data from a database. The results were visualized using a heatmap, which revealed the expression patterns of 12 R2R3-MYB genes in coconut tissues ( Figure 7 , Supplementary Table 5 ). Among these genes, three (CnMYB15, CnMYB68, CnMYB101) exhibited the highest expression across all four coconut tissues ( Figure 7 ). Conversely, genes with the lowest expression levels were predominantly observed in LSE, specifically CnMYB94 and CnMYB164 ( Figure 7 ). Notably, genes with higher expression levels were primarily observed in the EP tissues of all four coconut species ( Figure 7 ). Furthermore, the expression of CnMYB8 was found to be higher in all four tissues of the Green Dwarf (GD) species compared to other coconut varieties.

Moreover, we performed qPCR analysis to investigate the expression patterns of 12 candidate MYB genes (CnMYB8, CnMYB15, CnMYB27, CnMYB28, CnMYB61, CnMYB63, CnMYB68, CnMYB94, CnMYB101, CnMYB150, CnMYB153, CnMYB164) during different stages of fruit development (4M, 7M, and 10M). Among these genes, CnMYB28, CnMYB68, and CnMYB101 displayed higher expression levels compared to the other genes, especially CnMYB68 ( Figure 8 ). Notably, the expression level of CnMYB68 in the 10-month-old fruit peel of Red Dwarf (RD), Yellow Dwarf (YD), and Brown Dwarf (BD) coconuts was significantly higher than in other stages and Green Dwarf (GD) coconut (P<0.05). Additionally, the expression levels of CnMYB101 in RD, YD, and BD were higher than in GD, with the highest expression level observed in the 4-month-old coconut peel of RD, YD, and BD ( Figure 8 ). Conversely, the expressions of CnMYB8, CnMYB15, CnMYB27, CnMYB63, CnMYB150, CnMYB153, and CnMYB164 were relatively lower compared to the other genes.

MYB transcription factors have been extensively studied and recognized as one of the largest families of transcription factors in plants. Among all MYB subgroups, the R2R3-MYB protein comprises the highest number of members that can play a crucial role in various aspects of the secondary metabolism in plants (Song et al., 2021; Zhang et al., 2021). The R2R3-MYB family can control the diverse biological processes such as plant metabolism, growth and development (Li et al., 2023; Zhang et al., 2023). Additionally, it regulates plant responses to hormones and various types of stresses (Zhao et al., 2023). Several genomic-scale investigations have been concluded to identify MYB family members in various plant genomes (Li Y. et al., 2020; Liu et al., 2021; Zhao et al., 2023). The genomic identification of R2R3-MYB transcription factors has been conducted in various plant species, revealing their significant presence in plant genomes. For instance, extensive sequencing efforts have led to the identification of 138 R2R3-AtMYBs in Arabidopsis thaliana (Katiyar et al., 2012), 99 R2R3-DcMYB in Dendrobium catenatum (Zhang et al., 2021), 393 R2R3-TaMYBs in Triticum aestivum (Wei et al., 2020), 111 R2R3-StMYB in Solanum tuberosum (Li Y. et al., 2020) and 174 R2R3-NtMYB in Nicotiana tabacum (Yang et al., 2022). However, little is known about the MYB gene family in Cocos nucifera.

In this study, we identified 179 MYB members in Cocos nucifera through bioinformatics analysis, and characterized their phylogenetic relationships with A. thaliana MYBs. Our compresinsive analysis encompassed phylogenetic characteristics, physical and chemical properties, gene structure, chromosome location, tissue-specific expression patterns, and expression profiles during coconut fruit development stages. Our findings revealed that the MYB domain and other motifs of CnMYBs exhibit a high degree of conservation with those of Arabidopsis AtMYBs. Nonetheless, we also observed divergence between CnMYBs and AtMYBs, indicating a combination of conservation and diversity within plant MYBs.

Specifically focusing on the R2R3-MYB family, our analysis identified 124 members within the coconut genome. Phylogenetic analysis categorized these members into 26 subfamilies, with unified subfamily members displaying similar motifs and gene structures ( Figure 6 ). The oil palm genome was found to contain 152 R2R3-MYB family members, which were classified into 23 subfamilies (Zhou et al., 2020). Furthermore, the protein sequences of the 124 coconut R2R3-MYB family members exhitited varying lengths, ranging from 89 (CnMYB123) to 555 (CnMYB118) amino acids, with an average length of 314 amino acids. Their molecular masses also varied, spanning from 10.37 (CnMYB123) to 59.79 kDa (CnMYB118), with an average molecular mass of 34.84 kDa. Theoretical isoelectric points ranged from 4.28 (CnMYB90) to 10.06 (CnMYB107), with an average molecular mass of 6.99 ( Supplementary Table 1 ). In a separate investigation involving 124 longan R2R3-MYB family members, their protein sequences ranged from 134 to 661 amino acids, with an average length of 307 amino acids. The molecular masses ranged from 15.69 to 775.81 kDa, with an average of 34.68 kDa and theoretical isoelectric points ranged from 4.91 to 10.6 (Lv et al., 2023). Based on the findings of a previous study conducted by Zhang et al. 2021, our study provides further support for the hypothesis that the MYB gene family exhibits a significant level of evolutionary conservation among various plant species.

The pattern of gene structure is a useful tool for studying the evolutionary associations within a gene family. In our study, we identified 124 R2R3-CnMYB genes, which exhibited a range of exon numbers from 1 to 16 ( Figure 3 ). Most of the R2R3-CnMYBs, similar to those found in other plant species, consisted of three exons and two introns (Liu et al., 2014; Chen et al., 2022; Du et al., 2022). It is worth noting that within the same subfamily, the exon/intron patterns of R2R3- CnMYB genes showed remarkable, with the majority of genes having no more than two introns. This finding aligns with previous research that also observed the presence of a maximum of two introns in most plant R2R3-MYB genes (Sabir et al., 2022; Yin et al., 2022). Furthermore, we analyzed the motif compositions of CnMYB proteins. Our findings revealed that the majority of MYB proteins contained motifs 1, 2, 3, 4, 5, and 6 ( Figure 4 ), while most MYB proteins contained motifs 1, 3, 5, and 8 in Dendrobium catenatum (Zhang et al., 2021). The majority of MYB genes, included several motifs (1, 2, 3, 4, and 5) in Prunus avium (Sabir et al., 2022). These motifs were conserved within specific subgroups, and proteins in the same subgroup that share these motifs likely (Qin et al., 2023).

Prior research has implicated R2R3-MYB genes in essential functions in the anthocyanins, carotenoids and flavonoid biosynthesis in a variety of plants. For instance, in red-fleshed apple, MdMYB22 and MdMYB12 have been identified as key regulatorsof proanthocyanidins and flavonols biosynthesis (Wang et al., 2017). Moreover, PpMYB10 has been found to primarily control anthocyanin biosynthesis in the exocarp of Prunus persica (Ravaglia et al., 2013), while MdMYB1 is a major regulator of anthocyanin biosynthesis in red-skinned fruit (Zhang et al., 2019). In pear fruit, PbMYB12b has been found to positively regulate flavonol biosynthesis by enhancing the expression of PbCHSb and PbFLS (Zhai et al., 2019). In wolfberry, Lba11g0183 and Lba02g01219 have been identified as candidate genes involved in carotenoid biosynthesis (Yin et al., 2022). Additionally, in Elaeis guineensis, the VIR gene, encoding R2R3-MYB-like transcription factor with homology to Lilium LhMYB12 and similarity to Arabidopsis PRODUCTION OF ANTHOCYANIN PIGMENT1 (PAP1), regulate the heterogeneity of red and yellow fruit color (Singh et al., 2014). Genetic variations in the VIRESCENS gene have been associated with the conspicuousness of fruit colors in palms, indicating potential selection by frugivores (Wang et al., 2022). Notably, when comparing the VIRESCENS gene sequences of oil palm and date palms, no matching R2R3-MYB gene was found in coconut. In coconut, the VIRESCENS gene appears to be disrupted by the insertion of a highly repetitive sequence spanning 100 kb (Wang et al., 2022).

The influence of light exposure time on flavonoid and anthocyanin biosynthesis in fruits has been extensively studied by Premathilake et al. (2020). It was found that exposing pear fruit to light for a long duration up-regulated the expression of the R2R3-MYB DNA-binding protein PpMYB17 (Premathilake et al., 2020), resulting in higher anthocyanin biosynthesis (Alabd et al., 2022). In our study, we observed that the epicarp covered by sepals exhibited a lighter or even whiter color, which may be attributed to the obstruction of other coconut fruits within the clusters, which can result in variations in the overall pericarp color. It is speculated that the uneven coloration is influenced by light obstruction ( Supplementary Figure 1 ). We also identified higher expression levels of CnMYB68, CnMYB101, and CnMYB28 in various tissues and developmental stages across the four coconut species, as evident from both transcriptome data and quantitative data. However, further experiments are required to confirm their functional roles.