Genome-Wide Identification, Classification and Expression Analysis of the MYB Transcription Factor Family in Petunia

A lot of researches have been focused on the evolution and function of MYB transcription factors (TFs). For revealing the formation of petunia flower color diversity, MYB gene family in petunia was identified and analyzed. In this study, a total of 155 MYB genes, including 40 1R-MYBs, 106 R2R3-MYBs, 7 R1R2R3-MYBs and 2 4R-MYBs, have been identified in the Petunia axillaris genome. Most R2R3 genes contain three exons and two introns, whereas the number of PaMYB introns varies from 0 to 12. The R2R3-MYB members could be divided into 28 subgroups. Analysis of gene structure and protein motifs revealed that members within the same subgroup presented similar exon/intron and motif organization, further supporting the results of phylogenetic analysis. Genes in subgroup 10, 11 and 21 were mainly expressed in petal, not in vegetative tissues. Genes in subgroup 9, 19, 25 and 27 expressed in all tissues, but the expression patterns of each gene were different. According to the promoter analysis, five R2R3-MYB and two MYB-related genes contained MBSI cis-element, which was involved in flavonoid biosynthetic regulation. PaMYB100/DPL has been reported to positively regulate to pigmentation. However, although PaMYB82, PaMYB68 and Pa1RMYB36 contained MBSI cis-element, their function in flavonoid biosynthesis has not been revealed. Consistent with existing knowledge, PaMYBs in subgroup 11 had similar function to AtMYBs in subgroup 6, genes in which played an important role in anthocyanin biosynthesis. In addition, PaMYB1 and PaMYB40 belonged to P9 (S7) and were potentially involved in regulation of flavonoid synthesis in petunia vegetative organs. This work provides a comprehensive understanding of the MYB gene family in petunia and lays a significant foundation for future studies on the function and evolution of MYB genes in petunia.


Introduction
As one of the largest transcription factor (TFs) families in plants, MYB genes play an important role in plant biology, including development, environmental adaptation and metabolic regulation [1][2][3]. The first plant MYB-encoding gene was isolated from maize (Zea mays), named as COLORED1 (C1), indicating its function in anthocyanin biosynthesis [4]. Subsequently, more and more MYB genes were identified in plants, and the MYB superfamily is dramatically expanded in plants [5]. Currently, 22,032 MYB and 15,369 MYB-related sequences are available in the plant TF database (http://planttfdb.gaolab.org/, accessed on 1 February 2021) [6].
MYB transcription factor harbors a conserved MYB DNA binding domain (DBD) at the N terminus and a diverse C-terminal modulator region responsible for the regulatory activity of the protein. The MYB domain frequently consists of up to four imperfect repeats, namely 1R-MYB (MYB-related), R2R3-MYB, R1R2R3-MYB and 4R-MYB based on the number of MYB repeats [7]. MYB subfamily containing two repeats, characterized by R2R3-type MYB proteins, is the largest within the MYB family in plants [8,9]. In general, point (PI) and subcellular localizations of all PaMYB genes are shown in supplementary Tables S1 and S2. Peaxi162Scf00110g00004 was the smallest protein with 73 amino acids, because it obtained only one DBD domain. Peaxi162Scf00523g00004 was the largest protein with 1051 amino acids. Consequently, the protein MWs varied from 8.63 kDa (Peaxi162Scf00110g00004) to 116.60 kDa (Peaxi162Scf00523g00004), and the pIs ranged from 4.66 (Peaxi162Scf00944g00140) to 10.02 (Peaxi162Scf00118g00310). The subcellular localization of all the PaMYB proteins was the nucleus, indicating that PaMYBs were nucleoproteins. The sequence similarity between PaMYBs and AtMYBs ranged from 38% to 99%.

Analysis of the Gene Structure and Conserved MYB Domains in PaMYBs
According to the result of pre-screening, the sequence of 1R-MYB showed lowest similarity to other types of MYB genes. Thus, the analysis of gene structure and conserved domain was separated into two parts, the 1R-MYB part and the other MYBs part, respectively ( Figure 1 and Supplementary Figure S1).
In total, there are 10 conserved motifs identified in part of 2R-,3R-and 4R-MYB and part of 1R-MYB by MEME, respectively (Supplementary Tables S3 and S4). In order to better understanding the evolutionary relationships of MYB proteins in Petunia, an unrooted phylogenetic tree was constructed ( Figure 1A and supplementary Figure S1A). Only motif 3 existed in every 2R-, 3R-and 4R-MYB protein. Besides this, most of the genes with close evolutionary relationships shared the same motif compositions, suggesting functional similarity of MYB proteins in the same group. To investigate the MYB DBD domain sequence features of petunia, the frequency of amino acid sequences of R2-and R3-repeats of 106 R2R3-MYB proteins was performed by WebLogo and multiple alignment analysis ( Figure 2). R2-repeats contained motif 3, 5 and 1, and R3-repeats contained motif 4 and 2. In general, the basic regions of MYB domains had 108 basic residues, including the linker region. The R2 and R3 MYB repeats of the PaR2R3-MYBs contained characteristic amino acids, including three conserved Trp (W) residues identified in R2-repeat, while only the second and third Trp (W) were conserved in R3-repeat. The first Trp (W) in R3-repeat was usually replaced by F (69%)\I (13%) in most R2R3-MYBs. In R2-repeat, a total of 29 (out of 54) positions were occupied by a single residue in >70% of the proteins. For R3-repeat, the equivalent number was 30 (out of 54), indicating that R3-repeat was slightly more conserved than R2-repeat. The linker region between R2-and R3-repeats was highly conserved, and that four amino acids in the first half of the linker (LRPD, 53-56 amino acid) formed a highly conserved motif [38]. It is noteworthy that a similar motif exists in the linker region in tomato MYB proteins as well [39].
To understand the structure of PaMYB genes, the exon and intron organizations were studied by comparing the complementary DNA (cDNA) sequences with the corresponding genomic DNA sequences ( Figure 1C and Supplementary Figure S1C). According to the results of gene structure analysis, the number of introns varies from 0 to 12, and most genes in the same subgroup have similar exon/intron structures. Subgroup P1 possessed no introns and all the members in P5 had multiple introns. Most of the R2R3-MYBs in petunia were typical splicing of three exons and two introns (73 of the 106 R2R3-MYBs, 67%). Although the lengths of intron and exon are different, the structure of R2R3-MYB in the same subfamily is highly conservative. On average, the 1R-MYBs were disrupted by 2.1 introns, 2R-MYBs were disrupted by 1.97 introns, 3R-MYBs were disrupted by 6.7 introns and 4R-MYBs were disrupted by 4 introns.

Analysis of the Gene Structure and Conserved MYB Domains in PaMYBs
According to the result of pre-screening, the sequence of 1R-MYB showed lowest similarity to other types of MYB genes. Thus, the analysis of gene structure and conserved domain was separated into two parts, the 1R-MYB part and the other MYBs part, respectively ( Figure 1 and Supplementary Figure S1).  The phylogenetic tree was constructed with MEGA 10.1.8 using the neighbor-joining (NJ) method with 1000 bootstrap replicates based on a multiple alignment of amino acid sequences of MYB genes. The subgroup was classified and marked by alternated dark and light gray; (B) Protein motif. Schematic diagram of the conserved motifs in the MYB proteins, which were elucidated using MEME. Each motif is represented by a number in the colored box. The black lines represent the non-conserved sequences; (C) Gene structures of PaMYBs. The exons are represented by green box, the black lines connecting two exons represent introns and untranslated regions are indicated by yellow box.

Phylogenetic Analysis and Classification of R2R3-MYB Proteins in Petunia
R2R3-MYB subfamily occupied the largest proportion with 68.39% of the total petunia MYB TF family. In order to better understanding the feature of R2R3-MYB in Petunia, an unrooted phylogenetic tree was constructed depending on the full-length sequences of R2R3-MYB proteins from Petunia and Arabidopsis (Figure 3 and Supplementary Table S5). Relying on the group classification of R2R3-MYB in Arabidopsis, the PaR2R3-MYBs were subdivided into 28 subgroups (designated P1-P28 in this study) based on the topology of the tree [40]. Most petunia proteins were found to be shared among the R2R3-MYBs from Arabidopsis. In contrast, seven clades (P6, P10, P13, P22, P23, P26, P28) were not clustered together with Arabidopsis MYBs and several clades were found only containing one species (P11, P18).
repeat, a total of 29 (out of 54) positions were occupied by a single residue in >70% of the proteins. For R3-repeat, the equivalent number was 30 (out of 54), indicating that R3-repeat was slightly more conserved than R2-repeat. The linker region between R2-and R3repeats was highly conserved, and that four amino acids in the first half of the linker (LRPD, 53-56 amino acid) formed a highly conserved motif [38]. It is noteworthy that a similar motif exists in the linker region in tomato MYB proteins as well [39]. To understand the structure of PaMYB genes, the exon and intron organizations were studied by comparing the complementary DNA (cDNA) sequences with the corresponding genomic DNA sequences ( Figure 1C and Supplementary Figure S1C). According to the results of gene structure analysis, the number of introns varies from 0 to 12, and most genes in the same subgroup have similar exon/intron structures. Subgroup P1 possessed no introns and all the members in P5 had multiple introns. Most of the R2R3-MYBs in petunia were typical splicing of three exons and two introns (73 of the 106 R2R3-MYBs, 67%). Although the lengths of intron and exon are different, the structure of R2R3-MYB in the same subfamily is highly conservative. On average, the 1R-MYBs were disrupted by 2.1 introns, 2R-MYBs were disrupted by 1.97 introns, 3R-MYBs were disrupted by 6.7 introns and 4R-MYBs were disrupted by 4 introns.

Phylogenetic Analysis and Classification of R2R3-MYB Proteins in Petunia
R2R3-MYB subfamily occupied the largest proportion with 68.39% of the total petunia MYB TF family. In order to better understanding the feature of R2R3-MYB in Petunia, an unrooted phylogenetic tree was constructed depending on the full-length sequences of R2R3-MYB proteins from Petunia and Arabidopsis (Figure 3 and Supplementary Table S5). Relying on the group classification of R2R3-MYB in Arabidopsis, the PaR2R3-MYBs were subdivided into 28 subgroups (designated P1-P28 in this study) based on the topology of the tree [40]. Most petunia proteins were found to be shared among the R2R3-MYBs from Arabidopsis. In contrast, seven clades (P6, P10, P13, P22, P23, P26, P28) were not clustered together with Arabidopsis MYBs and several clades were found only containing one species (P11, P18).

Promoter Analysis of PaMYB Genes
For exploring the transcriptional regulation characteristics of PaMYBs, cis-elements were predicted using PlantCARE and marked on promoter regions ( Figure 4). The hormone and stress responsive elements were widely distributed in PaMYBs promoter. By

Promoter Analysis of PaMYB Genes
For exploring the transcriptional regulation characteristics of PaMYBs, cis-elements were predicted using PlantCARE and marked on promoter regions ( Figure 4). The hormone and stress responsive elements were widely distributed in PaMYBs promoter. By calculating the number of different cis-elements, MeJA-responsive element (CGTCA-motif and TGACG-motif) was the most frequent in PaMYBs promoter, followed by abscisic acid-responsive element (ABRE) and anaerobic induction (ARE) ( Table 1). Notably, the common cis-elements in 1R-and R2R3-MYB promoters were significantly different to that in 3R-and 4R-MYB promoters. For example, flavonoid biosynthetic regulation (MBSI) and wound-responsive elements (WUN-motif) only existed in 1R-and R2R3-MYB promoters. None of 4R-MYB genes contained low temperature-responsive element (LTR) and salicylic acid-responsive element (TCA-element and SARE). This result suggested that transcriptional regulation of different types of PaMYB gene was diverse, indicating the diversity of PaMYBs function. calculating the number of different cis-elements, MeJA-responsive element (CGTCA-motif and TGACG-motif) was the most frequent in PaMYBs promoter, followed by abscisic acid-responsive element (ABRE) and anaerobic induction (ARE) ( Table 1). Notably, the common cis-elements in 1R-and R2R3-MYB promoters were significantly different to that in 3R-and 4R-MYB promoters. For example, flavonoid biosynthetic regulation (MBSI) and wound-responsive elements (WUN-motif) only existed in 1R-and R2R3-MYB promoters. None of 4R-MYB genes contained low temperature-responsive element (LTR) and salicylic acid-responsive element (TCA-element and SARE). This result suggested that transcriptional regulation of different types of PaMYB gene was diverse, indicating the diversity of PaMYBs function.

Expression of PaMYB Genes in Different Tissues
Based on RNA-seq database (unpublished data) owned by our lab, the expressions of all the PaMYB genes in petal, leaf and root were extracted (supplementary Table S5). A percentage of 32.9% of PaMYB genes showed no expression or were undetected. Thirtytwo out of 155 PaMYB genes expressed in all the three tissues. Genes in subgroup 10, 11 and 21 were mainly expressed in petal, not in vegetative tissues. Genes in subgroup 9, 19, 25 and 27 expressed in all tissues, but the expression patterns of each gene were different. Among these three tissues, there are 27, 2 and 17 genes only expressed in petal, leaf and root, respectively. This may be caused by the different genetic background between varieties of sequencing and used in this work. To explore the specific functions of different subgroup of R2R3-MYB, we selected several PaMYB genes from each subgroup and determined their expression levels in petal, leaf and root by quantitative real-time polymerase chain reaction (qRT-PCR) ( Figure 5). The expression patterns of each gene were different. For example, PaMYB75 had the highest expression level in petal compared to leaf and root, while PaMYB52 was mostly expressed in leaf. On the other hand, the transcripts of PaMYB74, PaMYB59 and several genes were fewer in flower. This indicated that MYB genes in different subgroups may perform different functions during the growth and development of petunia.

Expression of PaMYB Genes in Different Tissues
Based on RNA-seq database (unpublished data) owned by our lab, the expressions of all the PaMYB genes in petal, leaf and root were extracted (supplementary Table S5). A percentage of 32.9% of PaMYB genes showed no expression or were undetected. Thirtytwo out of 155 PaMYB genes expressed in all the three tissues. Genes in subgroup 10, 11 and 21 were mainly expressed in petal, not in vegetative tissues. Genes in subgroup 9, 19, 25 and 27 expressed in all tissues, but the expression patterns of each gene were different. Among these three tissues, there are 27, 2 and 17 genes only expressed in petal, leaf and root, respectively. This may be caused by the different genetic background between varieties of sequencing and used in this work. To explore the specific functions of different subgroup of R2R3-MYB, we selected several PaMYB genes from each subgroup and determined their expression levels in petal, leaf and root by quantitative real-time polymerase chain reaction (qRT-PCR) ( Figure 5). The expression patterns of each gene were different. For example, PaMYB75 had the highest expression level in petal compared to leaf and root, while PaMYB52 was mostly expressed in leaf. On the other hand, the transcripts of PaMYB74, PaMYB59 and several genes were fewer in flower. This indicated that MYB genes in different subgroups may perform different functions during the growth and development of petunia.

Expression of PaMYB genes in Anthocyanin Biosynthesis
Based on the promoter analysis, the expression of seven genes that contained MBSI ciselement, which was involved in flavonoid biosynthetic regulation, were quantified during petal coloring process of white and purple colored petunia ( Figure 6). Two genes were not represented probably because that they were pseudogenes or not examined in these varieties. PaMYB100 (annotated as Deep purple, DPL) showed relative higher expression level than other genes, and positive related to the coloration. This is consistent to its previously report [34]. Although Pa1RMYB36 had a high expression, the differences in the coloring process and in different color plants were not significant. Except for flavonoid biosynthetic regulation element, drought-inducibility and SA-responsive element existed, indicating that Pa1RMYB36 should have a prominent function in some coloring process related to abiotic stress response. The mRNA level of PaMYB82 was only detected in bud. PaMYB68 was not expressed in blooming. This result indicated that their function in flavonoid biosynthesis was regulated by other environment or development signals. PaMYB genes in in petal, leaf and roots. The gene expression level of tissue with highest expression level was considered as 1. Blocks with colors indicate decreased (blue) or increased (red) transcript accumulation relative to the respective control.

Expression of PaMYB genes in Anthocyanin Biosynthesis
Based on the promoter analysis, the expression of seven genes that contained MBSI cis-element, which was involved in flavonoid biosynthetic regulation, were quantified during petal coloring process of white and purple colored petunia ( Figure 6). Two genes were not represented probably because that they were pseudogenes or not examined in these varieties. PaMYB100 (annotated as Deep purple, DPL) showed relative higher expression level than other genes, and positive related to the coloration. This is consistent to its previously report [34]. Although Pa1RMYB36 had a high expression, the differences in the coloring process and in different color plants were not significant. Except for flavonoid biosynthetic regulation element, drought-inducibility and SA-responsive element existed, indicating that Pa1RMYB36 should have a prominent function in some coloring process related to abiotic stress response. The mRNA level of PaMYB82 was only detected in bud. PaMYB68 was not expressed in blooming. This result indicated that their function in flavonoid biosynthesis was regulated by other environment or development signals. Because genes in S3, P25 (S4), S5, P11 (S6), P9 (S7) and P16 (S12) were related to primary and secondary metabolic processes in Arabidopsis, we examined the expression level of these genes in pigmentation and different colored petals (Figure 7). Four colors of petals were selected, e.g., white, pink, purple and deep purple. To determine the gene expression level in different pigmentation, the expression level of R2R3-MYBs was detected in different colored petals. Because they have different total anthocyanin contents and anthocyanin components. On the other hand, to determine the gene expression level in anthocyanin biosynthesis, the expression pattern of R2R3-MYB genes was analyzed by comparing the three development stages of white and purple flowers. The expression pattern of genes in the same subgroup was similar in pigmentation and different colored petals. PaMYB33/AN2 was belonged to P11 (S6), which was involved in flavonoid biosynthesis [41]. The expression level of PaMYB33/AN2 was elevated during pigmentation. Genes in S9 (P7) were related to the flavonoid synthesis in vegetative organs, such as AtMYB111 and AtMYB12 [42,43]. PaMYB1 and PaMYB40 were mainly expressed in the budding period and colorless petal, suggesting that they were involved in regulation of flavonoid Because genes in S3, P25 (S4), S5, P11 (S6), P9 (S7) and P16 (S12) were related to primary and secondary metabolic processes in Arabidopsis, we examined the expression level of these genes in pigmentation and different colored petals (Figure 7). Four colors of petals were selected, e.g., white, pink, purple and deep purple. To determine the gene expression level in different pigmentation, the expression level of R2R3-MYBs was detected in different colored petals. Because they have different total anthocyanin contents and anthocyanin components. On the other hand, to determine the gene expression level in anthocyanin biosynthesis, the expression pattern of R2R3-MYB genes was analyzed by comparing the three development stages of white and purple flowers. The expression pattern of genes in the same subgroup was similar in pigmentation and different colored petals. PaMYB33/AN2 was belonged to P11 (S6), which was involved in flavonoid biosynthesis [41]. The expression level of PaMYB33/AN2 was elevated during pigmentation. Genes in S9 (P7) were related to the flavonoid synthesis in vegetative organs, such as AtMYB111 and AtMYB12 [42,43]. PaMYB1 and PaMYB40 were mainly expressed in the budding period and colorless petal, suggesting that they were involved in regulation of flavonoid synthesis in petunia vegetative tissues. The expression pattern of genes in P16 (S12) was not significantly related to the pigmentation. We speculated that they may act as regulators in other secondary metabolic processes. synthesis in petunia vegetative tissues. The expression pattern of genes in P16 (S12) was not significantly related to the pigmentation. We speculated that they may act as regulators in other secondary metabolic processes.

Discussion
The MYB transcription factor (TF) family was one of the largest TF families in plants. Although many PaMYBs have been reported, it is necessary to build an integral MYB family in Petunia for promoting phylogenetic and functional cognition. On the other hand, knowledge of the functions of certain members should facilitate the confirmation of paralogous and orthologous functional relationships.

Discussion
The MYB transcription factor (TF) family was one of the largest TF families in plants. Although many PaMYBs have been reported, it is necessary to build an integral MYB family in Petunia for promoting phylogenetic and functional cognition. On the other hand, knowledge of the functions of certain members should facilitate the confirmation of paralogous and orthologous functional relationships.
In this study, the MYB members in Petunia axillaris genome were identified, including 106 R2R3-MYBs, 7 R1R2R3-MYBs, 2 4R-MYBs and 40 1R-MYB genes. R2R3-MYB subfamily is the most abundant subfamily in petunia, which is consistent with previous reports, with 126 R2R3-MYB members in A. thaliana [12], 121 in tomato [39]. MYB genes have undergone large-scale expansion in higher plants [44]. Restricted by the draft genome of petunia, we cannot localize the PaMYB genes to chromosomes. However, based on the number of PaMYBs, the MYB genes should be experienced gene expansion related to the genome duplications. For example, PaMYB/ASR was a duplication of the genomic fragment containing the other three R2R3-MYB genes with roles in pigmentation that were previously defined, the ANTHOCYANIN4-DEEP PURPLE-PURPLE HAZE (AN4-DPL-PHZ) cluster, through comparative, functional and phylogenic analysis of anthocyanin R2R3-MYB genes [35].
The number of intron and exon in the same group (subgroup) was similar. Interestingly, the most complex structure was 3R-MYB, not 4R-MYB. This result indicated that the number of conserved domains was not related to the complexity of gene structure. The motif composition and distribution orders in the same group (subgroup) were also similar, while the structures among the different groups significantly diverged. For example, only group 24 contained motif 10. The motif composition in group 16 was significantly different with other groups. Conserved motifs within the same TF family may play crucial roles in protein-specific functions. Three Trp (W) in R2-repeat and the second and third W in R3-repeat were conserved in petunia, as reported in plant MYB binding domain [45]. Like other plant R3-repeat, the first Trp (W) was usually replaced by other amino acid residue. W was replaced by F (69%)\I (13%) in most PaR2R3-MYBs, e.g., F/I in potato [9], F in soybean [13]. Except for the sequence conservation of R2-and R3-repeats, the sequence in the linkage region of PaMYBs were also highly conserved, which was consistent to tomato, explaining by the closely evolutionary relationship within Solanaceae [39].
As the most type of MYB genes, a phylogenetic tree of R2R3-MYB protein of Petunia and Arabidopsis was constructed for further explaining the evolutionary relationship of R2R3-MYB in petunia. R2R3-MYB in Petunia was divided into 28 subgroups referred to Arabidopsis [40], while seven of them were not clustered with Arabidopsis. In general, the orthologs clustered in a subgroup share similar gene structure and functions, indicating common evolutionary origins. AtMYB4 (subgroup 4) can repress the synthesis of sinapoylmalate [46]. In subgroup 25 of petunia, the clustered group of S4, PaMYB101/PhMYB4 functioned in the repression of C4H transcription, indirectly regulating the floral volatile signature of petunia [47]. This reflected similar functions in the same subgroup, also it revealed that PaMYB23, PaMYB53, PaMYB62 and PaMYB65 may also be involved in primary and secondary metabolism. In P21 (S19), PaMYB36/EOBI and PaMYB28/EOBII were the homologous genes of AtMYB21 and AtMYB24, respectively [48][49][50]. In Arabidopsis, MYB99 (belongs to S17) acts in a regulatory triad with MYB21 and MYB24. PaMYB3/ODORANT1 (belongs to P15), a putative ortholog of MYB99, could act with MYB21 and MYB24 to co-regulate petunia scent biosynthesis genes [51]. In S7 (P9), AtMYB12, AtMYB11 and At-MYB111 were significantly up-regulated by UV-B light to produce flavonols to protect from light damage [52]. PaMYB91/MYB-FL was the major determinant of differences in flavonol levels and UV absorption. The gain of UV absorbance in the derivation of P. axillaris from P. inflata was caused by upregulation of the MYB-FL promoter, whereas a frameshift mutation in PaMYB91/MYB-FL was responsible for the difference in UV absorbance between P. axillaris and P. exserta [53]. This result suggested the similar function within S7 (P9). In addition, genes in S7 (P9) were related to the flavonoid synthesis in vegetative tissues (stems, leaves and sepals) in response to stressors such as high light [42]. PaMYB1 and PaMYB40 were mainly expressed in the bud stage and colorless petal, suggesting that they were involved in regulation of flavonoid synthesis in petunia vegetative organs.
In Arabidopsis, the anthocyanin regulator MYB75/PAP1 and MYB90/PAP2 belonged to subgroup 6 (S6) [18,41], which was clustered with subgroup 11 (P11) in petunia, with five PaR2R3-MYB genes included. PaMYB33/AN2 played a central role in petal limb pigmentation, while PaMYB77/PaMYB78/AN4 was the regulator of anthocyanin in petal tube and anthers [32,33]. The mRNA level of PaMYB33/AN2 and PaMYB77/PaMYB78/AN4 were elevated in the process of pigmentation (Figure 7). PaMYB100/DPL regulated veinassociated anthocyanin pigmentation in the flower tube, while PaMYB81/PHZ determined light-induced anthocyanin accumulation on exposed petal surfaces (bud-blush) [34]. Similarity, the expression level in purple petal of PaMYB100/DPL and PaMYB81/PHZ were higher than that in white petal. Among these five genes, only promoter of PaMYB100/DPL contained flavonoid biosynthetic regulation (MBSI) element, suggesting that MBSI was not necessary for regulation of anthocyanin biosynthesis. Genes in subgroup 9 were related to the flavonoid synthesis in vegetative organs, such as AtMYB111 and AtMYB12 [42,43]. In petunia subgroup 7, PaMYB1 and PaMYB40 were mainly expressed in the bud stage and colorless petal, suggesting that they were involved in regulation of flavonoid synthesis in petunia vegetative organs. In addition, ANTHOCYANIN SYNTHESIS REGULATOR (ASR) was the novel anthocyanin regulators in the genome of the purple flowering P. inflata S6 wild accession, while it was absent in the white flowering species P. axillaris [35]. Apart from positive regulators of anthocyanin, negative regulators were also reported. PaMYB105/PhMYB27, a putative R2R3-MYB anthocyanin repressor that functions as part of the MYB-bHLH-WD40 (MBW) complex and represses transcription through its C-terminal EAR motif, was significantly up-regulated in white petal. Another anthocyanin repressor, PhMYBx (R3-MYB) was a competitive repressor of AN1, which was an essential component of the MBW activation complex for pigmentation [36]. However, MYBx was not identified in the white flowering P. axillaris in this work.
Based on cis-element analysis, seven gene promoters contained flavonoid biosynthetic regulation element, indicating their regulatory roles in flavonoid biosynthesis. Except for PaMYB100/DPL, another six genes were not reported as the regulator of flavonoid biosynthesis. In the most cases, R2R3-MYBs were the main regulators of flavonoid biosynthesis. In our study, a MYB-related gene, Pa1RMYB36 contained a flavonoid biosynthetic regulation element, suggesting that it was a potential regulator of flavonoid biosynthesis. Although its expression has little relationship to the coloration, it may play an important role in another flavonoid biosynthetic pathway. The regulator of grape berry skin, VvmybA1, was a MYBrelated gene. In addition, its homologs, VvMYBA2r, VlmybA1-1, VlmybA1-2, and VlmybA2, also regulated anthocyanin biosynthesis [54]. Through the identification, classification and evolutionary analysis of MYB gene family in Petunia, our result was supported by published data. Except for existing knowledge, many new MYB genes were explored and were of research value. Nevertheless, due to the limitation of the reference genome, there are still some MYB genes that have not been identified, so follow-up research will expand the genome information in order to obtain more comprehensive PaMYB family data.

Identification of Petunia MYB Proteins
To identify potential members of petunia MYB gene family, we performed multiple database searches. The amino acids sequences of the Arabidopsis MYB proteins were downloaded from the Arabidopsis Information Resource (TAIR, http://www.Arabidopsis. org/, accessed on 1 June 2020), and were employed as a query in BLAST searches against the Petunia axillaris genome database (https://www.sgn.cornell.edu/, accessed on 1 June 2020) [37]. All tentative consensus sequences and singlet sequences with an e-value of <1.0 were selected after the homology search. Any uncompleted sequences or somewhat different conserved domains were discarded from the data list. Then, the cDNA sequences, amino acid sequences and genomic sequences were extracted for further analysis. The theoretical isoelectric point (PI) and molecular weights (MW) of the PaMYB proteins were predicted by Compute PI/MW tool on the ExPASy server (http://web.expasy.org/ compute_pi/, accessed on 1 June 2020).

Phylogenetic Analysis
The alignments of full-length amino acid sequences of PaMYBs and AtMYBs were performed using ClustalX 2.1 with default settings [55]. An unrooted phylogenetic tree was constructed based on the alignments using MEGA 10.1.8 with the Neighbor-Joining (NJ) method [56]. The parameters used in the tree construction were JTT model plus gamma-distributed rates determined by ProtTest 3.0 [57] and 1000 bootstraps. The pairwise gap deletion mode was used to ensure that the more divergent C-terminal domains could contribute to the topology of the NJ tree. We further performed a NJ phylogenetic tree of the R2R3-MYB protein sequences of Petunia and Arabidopsis for classifying PaR2R3-MYB.

Conserved Motif Identification and Gene Structure Analysis
The deduced amino acid sequences of the PaMYBs were analyzed by Multiple Em for Motif Elicitation (MEME; http://memesuite.org/tools/meme, accessed on 1 July 2020) for motif analysis [58]. To identify conserved motifs in these sequences, selection of the maximum number of motifs was set to 10 and optimum width of motif from 6 to 200. The motif with an e-value less than 1e −10 was retained for further analysis. The amino acid sequence in PaMYB motifs were manually adjusted using Weblogo online program (http://weblog.berkeley.edu/logo.cgi, accessed on 1 July 2020) to display the characteristics of the MYB domains [59]. The exon-intron structure of these genes was graphically displayed by the Gene Structure Display Server [60] based on the genomic sequences and the corresponding coding sequences (CDS).

Promoter Cis-Element Analysis
The promoter sequences (2000 bp upstream of start codon) of PaMYB genes were collected. The cis-acting elements are predicted in PlantCARE [61]. Statistical analysis and data visualization were done by Excel and TBtools [62].

Plant Materials and Growth Condition
Four varieties of Petunia × hybrida inbred lines with different petal colors were used in this research. These inbred lines showed different petal colors. The plants were grown in a chamber with a 16-h light/8-h dark cycle under 25 • C. The petal, leaf and root of purple flowering variety were used as samples for tissue specific expression detection. To determine the gene expression level in anthocyanin biosynthesis, the petals were taken in different developmental stages of white and purple flowering petunia, including budding period (S1), expansion period (S2) and blooming period (S3). To determine the gene expression level in different pigmentation, four colored petals in S3 were employed. Three biological replicates were conducted for each group.

Quantitative RT-PCR Analysis
RNA of all the samples was extracted using plant RNA extraction kit (TaKaRa Bio Inc., Dalian, China) and reversed to single strand cDNA using Prime-Script RT reagent kit (TaKaRa Bio Inc., Dalian, China) according to the instructions. qRT-PCR was performed by SYBR Premix EX Taq II (TaKaRa Bio Inc., Dalian, China) on the Eppendorf Mastercycler ep Realplex system. Primers are listed in Supplementary Table S6. The gene expression quantification was calculated using the 2 −∆∆Ct method.

Conclusions
In this work, 155 MYB genes in P. axillaris were identified, including 106 R2R3-MYB, 7 3R-MYB, 2 4R-MYB and 40 MYB-related genes. The gene structure and motif composition were highly conserved in petunia MYB family. Through analysis of evolutionary relationship of R2R3-MYB genes between Arabidopsis and Petunia, genes in the same subgroup shared similar functions, which was supported by published data. We also obtained several unknown MYB genes had potential function in anthocyanin biosynthesis by means of promoter analysis. This work provides a comprehensive understating of the MYB gene family in Petunia and lays a significant foundation for future studies on the function and evolution of MYB genes in petunia.