Next Article in Journal
Vitamin D and Vitamin D-Binding Protein in Health and Disease
Previous Article in Journal
Aluminum and Fluoride Stresses Altered Organic Acid and Secondary Metabolism in Tea (Camellia sinensis) Plants: Influences on Plant Tolerance, Tea Quality and Safety
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 24
Issue 5
10.3390/ijms24054641
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide Analysis of MYB Transcription Factors and Screening of MYBs Involved in the Red Color Formation in Rhododendron delavayi

by
Fenfang Long
,
Hairong Wu
,
Huie Li
*,
Weiwei Zuo
and
Qian Ao
College of Agriculture, Guizhou University, Guiyang 550025, China
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2023, 24(5), 4641; https://doi.org/10.3390/ijms24054641
Submission received: 12 December 2022 / Revised: 23 February 2023 / Accepted: 24 February 2023 / Published: 28 February 2023
(This article belongs to the Section Molecular Plant Sciences)
Browse Figures
Versions Notes

Abstract

:
Flower color is one of the crucial traits of ornamental plants. Rhododendron delavayi Franch. is a famous ornamental plant species distributed in the mountain areas of Southwest China. This plant has red inflorescence and young branchlets. However, the molecular basis of the color formation of R. delavayi is unclear. In this study, 184 MYB genes were identified based on the released genome of R. delavayi. These genes included 78 1R-MYB, 101 R2R3-MYB, 4 3R-MYB, and 1 4R-MYB. The MYBs were divided into 35 subgroups using phylogenetic analysis of the MYBs of Arabidopsis thaliana. The members of the same subgroup in R. delavayi had similar conserved domains and motifs, gene structures, and promoter cis-acting elements, which indicate their relatively conserved function. In addition, transcriptome based on unique molecular identifier strategy and color difference of the spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex were detected. Results showed significant differences in the expression levels of R2R3-MYB genes. Weighted co-expression network analysis between transcriptome and chromatic aberration values of five types of red samples showed that the MYBs were the most important TFs involved in the color formation, of which seven were R2R3-MYB, and three were 1R-MYB. Two R2R3-MYB (DUH019226.1 and DUH019400.1) had the highest connectivity in the whole regulation network, and they were identified as hub genes for red color formation. These two MYB hub genes provide references for the study of transcriptional regulation of the red color formation of R. delavayi.

1. Introduction

Anthocyanins are flavonoids and the main pigments in petals and fruits. These flavonoids can give plants different colors to attract insects and help spread pollen and seeds to promote plant reproduction [1]. Anthocyanins also have anti-oxidation and anti-inflammatory activities [2].
Anthocyanin biosynthesis is controlled by numerous genes, such as transcription factor (TF) and structural genes. Anthocyanidin 3-O-glucosyltransferase (UDPGT, in short UGT) is a structural gene and plays a role in the anthocyanin biosynthesis of plants, such as Capsicum annuum [3], Morella rubra [4], and Vaccinium spp. [5]. The expression of structural genes is often regulated by TFs. MYB is one of the largest families among plant TFs. This TF can bind to cis elements in the promoters of target genes to mainly regulate plant anthocyanin biosynthesis, development, and stress resistance [6]. The N-terminus of MYB has a highly conserved DNA-binding domain (DBD). The domain contains three irregular repeats, each of which forms three A-helices. The other two helices construct a Helix-turn-Helix (HTH) structure with three evenly spaced tryptophan (Trp) residues. Hydrophobic nuclei are formed in the three-dimensional HTH structure [7]. Based on the number of repeats of domains, MYBs can be classified into 1R-MYB, R2R3-MYB, R3-MYB, and 4R-MYB [8]. Among the four types of MYB, R2R3-MYB plays the most central role in anthocyanin synthesis [9].
TFs activate/inhibit the genes that encode various catalytic enzymes to regulate anthocyanin biosynthesis [10]. MYB alone, or by forming a MBW complex with basic helix-loop-helix (bHLH) and WD40, regulates anthocyanin biosynthesis by binding to specific regions of the target gene promoter [11]. The function of R2R3-MYB in regulating anthocyanin biosynthesis has been reported in numerous species. In apple fruit, MdMYB306 interacts with MdMYB17 and MdbHLH33 through its N-terminal to inhibit anthocyanin synthesis and affect fruit color [12]. In Pistacia chinensis, PcMYB113 makes the leaves red by promoting anthocyanin biosynthesis during autumn leaf coloration [13]. MYB also regulates the anthocyanin biosynthesis in peony [14], pear [15], pomegranate [16], and kiwi fruit [17].
High-throughput RNA sequencing (RNA-Seq) can generate information-rich sequence and expression data sets to characterize the abundance of related genes in plants [18,19]. RNA-Seq using a unique molecular identifier (UMI) strategy can eliminate quantitative interference by polymerase chain reaction (PCR) amplification preferences, leading to more accurate expression results [20,21]. Hub genes associated with a certain trait can be identified using the weighted gene co-expression network analysis (WGCNA) of RNA-Seq genes based on the UMI strategy associated with trait phenotypic data [22,23,24].
Rhododendron delavayi is a small evergreen tree belonging to the subgenus Hymenanthes (Rhododendron, Ericaceae). This species is widely distributed throughout Southwest China with large inflorescences, and it is an important dominant species in the nature reserve. As one of the world’s famous ornamental plants, the inflorescence of R. delavayi consists of more than 10 small flowers. The flowers and branchlets of this species are all red and contain high concentrations of anthocyanins [25]. Previous studies have shown that overexpression of RdDFR1, Rd3GT1, and Rd3GT6 changed the transgenic tobacco flower color from light pink to dark pink [26,27]. However, no TF genes have been reported except for the structural genes related to the anthocyanin biosynthesis of R. delavayi.
Therefore, in view of the vital regulatory role of MYBs in the anthocyanin biosynthesis of other plants, in this study, the MYB TF family was analyzed based on the released R. delavayi genome [28]. In addition, the expression profiles of the members of this family were detected by RNA-Seq using the UMI strategy on five types of colorful samples, namely, spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex. With the WGCNA of the transcriptome and color difference values, the hub genes related to the red color formation of R. delavayi were identified. The data will provide an important reference for the study of the transcriptional regulation mechanism of the color formation of R. delavayi.

2. Results

2.1. Genome-Wide Identification of the Members of R. delavayi MYB Gene Family

The whole released genome of R. delavayi [28] was analyzed to identify the MYB gene family. A total of 202 MYB genes were selected by BLAST against the genome by using the protein sequences of 168 MYB genes in Arabidopsis thaliana. In addition, 292 MYB genes were hit by using the Hidden Markov Model (HMM) to query the MYB DBD (PF00249). Combining the blast results of the two different methods, a total of 184 MYB members were finally identified, including 78 1R-MYB, 101 R2R3-MYB (2R-MYB), 4 R1R2R3-MYB (3R-MYB), and one 4R-MYB protein (Figure 1) as analyzed against Pfam, InterPro, and Prosite databases. The isoelectric point of the identified MYBs ranged from 4.45 (DUH002843.1) to 10.42 (DUH004390.1), and the predicted molecular weights ranged from 7412.46 kDa (DUH004390.1) to 179,531.3 kDa (DUH000629.1) (Table S1). In addition, subcellular localization prediction of MYBs showed that 175 of the 184 MYB proteins were located in the nucleus, five in the cytoplasm, and four in mitochondria (Table S1).

2.2. Analysis of the MYB Conserved Motif and Gene Structure of R. delavayi

To study the sequence characteristics of the conserved MYB DBD of R. delavayi, we performed multiple sequence alignments on their amino acid sequences. The MYB repeats on R2 and R3 of R. delavayi contained characteristic amino acids (Figure 2). Three highly conserved Trp residues were found at positions 5, 25, and 45 in the R2 repeats of R2R3-MYB, whereas two Trp residues were found at positions 28 and 47 in R3 (Figure 2). The phylogenetic tree was constructed by the ML method using the 184 MYB protein sequences of R. delavayi. The MYB family was clustered into different branches (Figure 3). Ten conserved motifs of MYB proteins were identified by the MEME program. Among all MYBs, R2R3-MYB, 3R-MYB, and 4R-MYB genes had high similarity in the amino acid sequences (Figure 3A). The 1R-MYB protein contained motifs 3 and 4. Other types of MYB proteins basically contained motifs 1, 2, 3, and 5 (Figure 3A). The MYB proteins of R. delavayi in the same group had similar motif types and quantities, which indicates that these proteins may have similar functions.
To further measure the structural characteristics of R. delavayi MYB genes, we further analyzed their introns/exons. The most homologous genes in a group shared the same gene structure layout and the number of exons and introns (Figure 3B). All 3R-MYB contained seven exons, whereas 4R-MYB contained 10 exons. R2R3-MYB contained one–four exons, 68% of which contained conserved gene structure with three exons and two introns. In addition, 1R-MYB contained 1–13 exons, 51% of which had more than four exons. Some introns of the genes were large and accounted for a big part of the genes. Compared with CDS sequences, most MYB genes are with or without two untranslated regions (UTRs). The UTR sequences are less similar, and their diversity is expressed in terms of gene length, which indicates that UTR sequences are less conserved than the coding regions.

2.3. Phylogenetic Analysis of the MYB Family of R. delavayi

To study the evolutionary relationship of MYB genes in R. delavayi, we downloaded 132 protein sequences of A. thaliana MYBs from TAIR, together with 184 protein sequences of R. delavayi MYBs, and constructed a ML phylogenetic tree. The MYB members of R. delavayi were divided into 35 subgroups and named C1 to C35 based on sequence similarity and topological structure, according to the classification of the same family in A. thaliana. The clades (S1–S25 subgroups) were labeled based on the evolutionary relationship in A. thaliana in the tree (Figure 4). The phylogenetic tree showed that most of the R. delavayi MYBs were clustered with R2R3-MYB of A. thaliana. Subgroup C31 contained a maximum of 19 members, whereas subgroup C20 had none. Four subgroups, namely, C21, C31, C34, and C35, were not grouped with the MYB of A. thaliana. The results suggest that they may have undergone gene gain or loss events during evolution.

2.4. Analysis of R. delavayi MYB Genes Promoters Cis-Acting Elements

Cis-acting elements in the promoter region of R. delavayi MYB genes can be divided into three categories: plant hormone, biological/abiotic stress, and plant growth and development response elements (Figure 5 and Table S2). Plant hormone response elements, such as Me-JA response elements (e.g., CGTCA and TGACG motifs), abscisic acid response elements (ABRE), gibberellin response elements (e.g., P-box), and salicylic acid response elements (e.g., TCA element), were widely distributed in MYB promoters. The Me-JA elements accounted for the highest proportion (40%). Cis elements involved in abiotic stress responses included anaerobic inducible elements (ARE), low-temperature response elements (LTR), and drought-inducible elements (e.g., MBS). Among the growth and development response elements, the light (e.g., G-Box) and auxin response elements (e.g., TGA element) accounted for 77% and 11%, respectively. The defense and stress response elements (e.g., TC-rich repeats) accounted for 9%, but the trauma response elements (e.g., WUN motif) accounted for 1%. Flavonoid biosynthesis response elements (e.g., MBSI) accounted for 2%. All cis elements appeared most frequently in the promoters of 2R-MYB, and different subgroups of MYB promoters had distinct cis element types. One TCA element and one WUN motif were found in 1R-MYB and 2R-MYB promoters, respectively, which indicates the functional diversity of MYB in R. delavayi. Moreover, a total of 101,112 TF-binding sites and 223 potentially interacting TFs described previously in A. thaliana were revealed using the PlantPan 3.0 (Table S3). MYB binding sites were identified at all MYB promoters examined. There were also several binding sites of TFs other than the MYB family on the examined promoters, such as the binding sites of bZIP, bHLH, and WRKY.

2.5. Transcriptomic Analysis

To determine the MYB genes that play a role in the color formation of R. delavayi, we performed RNA-seq analysis based on the UMI strategy using spotted petals, unspotted petals, spotted throat, unspotted throat of flower, and branchlet cortex (Figure 6A). The produced raw data of the transcriptome have been uploaded to the Sequence Reading Archive database (http://www.ncbi.nlm.nih.gov/sra/, accessed on 3 December 2022) with entry number PRJNA907866. After the removal of low-quality reads and adapters, an average of 23,879,787 clean reads were obtained per sample. The percentage of bases with a phred value greater than 20 in the total bases (QC20), the percentage of bases with a phred value greater than 30 in the total bases (QC30), and GC percentage were greater than 97%, 92%, and 47%, respectively (Table S4). In these clean reads, 41,484,380 (90.77%) were mapped to the R. delavayi reference genome (Table S5). Fragments per kilobase of exon models per million mapped fragments (FPKM) were used to determine the expression level of unigenes. Based on the FPKM values, principal component analyses (PCA) and Pearson correlation coefficients showed high correlation coefficients between three biological replicates (Figure 6B,C, respectively).

2.6. Differentially Expressed Genes (DEGs), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis

The DEGs and shared genes between samples were analyzed. Comparative analysis showed 6736 DEGs between MY-1 and MY-2 (3090 up-regulated genes and 3646 down-regulated genes) (Figure 7A); 5722 DEGs between MY-3 and MY-4 (2841 upregulated and 2881 downregulated) (Figure 7A); 15,847 DEGs between MY-5 and MY-2 (8160 upregulated and 7687 downregulated) (Figure 7A); 15,806 DEGs between MY-5 and MY-4 (8104 upregulated and 7702 downregulated) (Figure 7A). In addition, the Venn diagram showed the following: 1160 genes shared expressions in MY-1 and MY-2; 940 genes shared expressions in MY-3 and MY-4; 1586 genes shared expressions in MY-5 and MY-2; 1538 genes shared expressions in MY-5 and MY-4 (Figure 7B). GO and KEGG enrichment analyses of different sample combinations revealed that these DEGs are related to various metabolic and biosynthetic pathways (Figure S1 and Figure Figure 7C), and the DEGs of all four combinations were enriched in the pathways encoding phenylpropanoid biosynthesis (ko00940), flavonoid biosynthesis (ko00941), and anthocyanin biosynthesis (ko00942). Thus, DEGs may be important for anthocyanin biosynthesis in R. delavayi.

2.7. Expression Analyses of R2R3-MYB Genes in R. delavayi

The MYB TFs play an important role in anthocyanin biosynthesis, especially R2R3-MYB. To evaluate the expression pattern of R2R3-MYB genes in anthocyanin biosynthesis in different tissues of R. delavayi, we evaluated the expression levels of R2R3-MYB genes based on the transcriptome data of spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex. A total of 91 R2R3-MYB genes were identified, and 85 of them showed differential expression in five types of red samples (Figure 8). Among these genes, 10 were significantly expressed in four tissues, including spotted petals, non-spotted petals, spotted larynx, and non-spotted larynx, whereas 34 were significantly expressed only in the branchlet cortex. The expression of certain genes was limited to certain tissues, such as MYB genes DUH010372.1 and DUH007487.1 which were only significantly expressed in unspotted petals and spotted throat, respectively. Heat maps of R2R3-MYB expression in R. delavayi showed that most MYB genes were expressed differently in the five types of samples, and several genes were expressed only in the branchlet cortex.

2.8. Analysis of Co-Expression Network

Tissue colors were determined using the CIE 1976 (L*a*b*) system [29], that is, the larger the value of L*, the brighter the color, and the smaller the value, the darker the color; the larger the value of a*, the more reddish, and the smaller the value, the more greenish the tissue; the larger the value of b*, the more yellowish, and the smaller the value, the more bluish the tissue [30,31]. To screen hub genes regulating the red-color formation of R. delavayi, the correlation of color indicators of L*, a*, and b* values (Table 1) and the gene expression levels based on the RNA-Seq data of the five types of red samples were analyzed using WGCNA. A total of 19 co-expression modules were detected in the screened unigenes (Figure 9). Modules ranged in size from 45 unigenes (gray module) to 5225 unigenes (turquoise module). For the a* value, the a* indicator represents the transition from red to green [30,31], which was positively correlated with the MEturquoise module, with the highest correlation and strongest significance (r = 0.9, p = 7 × 10−6). A total of 23 highly connective genes (degree ≥ 4962), including 10 MYBs, six bHLHs, one WD40 gene, two MYB-relate genes, and four structural genes (UDPGT) were screened in the module (Figure 10). Among the 10 MYBs, seven were R2R3-MYB genes (DUH019226.1, DUH019400.1, DUH013179.1, DUH007661.1, DUH003273.1, DUH022833.1, and DUH022115.1), and three were 1R-MYB genes (DUH007734.1, DUH012270.1, and DUH021327.1). These results indicate that the MYBs, bHLHs, WD40, and UDPGT identified are probably involved in the red color formation of R. delavayi. In addition, the co-expression network showed that two R2R3-MYBs (DUH019226.1 and DUH019400.1) had the highest connectivity in the whole regulation network and were identified as hub genes for red color formation. The identified bHLH, WD40, and UDPGT genes may also be involved in red color formation in R. delavayi.

2.9. Quantitative Reverse Transcription qRT-PCR Analysis of Gene Expression

Based on the co-expression network analysis, nine genes were randomly selected to further verify the accuracy of transcriptome results. Fold changes in the different tissues relative to the branchlet cortex were analyzed using qRT-PCR. The results showed that all genes detected by qRT-PCR had expression patterns similar to those of transcriptome data, and the relative expression was significantly different in petal and vegetative tissues (Figure 11). UDPGT (DUH009790.1) and MYB (DUH013179.1) were differentially expressed in all five types of red samples. Meanwhile, MYB (DUH019226.1, DUH012270.1, and DUH003273.1), bHLH (DUH011332.1), and WD40 (DUH022768.1) had high relative expression levels in the branchlet cortex. However, MYB-related genes (DUH031703.1) and MYB (DUH019400.1) had high expression in flower tissues, including spotted petals, unspotted petals, spotted throat, and unspotted throat. These qRT-PCR results further confirmed the gene repression reliability of the transcriptome.

3. Discussion

MYB TFs are one of the largest gene families in plants and play a crucial role in anthocyanin biosynthesis. This family has been systematically studied in A. thaliana, rice, tomato, petunia, and kiwi fruit [8,32,33,34,35,36]. A total of 184 members of the MYB family were identified in R. delavayi. This family has more members in R. delavayi than in rice (155) and A. thaliana (168), whereas the number was lower than those in peach (256), apple (229), and other plant species [37,38]. Overall, the number is greater than 100 for different species, which indicates that this family is greatly amplified in higher plants [39]. In R. delavayi, the R2R3-MYB subfamily is the most abundant in the MYB family, consistent with previous studies that reported 110 members of R2R3-MYB in rice [40] and 126 members in A. thaliana [9]. In this study, the majority of R2R3-MYB genes (68%) in R. delavayi contained a conserved gene structure with three exons and two introns, which are also found in other plant species [41]. Most homologous genes in the same group of the family have a common exon/intron structure and exon number, which indicates a common evolutionary history. Besides, some introns of the genes were large and accounted for a big part of the genes, the length of introns sometimes are also related to some process of adaptation in evolution [42].
Phylogenetic analysis revealed clustered homologous genes in the same clades and subclades, which often exhibit similar functions [43]. The phylogenetic trees constructed from 132 A. thaliana MYBs and 184 R. delavayi MYBs further illustrated the evolutionary relationship of R. delavayi MYBs. Based on the classification of MYB in A. thaliana [9], the MYB members of R. delavayi were divided into 35 subgroups (C1–C35) (Figure 4), among which four subgroups (C21, C31, C34, and C35) of R. delavayi did not cluster with the MYBs of A. thaliana. This result indicated that the six subgroups may have acquired or lost genes during evolution compared with A. thaliana. In addition, most of the clades had different numbers of A. thaliana MYBs and R. delavayi MYBs, and genes from the same clades may share a common evolutionary process. Based on the reported classification and function of MYBs in A. thaliana, phylogenetic trees would be helpful to predict the function of R. delavayi MYB. R2R3-MYBs of A. thaliana in the subclades of S3, S4, S5, S6, and S7 regulated phenylpropanoid biosynthesis [8]. This result indicated that R. delavayi MYBs clustered in the same subclades may also be involved in regulating phenylpropanoid biosynthesis.
Transcriptome analysis had been widely used to study gene regulation in horticultural traits [43,44,45]. In this study, KEGG analysis showed that the petals and flower throat had enriched genes in the pathway of phenylpropanoid and flavonoid biosynthesis. In addition, R2R3-MYB TFs in the reported plant species, such as A. thaliana [8], strawberry [46], and kiwi fruit [47], are involved in the activation of structural genes related to phenylpropanoid biosynthesis pathway and regulation of biosynthesis and accumulation of anthocyanins. Similarly, in this study, 86 R2R3-MYB TFs of R. delavayi were differentially expressed in the differently colored tissues. Ten R2R3-MYB genes were significantly expressed in the tissues from spotted petals, unspotted petals, spotted throat, and unspotted throat. Several R2R3-MYB genes were expressed only in specific tissues, which suggests that these R2R3-MYB genes play a key role in various aspects of color formation or the development of R. delavayi [34]. In addition, an analysis of cis-acting elements in R. delavayi showed that almost all MYB genes contained cis-acting elements, which respond to hormone signal transduction, abiotic stress, and plant growth and development. This result is consistent with previous studies on Petunia [34] and Actinidia, in which MYB promoters also contained these cis-acting elements [47].
TFs play an important role in color formation by regulating the temporal and spatial expression of structural genes [48]. MYB can form MBW complexes with bHLH and WD40 to activate or inhibit the transcription of target genes to regulate anthocyanin synthesis [48,49]. MdMYB10 affects the anthocyanin biosynthesis of Malus spectabilis by regulating the expression of ANS [50]. PsMYB10.2 promotes anthocyanin accumulation in Prunus salicina fruits by activating PsUFGT and PsGST [51]. GbBM of Gossypium barbadense encoding MYB113 directly targets the promoter of four flavonoid biosynthesis genes and positively regulates the development of petal spots [52]. The interaction between AcMYB10 and AcbHLH5 in A. chinensis promotes the expression of AaF3H and AaF3G, respectively, and this process is conducive to anthocyanin accumulation in the flesh [53]. These findings suggest that MYBs are also widely involved in the regulation of anthocyanin biosynthesis in non-model plant species. In addition, WGCNA showed that DcMYB113 of Daucus carota activates the expressions of DcbHLH3 and anthocyanin biosynthesis-related structural genes [54]. Eleven genes (e.g., PaMYB10 and seven structural genes) were identified in the Prunus armeniaca co-expression network. PaMYB10, MdMYB10, and PsMYB10 belong to the anthocyanin-associated R2R3-MYB branch, showing high homology [55]. Similar to previous reports, GbBM and DcMYB113 belong to the S6 clade of the R2R3 MYB. In this study, WGCNA showed that 10 MYBs, six bHLHs, one WD40, two MYB-related genes, and four UDPGTs genes constructed the co-expression network behind the coloration in R. delavayi. Among these genes, two R2R3-MYB TF genes (DUH019400.1 and DUH019226.1) were identified as hub genes associated with anthocyanin biosynthesis in R. delavayi. qRT-PCR analysis revealed that the MYB TF DUH019226.1 was significantly expressed in the branchlet cortex, and DUH019400.1 was significantly expressed in the tissues of the remaining four petals. DUH019226.1 in the phylogenetic tree clustered into the S4 subgroup related to anthocyanin synthesis [5]. DUH019400.1 was clustered into the S19 subgroup, which controls another development or function [48]. The results indicate the possible gene acquisition or loss events during the evolution of R. delavayi MYBs.

4. Materials and Methods

4.1. Plant Materials and Color Determination

Plant samples of R. delavayi were collected from Baili Rhododendron Reserve, Guizhou Province, China. The samples were obtained from five types of red samples, namely, spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex. The samples were randomly collected with three biological replicates. After sampling, the samples were immediately frozen in liquid nitrogen and stored at −80 °C. Color deference was represented using the CIE 1976 (L*a*b*) system, the color values of all samples were determined by a portable color spectrometer (CHN Spe, Hangzhou, China) before sampling, and the values of L*, a*, and b* were determined with three replicates.

4.2. Identification and Sequence Analysis of the MYB Gene Family

Whole released genome sequence of R. delavayi was downloaded from the Rhododendron Plant Genome Database (RPGD) (http://bioinfor.kib.ac.cn/RPGD/index.html, accessed on 1 October 2022) [28]. The protein sequences of 168 MYBs of A. thaliana were downloaded from the Plant Transcription Factor Database (http://planttfdb.gao-lab.org/, accessed on 3 October 2022) [8]. All MYBs of R. delavayi were identified using blastp (E < 1 × 10−5) with MYBs of A. thaliana queries, which is available from the NCBI (http://www.ncbi.nlm.nih.gov/, accessed on 3 October 2022). The DBD (PF00249) of MYB was downloaded from the Pfam database (http://pfam.sanger.ac.uk/, accessed on 4 October 2022) and used as a query. HMMER (http:/hmmer.janelia.org/, accessed on 4 October 2022) version 3.1b2 software was used to search and identify the MYB candidate sequence in the database of R. delavayi, and the cut-off of E value was set to <1.0 × 10−5 [56]. Then, the intersection of results from the two different methods elucidated common elements. To verify the reliability of the results, candidate MYB proteins were further analyzed using Pfam, InterPro (http://www.ebi.ac.uk/interpro/, accessed on 5 October 2022), and PROSITE [57] (https://prosite.expasy.org/, accessed on 5 October 2022). Sequences lacking the central or entire MYB DBDs were filtered out. The selected MYB proteins of R. delavayi molecules (molecular weight and isoelectric point) were calculated through Expasy (https://web.expasy.org/compute_pi/, accessed on 7 October 2022). Subcellular localization of MYBs was predicted using Cello (http://cello.life.nctu.edu.tw, accessed on 7 October 2022) [58].

4.3. Analysis of the MYB Conserved Motif and Gene Structure

The gene structure of MYBs was analyzed using TBtools software v1.0987663 (https://github.com/CJ-Chen/TBtools/releases, accessed on 12 October 2022). MYB exon/intron was visualized using the TBtools [59]. The conserved motif of MYB was determined by MEME (https://meme-suite.org/meme/, accessed on 12 October 2022) [60], with an optimal motif length of 6–200 residues and a search number of 10. The conservative DBD characteristics of MYB were analyzed using WebLogo (http://weblogo.berkeley.edu/logo.cgi, accessed on 1 October 2022) [61]. All software was used under default parameter settings.

4.4. MYB Sequence Alignment and Phylogenetic Analysis

Muscle was used for multiple sequence alignment of MYBs. MEGAX v7.0.14 was used to explore the evolutionary relationship between the MYB sequences of A. thaliana as a reference sequence [62]. A phylogenetic tree (ML) was constructed using Muscle to compare the full lengths of MYB amino acid sequences (184 R. delavayi MYBs and 132 A. thaliana MYBs). Parameters were set to WAG + G, paired deletion, and bootstrap analysis with 1000 bootstrap replicates. R. delavayi MYBs were classified based on the phylogenetic relationship with the corresponding classification of MYBs in A. thaliana.

4.5. Analysis of Promoter Cis-Acting Elements of MYB

The upstream 2000 bp sequence of MYB CDSs of R. delavayi was extracted from the RPGD as promoters [28]. Putative binding sites for TF candidates were determined using the PlantPAN3.0 server (http://plantpan.itps.ncku.edu.tw/, accessed on 9 February 2023) and A. thaliana as references [63].Cis elements of the promoters were predicted using PlantCARE [64] to identify the three main types of cis-acting elements on the promoter region: phytohormone response, including Me-JA response element (CGTCA and TGACG motifs), ABRE, gibberellin response element (P-box), and salicylic acid response element (TCA element); biotic/abiotic stress response, including ARE, LTR, and drought induction element (MBS); plant growth and development, including light response element (G-Box), growth hormone response element (TGA element), defense and stress response element (TC-rich repeats), trauma response element (WUN-motif), and flavonoid biosynthesis (MBSI). Statistical analysis and cis-element visualization were performed using Excel and TBtools, respectively [59].

4.6. Transcriptome Analysis of Five Types of Red Samples

Total RNA was extracted from samples using the RNAprep Pure Polysaccharide Polyphenol Plant Kit (TianGen, Beijing, China). The qualified RNA samples were sent to Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). The libraries were sequenced on Illumina Novaseq 6000 Platform and generated 150 nt pair-end reads. Raw reads were first processed through in-house perl scripts. Clean reads were obtained by removing reads containing adapter and/or ploy-N, and low-quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. UMI sequences on each read were identified with UMI-tools (1.0.0) [65]. Genome of R. delavayi (http://bioinfor.kib.ac.cn/RPGD/index.html, accessed on 1 October 2022) [28,66] was used as reference to assemble the transcriptome using HISAT2 v2.0.4 [67]. The number of reads mapped to each gene was calculated using HTSeq v0.6.1 [68]. The FPKM of each gene was calculated based on the gene length and counts mapped to the corresponding gene.

4.7. Analysis of DEGs, GO, and KEGG

DESeq R package (1.10.1) was used to analyze the differential expression between two conditions/groups (two biological replicates per condition). Corrected p-value of 0.005 and log2 (fold change) of 1 were set as the threshold for significance. GO analysis of DEGs was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with corrected p-value less than 0.05 were considered significantly enriched by DEGs. In KEGG enrichment analysis, KOBAS 3.0 (http://kobas.cbi.pku.edu.cn/index.php, accessed on 20 May 2022 ) database was used to test the statistical enrichment of DEGs.

4.8. Expression Analyses of R2R3-MYB Genes in R. delavayi

The transcript abundances of R. delavayi R2R3-MYB genes were calculated with FPKM values. Heatmaps of R2R3-MYB gene expression data were generated using TBtools v1.0987663 [59]. Gene expression levels were expressed by log2 FPKM values.

4.9. Analysis of Weight Co-Expression Network

WGCNA was performed using the R v1.63 software package. Genes with coefficient of variation less than 0.5 were removed, and cluster trees were constructed based on the correlation and gene expression levels. In the WGCNA network, the soft threshold was set to 5 (Figure S2), the module size was set to at least 30, and the merge cut height value was set to 0.25 (Figure S3). The correlation between the modules and the chromatic aberration values of the five types of red samples was analyzed for all genes in each module. The significant trait correlation module was identified based on the high correlation values and significance (p-value). Finally, Cytoscape v3.9.1 [69] was used to display the weight co-expression network of the selected modules, calculate the degree values of the whole network, and screen the hub genes with high connectivity.

4.10. qRT-PCR Verification

Expression patterns of nine selected genes were verified by qRT-PCR. The qualified RNA (1 µg) was reverse transcribed using TRUEscript RT Kit with gDNA clearance (Aidlab, Beijing, China). Specific primers were designed by Primer v6.0 (ref), and 18S was used as an internal control [70] (Table S6). Three biological replicates and three technological replicates of qRT-PCR were performed for verification. Each reaction contained 10 µL SYBR qPCR Mix (Aidlab, China), 1 µL cDNA, 10.5 µL distilled water, 0.5 µL forward primers, and 0.5 µL reverse primers under the following conditions: 94 °C for 2 min, 40 cycles; 94 °C for 30 s; 60 °C for 20 s; 72 °C for 30 s; 4 °C. The relative expression level of genes was calculated using 2−ΔΔct.

5. Conclusions

In this study, 184 MYB genes were identified in the R. delavayi genome and classified into 35 subgroups based on phylogenetic analyses. Bioinformatics analysis was performed including, conserved domains and motifs, gene structures, promoter cis-acting elements, and transcriptomic profiling. In addition, the expression levels of the R2R3-MYB gene were significantly different in five types of red samples. WGCNA between transcriptome and chromatic aberration values showed that 10 MYBs, six bHLHs, one WD40, two MYB-related genes, and four structural genes UDPGTs constructed the co-expression network behind the coloration in R. delavayi. The MYBs were the most important transcription factors involved in the red color formation, among them, two R2R3-MYB (DUH019226.1 and DUH019400.1) were identified as hub genes. Furthermore, the qRT-PCR validated the results of the transcriptome. This work has enhanced a comprehensive understating of the MYB gene family in R. delavayi and provided references for the study of transcriptional regulation of the red color formation of R. delavayi.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms24054641/s1.

Author Contributions

F.L. and H.L. designed the experiments scheme. F.L., H.W., W.Z. and Q.A. performed sampling and experiment, F.L. wrote the original manuscript, H.L refined the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (32260415).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

R. delavayi genome data can be searched in open access RPGD database (http://bioinfor.kib.ac.cn/RPGD/index.html, accessed on 1 October 2022). The raw data of transcriptome has been uploaded to the Sequence Reading Archive database (http://www.ncbi.nlm.nih.gov/sra/, accessed on 3 December 2022) with the entry number PRJNA907866. Additional data presented in this study are presented in the supplementary materials.

Acknowledgments

We appreciated the research team that released the genome of R. delavayi.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Winkel-Shirley, B. Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology. Plant Physiol. 2001, 126, 485–493. [Google Scholar] [CrossRef] [PubMed]
  2. Liu, J.; Zhou, H.; Song, L.; Yang, Z.; Qiu, M.; Wang, J.; Shi, S. Anthocyanins: Promising natural products with diverse pharmacological activities. Molecules 2021, 26, 3807. [Google Scholar] [CrossRef] [PubMed]
  3. Zhou, Y.; Mumtaz, M.A.; Zhang, Y.; Shu, H.; Hao, Y.; Lu, X.; Cheng, S.; Zhu, G.; Wang, Z. Response of anthocyanin accumulation in pepper (Capsicum annuum) fruit to light days. Int. J. Mol. Sci. 2022, 23, 8357. [Google Scholar] [CrossRef]
  4. Ren, C.; Cao, Y.; Xing, M.; Guo, Y.; Li, J.; Xue, L.; Sun, C.; Xu, C.; Chen, K.; Li, X. Genome-wide analysis of UDP-glycosyltransferase gene family and identification of members involved in flavonoid glucosylation in Chinese bayberry (Morella rubra). Front. Plant Sci. 2022, 13, 998985. [Google Scholar] [CrossRef] [PubMed]
  5. Lin, Y.; Wang, Y.; Li, B.; Tan, H.; Li, D.; Li, L.; Liu, X.; Han, J.; Meng, X. Comparative transcriptome analysis of genes involved in anthocyanin synthesis in blueberry. Plant Physiol. Biochem. 2018, 127, 561–572. [Google Scholar] [CrossRef] [PubMed]
  6. Wang, X.; Niu, Y.; Zheng, Y. Multiple functions of MYB transcription factors in abiotic stress responses. Int. J. Mol. Sci. 2021, 22, 6125. [Google Scholar] [CrossRef]
  7. Ogata, K.; Kanei-Ishii, C.; Sasaki, M.; Hatanaka, H.; Nagadoi, A.; Enari, M.; Nakamura, H.; Nishimura, Y.; Ishii, S.; Sarai, A. The cavity in the hydrophobic core of MYB DNA-binding domain is reserved for DNA recognition and trans-activation. Nat. Struct. Biol. 1996, 3, 178–187. [Google Scholar] [CrossRef]
  8. Dubos, C.; Stracke, R.; Grotewold, E.; Weisshaar, B.; Martin, C.; Lepiniec, L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010, 15, 573–581. [Google Scholar] [CrossRef]
  9. Stracke, R.; Werber, M.; Weisshaar, B. The R2R3-MYB gene family in Arabidopsis thaliana. Curr. Opin. Plant Biol. 2001, 4, 447–456. [Google Scholar] [CrossRef]
  10. Wiltshire, E.J.; Eady, C.C.; Collings, D.A. Induction of anthocyanin in the inner epidermis of red onion leaves by environmental stimuli and transient expression of transcription factors. Plant Cell Rep. 2017, 36, 987–1000. [Google Scholar] [CrossRef]
  11. Sun, X.; Zhang, Z.; Li, J.; Zhang, H.; Peng, Y.; Li, Z. Uncovering hierarchical regulation among MYB-bHLH-WD40 proteins and manipulating anthocyanin pigmentation in rice. Int. J. Mol. Sci. 2022, 23, 8203. [Google Scholar] [CrossRef] [PubMed]
  12. Wang, S.; Zhang, Z.; Li, L.-X.; Wang, H.-B.; Zhou, H.; Chen, X.-S.; Feng, S.-Q. Apple MdMYB306-like inhibits anthocyanin synthesis by directly interacting with MdMYB17 and MdbHLH33. Plant J. 2022, 110, 1021–1034. [Google Scholar] [CrossRef] [PubMed]
  13. Song, X.; Yang, Q.; Liu, Y.; Li, J.; Chang, X.; Xian, L.; Zhang, J. Genome-wide identification of Pistacia R2R3-MYB gene family and function characterization of PcMYB113 during autumn leaf coloration in Pistacia chinensis. Int. J. Biol. Macromol. 2021, 192, 16–27. [Google Scholar] [CrossRef]
  14. Luan, Y.; Tang, Y.; Wang, X.; Xu, C.; Tao, J.; Zhao, D. Tree peony R2R3-MYB transcription factor PsMYB30 promotes petal blotch formation by activating the transcription of the anthocyanin synthase gene. Plant Cell Physiol. 2022, 63, 1101–1116. [Google Scholar] [CrossRef]
  15. Feng, S.; Xu, Y.; Yang, L.; Sun, S.; Wang, D.; Chen, X. Genome-wide identification and characterization of R2R3-MYB transcription factors in pear. Sci. Hortic. 2015, 197, 176–182. [Google Scholar] [CrossRef]
  16. Khaksar, G.; Tabatabaei, B.E.S.; Arzani, A.; Ghobadi, C.; Ebrahimie, E. Functional analysis of a pomegranate (Punica granatum L.) MYB Transcription factor involved in the regulation of anthocyanin biosynthesis. Iran. J. Biotechnol. 2015, 13, 17–25. [Google Scholar] [CrossRef]
  17. Wang, W.-q.; Moss, S.M.A.; Zeng, L.; Espley, R.V.; Wang, T.; Kui, L.-W.; Fu, B.-l.; Schwinn, K.E.; Allan, A.C.; Yin, X.-r. The red flesh of kiwifruit is differentially controlled by specific activation-repression systems. New Phytol. 2022, 235, 630–645. [Google Scholar] [CrossRef]
  18. Gongora-Castillo, E.; Buell, C.R. Bioinformatics challenges in de novo transcriptome assembly using short read sequences in the absence of a reference genome sequence. Nat. Prod. Rep. 2013, 30, 490–500. [Google Scholar] [CrossRef]
  19. Zhang, H.; He, L.; Cai, L. Transcriptome sequencing: RNA-Seq. In Methods in Molecular Biology; Springer: Berlin/Heidelberg, Germany, 2018; Volume 1754, pp. 15–27. [Google Scholar] [CrossRef]
  20. Belouah, I.; Benard, C.; Denton, A.; Blein-Nicolas, M.; Balliau, T.; Teyssier, E.; Gallusci, P.; Bouchez, O.; Usadel, B.; Zivy, M.; et al. Transcriptomic and proteomic data in developing tomato fruit. Data Brief 2020, 28, 105015. [Google Scholar] [CrossRef]
  21. Li, H.; Hu, Y.; Gao, C.; Guo, Q.; Deng, Q.; Nan, H.; Yang, L.; Wei, H.; Qiu, J.; Yang, L. Integrated SMRT technology with UMI RNA-Seq reveals the hub genes in stamen petalody in Camellia oleifera. Forests 2021, 12, 749. [Google Scholar] [CrossRef]
  22. Garcia-Ruiz, S.; Gil-Martinez, A.L.; Cisterna, A.; Jurado-Ruiz, F.; Reynolds, R.H.; Cookson, M.R.; Hardy, J.; Ryten, M.; Botia, J.A. CoExp: A web tool for the exploitation of co-expression networks. Front. Genet. 2021, 12, 630187. [Google Scholar] [CrossRef]
  23. Sanchez-Baizan, N.; Ribas, L.; Piferrer, F. Improved biomarker discovery through a plot twist in transcriptomic data analysis. BMC Biol. 2022, 20, 1. [Google Scholar] [CrossRef]
  24. DiLeo, M.V.; Strahan, G.D.; den Bakker, M.; Hoekenga, O.A. Weighted Correlation Network Analysis (WGCNA) Applied to the Tomato Fruit Metabolome. PLoS ONE 2011, 6, 10. [Google Scholar] [CrossRef] [PubMed]
  25. Du, H.; Lai, L.; Wang, F.; Sun, W.; Zhang, L.; Li, X.; Wang, L.; Jiang, L.; Zheng, Y. Characterisation of flower colouration in 30 Rhododendron species via anthocyanin and flavonol identification and quantitative traits. Plant Biol. 2018, 20, 121–129. [Google Scholar] [CrossRef]
  26. Sun, W.; Zhou, N.; Wang, Y.; Sun, S.; Zhang, Y.; Ju, Z.; Yi, Y. Characterization and functional analysis of RdDFR1 regulation on flower color formation in Rhododendron delavayi. Plant Physiol. Biochem. 2021, 169, 203–210. [Google Scholar] [CrossRef] [PubMed]
  27. Sun, W.; Sun, S.; Xu, H.; Wang, Y.; Chen, Y.; Xu, X.; Yi, Y.; Ju, Z. Characterization of two key flavonoid 3-O-glycosyltransferases involved in the formation of flower color in Rhododendron delavayi. Front. Plant Sci. 2022, 13, 863482. [Google Scholar] [CrossRef]
  28. Zhang, L.; Xu, P.; Cai, Y.; Ma, L.; Li, S.; Li, S.; Xie, W.; Song, J.; Peng, L.; Yan, H.; et al. The draft genome assembly of Rhododendron delavayi Franch. var. delavayi. Gigascience 2017, 6, 10. [Google Scholar] [CrossRef] [PubMed]
  29. Kuehni, R.G. Color-tolerance data and the tentative CIE 1976 L a b formula. J. Opt. Soc. Am. 1976, 66, 497–500. [Google Scholar] [CrossRef]
  30. Zhao, Y.; Li, A.; Qi, S.; Su, K.; Guo, Y. Identification of candidate genes related to anthocyanin biosynthesis in red sarcocarp hawthorn (Crataegus pinnatifida). Sci. Hortic. 2022, 298, 110987. [Google Scholar] [CrossRef]
  31. Lin, H.-H. A systematic approach of predicting color preference on the basis of gray relational grade. Color Res. Appl. 2019, 44, 194–204. [Google Scholar] [CrossRef]
  32. Katiyar, A.; Smita, S.; Lenka, S.K.; Rajwanshi, R.; Chinnusamy, V.; Bansal, K.C. Genome-wide classification and expression analysis of MYB transcription factor families in rice and Arabidopsis. BMC Genom. 2012, 13, 544. [Google Scholar] [CrossRef]
  33. Wong, G.R.; Latif, S.N.F.B.A.; Mazumdar, P. Genome-wide investigation and comparative expression profiling reveal R2R3-MYB genes involved in sclerotinia defence in tomato. Physiol. Mol. Plant Pathol. 2022, 121, 101873. [Google Scholar] [CrossRef]
  34. Chen, G.; He, W.; Guo, X.; Pan, J. Genome-wide identification, classification and expression analysis of the MYB transcription factor family in Petunia. Int. J. Mol. Sci. 2021, 22, 4838. [Google Scholar] [CrossRef] [PubMed]
  35. Li, W.; Ding, Z.; Ruan, M.; Yu, X.; Peng, M.; Liu, Y. Kiwifruit R2R3-MYB transcription factors and contribution of the novel AcMYB75 to red kiwifruit anthocyanin biosynthesis. Sci. Rep. 2017, 7, 16861. [Google Scholar] [CrossRef] [PubMed]
  36. Brendolise, C.; Espley, R.V.; Lin-Wang, K.; Laing, W.; Peng, Y.; McGhie, T.; Dejnoprat, S.; Tomes, S.; Hellens, R.P.; Allan, A.C. Multiple copies of a simple MYB-binding site confers trans-regulation by specific flavonoid-related R2R3 MYBs in Diverse Species. Front. Plant Sci. 2017, 8, 1864. [Google Scholar] [CrossRef] [PubMed]
  37. Cao, Z.-H.; Zhang, S.-Z.; Wang, R.-K.; Zhang, R.-F.; Hao, Y.-J. Genome wide analysis of the apple MYB transcription factor Family Allows the Identification of MdoMYB121 Gene Confering Abiotic Stress Tolerance in Plants. PLoS ONE 2013, 8, 7. [Google Scholar] [CrossRef] [PubMed]
  38. Zhang, C.; Ma, R.; Xu, J.; Yan, J.; Guo, L.; Song, J.; Feng, R.; Yu, M. Genome-wide identification and classification of MYB superfamily genes in peach. PLoS ONE 2018, 13, 6. [Google Scholar] [CrossRef] [PubMed]
  39. Feller, A.; Machemer, K.; Braun, E.L.; Grotewold, E. Evolutionary and comparative analysis of MYB and bHLH plant transcription factors. Plant J. 2011, 66, 94–116. [Google Scholar] [CrossRef]
  40. Kang, L.; Teng, Y.; Cen, Q.; Fang, Y.; Tian, Q.; Zhang, X.; Wang, H.; Zhang, X.; Xue, D. Genome-wide identification of R2R3-MYB transcription factor and expression analysis under abiotic stress in rice. Plants 2022, 11, 1928. [Google Scholar] [CrossRef]
  41. Zhu, K.; Fan, P.; Mo, Z.; Tan, P.; Feng, G.; Li, F.; Peng, F. Identification, expression and co-expression analysis of R2R3-MYB family genes involved in graft union formation in pecan (Carya illinoinensis). Forests 2020, 11, 917. [Google Scholar] [CrossRef]
  42. Sabir, I.A.; Manzoor, M.A.; Shah, I.H.; Liu, X.; Zahid, M.S.; Jiu, S.; Wang, J.; Abdullah, M.; Zhang, C. MYB transcription factor family in sweet cherry (Prunus avium L.): Genome-wide investigation, evolution, structure, characterization and expression patterns. BMC Plant Biol. 2022, 22, 1. [Google Scholar] [CrossRef]
  43. Rosinski, J.A.; Atchley, W.R. Molecular evolution of the MYB family of transcription factors: Evidence for polyphyletic origin. J. Mol. Evol. 1998, 46, 74–83. [Google Scholar] [CrossRef] [PubMed]
  44. Li, J.; He, Y.-J.; Zhou, L.; Liu, Y.; Jiang, M.; Ren, L.; Chen, H. Transcriptome profiling of genes related to light-induced anthocyanin biosynthesis in eggplant (Solanum melongena L.) before purple color becomes evident. BMC Genom. 2018, 19, 201. [Google Scholar] [CrossRef]
  45. Ye, L.-J.; Moller, M.; Luo, Y.-H.; Zou, J.-Y.; Zheng, W.; Wang, Y.-H.; Liu, J.; Zhu, A.-D.; Hu, J.-Y.; Li, D.-Z.; et al. Differential expressions of anthocyanin synthesis genes underlie flower color divergence in a sympatric Rhododendron sanguineum complex. BMC Plant Biol. 2021, 21, 1. [Google Scholar] [CrossRef]
  46. Liu, H.; Xiong, J.-s.; Jiang, Y.-t.; Wang, L.; Cheng, Z.-m. Evolution of the R2R3-MYB gene family in six rosaceae species and expression in woodland strawberry. J. Integr. Agric. 2019, 18, 2753–2770. [Google Scholar] [CrossRef]
  47. Yu, M.; Man, Y.; Wang, Y. Light- and temperature-induced expression of an R2R3-MYB gene regulates anthocyanin biosynthesis in red-fleshed kiwifruit. Int. J. Mol. Sci. 2019, 20, 5228. [Google Scholar] [CrossRef] [PubMed]
  48. Yin, X.; Wang, T.; Zhang, M.; Zhang, Y.; Irfan, M.; Chen, L.; Zhang, L. Role of core structural genes for flavonoid biosynthesis and transcriptional factors in flower color of plants. Biotechnol. Biotechnol. Equip. 2021, 35, 1214–1229. [Google Scholar] [CrossRef]
  49. Li, Y.; Shan, X.; Gao, R.; Han, T.; Zhang, J.; Wang, Y.; Kimani, S.; Wang, L.; Gao, X. MYB repressors and MBW activation complex collaborate to fine-tune flower coloration in Freesia hybrida. Commun. Biol. 2020, 3, 1. [Google Scholar] [CrossRef] [PubMed]
  50. Meng, J.; Gao, Y.; Han, M.; Liu, P.; Yang, C.; Shen, T.; Li, H. In vitro anthocyanin induction and metabolite analysis in Malus spectabilis leaves under low nitrogen conditions. Hortic. Plant J. 2020, 6, 284–292. [Google Scholar] [CrossRef]
  51. Fang, Z.-Z.; Kui, L.-W.; Zhou, D.-R.; Lin, Y.-J.; Jiang, C.-C.; Pan, S.-L.; Espley, R.V.; Andre, C.M.; Ye, X.-F. Activation of PsMYB10.2 transcription causes anthocyanin accumulation in flesh of the red-fleshed mutant of ‘Sanyueli’ (Prunus salicina Lindl.). Front. Plant Sci. 2021, 12, 1167. [Google Scholar] [CrossRef]
  52. Abid, M.A.; Wei, Y.; Meng, Z.; Wang, Y.; Ye, Y.; Wang, Y.; He, H.; Zhou, Q.; Li, Y.; Wang, P.; et al. Increasing floral visitation and hybrid seed production mediated by beauty mark in Gossypium hirsutum. Plant Biotechnol. J. 2022, 20, 1274–1284. [Google Scholar] [CrossRef] [PubMed]
  53. Yu, M.; Man, Y.; Lei, R.; Lu, X.; Wang, Y. Metabolomics study of flavonoids and anthocyanin-related gene analysis in kiwifruit (Actinidia chinensis) and kiwiberry (Actinidia arguta). Plant Mol. Biol. Rep. 2020, 38, 353–369. [Google Scholar] [CrossRef]
  54. Xu, Z.-S.; Yang, Q.-Q.; Feng, K.; Yu, X.; Xiong, A.-S. DcMYB113, a root-specific R2R3-MYB, conditions anthocyanin biosynthesis and modification in carrot. Plant Biotechnol. J. 2020, 18, 1585–1597. [Google Scholar] [CrossRef] [PubMed]
  55. Xi, W.; Feng, J.; Liu, Y.; Zhang, S.; Zhao, G. The R2R3-MYB transcription factor PaMYB10 is involved in anthocyanin biosynthesis in apricots and determines red blushed skin. BMC Plant Biol. 2019, 19, 287. [Google Scholar] [CrossRef] [PubMed]
  56. Yoon, B.-J. Hidden markov models and their applications in biological sequence analysis. Curr. Genom. 2009, 10, 402–415. [Google Scholar] [CrossRef]
  57. Sigrist, C.J.A.; Cerutti, L.; de Castro, E.; Langendijk-Genevaux, P.S.; Bulliard, V.; Bairoch, A.; Hulo, N. PROSITE, a protein domain database for functional characterization and annotation. Nucleic Acids Res. 2010, 38, D161–D166. [Google Scholar] [CrossRef]
  58. Bernstein, M.N.; Ma, Z.; Gleicher, M.; Dewey, C.N. CellO: Comprehensive and hierarchical cell type classification of human cells with the cell ontology. iScience 2021, 24, 1. [Google Scholar] [CrossRef]
  59. Chen, C.; Chen, H.; Zhang, Y.; Thomas, H.R.; Frank, M.H.; He, Y.; Xia, R. TBtools: An integrative toolkit developed for interactive analyses of big biological data. Mol. Plant 2020, 13, 1194–1202. [Google Scholar] [CrossRef]
  60. Bailey, T.L.; Johnson, J.; Grant, C.E.; Noble, W.S. The MEME suite. Nucleic Acids Res. 2015, 43, W39–W49. [Google Scholar] [CrossRef]
  61. Crooks, G.E.; Hon, G.; Chandonia, J.-M.; Brenner, S.E. WebLogo: A sequence logo generator. Genome Res. 2004, 14, 1188–1190. [Google Scholar] [CrossRef]
  62. Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 2016, 33, 1870–1874. [Google Scholar] [CrossRef] [PubMed]
  63. Chow, C.-N.; Lee, T.-Y.; Hung, Y.-C.; Li, G.-Z.; Tseng, K.-C.; Liu, Y.-H.; Kuo, P.-L.; Zheng, H.-Q.; Chang, W.-C. PlantPAN3.0: A new and updated resource for reconstructing transcriptional regulatory networks from ChIP-seq experiments in plants. Nucleic Acids Res. 2019, 47, D1155–D1163. [Google Scholar] [CrossRef] [PubMed]
  64. Rombauts, S.; Dehais, P.; Van Montagu, M.; Rouze, P. PlantCARE, a plant cis-acting regulatory element database. Nucleic Acids Res. 1999, 27, 295–296. [Google Scholar] [CrossRef] [PubMed]
  65. Liu, N.; Zhang, L.; Zhou, Y.; Tu, M.; Wu, Z.; Gui, D.; Ma, Y.; Wang, J.; Zhang, C. The Rhododendron Plant Genome Database (RPGD): A comprehensive online omics database for Rhododendron. BMC Genom. 2021, 22, 376. [Google Scholar] [CrossRef]
  66. Smith, T.; Heger, A.; Sudbery, I. UMI-tools: Modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res. 2017, 27, 491–499. [Google Scholar] [CrossRef]
  67. Kim, D.; Paggi, J.M.; Park, C.; Bennett, C.; Salzberg, S.L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 2019, 37, 8. [Google Scholar] [CrossRef]
  68. Liao, Y.; Smyth, G.K.; Shi, W. Feature counts: An efficient general-purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930. [Google Scholar] [CrossRef]
  69. Lopes, C.T.; Franz, M.; Kazi, F.; Donaldson, S.L.; Morris, Q.; Bader, G.D. Cytoscape Web: An interactive web-based network browser. Bioinformatics 2010, 26, 2347–2348. [Google Scholar] [CrossRef]
  70. Wang, Y.; Tian, R.M.; Gao, Z.M.; Bougouffa, S.; Qian, P.-Y. Optimal eukaryotic 18S and universal 16S/18S ribosomal RNA primers and their application in a study of symbiosis. PLoS ONE 2014, 9, 3. [Google Scholar] [CrossRef]
Ijms 24 04641 g001 550
Figure 1. Conserved MYB domain in MYBs of Rhododendron delavayi. (A) Phylogenetic tree trunk of 184 R. delavayi MYB genes. It was constructed by maximum likelihood (ML) method based on multiple sequence alignment of 184 amino acid sequences of R. delavayi MYB gene. Different colored circles indicate each node of the phylogenetic tree. (B) Branches on each node of the trunk of the phylogenetic tree. The MYB domain was divided into 1R-MYB, R2R3-MYB, 3R-MYB, and 4R-MYB and marked with a colored background. (C) Conserved domain composition of these genes. The colored squares indicate the different structural domains.
Figure 1. Conserved MYB domain in MYBs of Rhododendron delavayi. (A) Phylogenetic tree trunk of 184 R. delavayi MYB genes. It was constructed by maximum likelihood (ML) method based on multiple sequence alignment of 184 amino acid sequences of R. delavayi MYB gene. Different colored circles indicate each node of the phylogenetic tree. (B) Branches on each node of the trunk of the phylogenetic tree. The MYB domain was divided into 1R-MYB, R2R3-MYB, 3R-MYB, and 4R-MYB and marked with a colored background. (C) Conserved domain composition of these genes. The colored squares indicate the different structural domains.
Ijms 24 04641 g001
Ijms 24 04641 g002 550
Figure 2. R2 and R3 MYB repeats of R2R3-MYBs are highly conserved in Rhododendron delavayi. R2 (A) and R3 (B) MYB repeat sequences were based on multiple alignment analysis of the R2R3-MYB structural domain of all R. delavayi. Bits indicate the relative frequency of amino acids. Asterisks indicate conserved Trp residues in the MYB domain. The amino acids are arranged from the N terminal to the C terminal.
Figure 2. R2 and R3 MYB repeats of R2R3-MYBs are highly conserved in Rhododendron delavayi. R2 (A) and R3 (B) MYB repeat sequences were based on multiple alignment analysis of the R2R3-MYB structural domain of all R. delavayi. Bits indicate the relative frequency of amino acids. Asterisks indicate conserved Trp residues in the MYB domain. The amino acids are arranged from the N terminal to the C terminal.
Ijms 24 04641 g002
Ijms 24 04641 g003 550
Figure 3. Motif composition and gene structure of Rhododendron delavayi MYB genes. (A) Phylogenetic tree trunk of 184 R. delavayi MYB genes. Different colored circles indicate each node of the phylogenetic tree. (B) The MYB structural domains on each node branch were divided and marked with colored backgrounds. (C) Conservation motifs of the MYB proteins. Each motif represents a number in the colored box, and a black line represents the non-conserved sequence. (D) Gene structure of MYB. Green, yellow, and black colors represent the UTR, exon, and intron, respectively. The scale bar at the bottom indicates the length of motifs and genes, respectively.
Figure 3. Motif composition and gene structure of Rhododendron delavayi MYB genes. (A) Phylogenetic tree trunk of 184 R. delavayi MYB genes. Different colored circles indicate each node of the phylogenetic tree. (B) The MYB structural domains on each node branch were divided and marked with colored backgrounds. (C) Conservation motifs of the MYB proteins. Each motif represents a number in the colored box, and a black line represents the non-conserved sequence. (D) Gene structure of MYB. Green, yellow, and black colors represent the UTR, exon, and intron, respectively. The scale bar at the bottom indicates the length of motifs and genes, respectively.
Ijms 24 04641 g003
Ijms 24 04641 g004 550
Figure 4. Phylogenetic tree of MYB proteins between Arabidopsis thaliana and Rhododendron delavayi (S). ML phylogeny was determined by MEGA7.0. Phylogenetic tree was divided into 35 main groups. Different colored clusters are labeled outside the tree, whereas their subgroups are marked with the corresponding color backgrounds. The blue circles on the branches are bootstrap and the numbers are branch lengths.
Figure 4. Phylogenetic tree of MYB proteins between Arabidopsis thaliana and Rhododendron delavayi (S). ML phylogeny was determined by MEGA7.0. Phylogenetic tree was divided into 35 main groups. Different colored clusters are labeled outside the tree, whereas their subgroups are marked with the corresponding color backgrounds. The blue circles on the branches are bootstrap and the numbers are branch lengths.
Ijms 24 04641 g004
Ijms 24 04641 g005a 550Ijms 24 04641 g005b 550
Figure 5. Cis-acting elements in the MYB promoter region of Rhododendron delavayi. (A) MYB gene promoter region regulatory elements. Boxes of different colors represent various cis elements and are distributed in different positions among 184 MYB genes. Axes represent gene length. (B) Statistics of cis elements in MYB promoter involved in plant growth and development, phytohormone response, and abiotic stress response.
Figure 5. Cis-acting elements in the MYB promoter region of Rhododendron delavayi. (A) MYB gene promoter region regulatory elements. Boxes of different colors represent various cis elements and are distributed in different positions among 184 MYB genes. Axes represent gene length. (B) Statistics of cis elements in MYB promoter involved in plant growth and development, phytohormone response, and abiotic stress response.
Ijms 24 04641 g005aIjms 24 04641 g005b
Ijms 24 04641 g006 550
Figure 6. Basic transcriptome data of five types of red samples of R. delavayi spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex. (A) Phenotypes of five types of samples of R. delavayi: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The enclosed part represents the sampling site. (B) PCA analysis of the five types of samples; (C) Pearson correlation coefficient of the five types of samples.
Figure 6. Basic transcriptome data of five types of red samples of R. delavayi spotted petals, unspotted petals, spotted throat, unspotted throat, and branchlet cortex. (A) Phenotypes of five types of samples of R. delavayi: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The enclosed part represents the sampling site. (B) PCA analysis of the five types of samples; (C) Pearson correlation coefficient of the five types of samples.
Ijms 24 04641 g006
Ijms 24 04641 g007 550
Figure 7. RNA-Seq analysis of R. delavayi. (A) Histograms of DEGs obtained by comparison of four pairs of samples, namely, MY-1 vs MY-2, MY-3 vs MY-4, MY-5 vs MY-2, and MY-5 vs MY-4, respectively. (B) Venn diagram of DEGs in the four combinations. (C). KEGG enrichment analysis of MY-1 vs MY-2 DEGs.
Figure 7. RNA-Seq analysis of R. delavayi. (A) Histograms of DEGs obtained by comparison of four pairs of samples, namely, MY-1 vs MY-2, MY-3 vs MY-4, MY-5 vs MY-2, and MY-5 vs MY-4, respectively. (B) Venn diagram of DEGs in the four combinations. (C). KEGG enrichment analysis of MY-1 vs MY-2 DEGs.
Ijms 24 04641 g007
Ijms 24 04641 g008 550
Figure 8. Expression heat maps of 86 R2R3-MYB genes in different samples of R. delavayi: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The genes on the right indicate the differentially expressed R2R3-MYB in the five types of red samples. Blue rectangles indicate low expression levels, whereas red rectangles indicate high expression levels. Generated hierarchical clustering tree maps and heat maps based on log2 values of FPKM.
Figure 8. Expression heat maps of 86 R2R3-MYB genes in different samples of R. delavayi: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The genes on the right indicate the differentially expressed R2R3-MYB in the five types of red samples. Blue rectangles indicate low expression levels, whereas red rectangles indicate high expression levels. Generated hierarchical clustering tree maps and heat maps based on log2 values of FPKM.
Ijms 24 04641 g008
Ijms 24 04641 g009 550
Figure 9. Correlation analysis of WGCNA module with chromatic aberration. Each column corresponds to L*, a*, and b* values, and each row represents a module. Chromatic aberration values and module correlation coefficients are shown within each module. p-values are shown in parentheses.
Figure 9. Correlation analysis of WGCNA module with chromatic aberration. Each column corresponds to L*, a*, and b* values, and each row represents a module. Chromatic aberration values and module correlation coefficients are shown within each module. p-values are shown in parentheses.
Ijms 24 04641 g009
Ijms 24 04641 g010 550
Figure 10. Visualization of hub genes in the MEturquoise module related to the a* value. Different color lumps represent distinct families, with a larger shape for a greater degree and a darker line for a greater weight.
Figure 10. Visualization of hub genes in the MEturquoise module related to the a* value. Different color lumps represent distinct families, with a larger shape for a greater degree and a darker line for a greater weight.
Ijms 24 04641 g010
Ijms 24 04641 g011 550
Figure 11. qRT-PCR verification of the transcriptome result. Expression levels of nine genes in different tissues: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The 18s gene was applied as a reference gene. Relative gene expression was calculated with MY-5 as control. Each column represents the mean ± standard error of the expression of three independent experiments with three replications.
Figure 11. qRT-PCR verification of the transcriptome result. Expression levels of nine genes in different tissues: MY-1, spotted petals; MY-2, unspotted petals; MY-3, spotted throat; MY-4, unspotted throat; MY-5, branchlet cortex. The 18s gene was applied as a reference gene. Relative gene expression was calculated with MY-5 as control. Each column represents the mean ± standard error of the expression of three independent experiments with three replications.
Ijms 24 04641 g011
Table
Table 1. Chromatic aberration values of different tissues of R. delavayi.
Table 1. Chromatic aberration values of different tissues of R. delavayi.
SamplesL*a*b*
MY-132.26 ± 1.37 b57.90 ± 3.98 b4.37 ± 1.39 b
MY-236.25 ± 0.29 a64.23 ± 2.27 a7.86 ± 0.61 a
MY-329.28 ± 1.53 c39.83 ± 4.39 c1.99 ± 1.12 c
MY-437.79 ± 1.13 a53.02 ± 1.10 b3.18 ± 0.76 bc
MY-532.08 ± 1.14 b10.23 ± 2.74 d8.58 ± 1.30 a
Significant differences (p < 0.05) are indicated by different letters.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Long, F.; Wu, H.; Li, H.; Zuo, W.; Ao, Q. Genome-Wide Analysis of MYB Transcription Factors and Screening of MYBs Involved in the Red Color Formation in Rhododendron delavayi. Int. J. Mol. Sci. 2023, 24, 4641. https://doi.org/10.3390/ijms24054641

AMA Style

Long F, Wu H, Li H, Zuo W, Ao Q. Genome-Wide Analysis of MYB Transcription Factors and Screening of MYBs Involved in the Red Color Formation in Rhododendron delavayi. International Journal of Molecular Sciences. 2023; 24(5):4641. https://doi.org/10.3390/ijms24054641

Chicago/Turabian Style

Long, Fenfang, Hairong Wu, Huie Li, Weiwei Zuo, and Qian Ao. 2023. "Genome-Wide Analysis of MYB Transcription Factors and Screening of MYBs Involved in the Red Color Formation in Rhododendron delavayi" International Journal of Molecular Sciences 24, no. 5: 4641. https://doi.org/10.3390/ijms24054641

APA Style

Long, F., Wu, H., Li, H., Zuo, W., & Ao, Q. (2023). Genome-Wide Analysis of MYB Transcription Factors and Screening of MYBs Involved in the Red Color Formation in Rhododendron delavayi. International Journal of Molecular Sciences, 24(5), 4641. https://doi.org/10.3390/ijms24054641

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop