Genome-Wide Characterization of R2R3-MYB Transcription Factors in Pitaya Reveals a R2R3-MYB Repressor HuMYB1 Involved in Fruit Ripening through Regulation of Betalain Biosynthesis by Repressing Betalain Biosynthesis-Related Genes

The MYB (myeloblastosis) superfamily constitutes one of the most abundant transcription factors (TFs) regulating various biological processes in plants. However, the molecular characteristics and functions of MYB TFs in pitaya remain unclear. To date, no genome-wide characterization analysis of this gene family has been conducted in the Cactaceae species. In this study, 105 R2R3-MYB members were identified from the genome data of Hylocereus undatus and their conserved motifs, physiological and biochemical characteristics, chromosome locations, synteny relationship, gene structure and phylogeny were further analyzed. Expression analyses suggested that three up-regulated HuMYBs and twenty-two down-regulated HuMYBs were probably involved in fruit ripening of pitaya. Phylogenetic analyses of R2R3-MYB repressors showed that seven HuMYBs (HuMYB1, HuMYB21, HuMYB48, HuMYB49, HuMYB72, HuMYB78 and HuMYB101) were in clades containing R2R3-MYB repressors. HuMYB1 and HuMYB21 were significantly down-regulated with the betalain accumulation during fruit ripening of ‘Guanhuahong’ pitaya (H. monacanthus). However, only HuMYB1 had R2 and R3 repeats with C1, C2, C3 and C4 motifs. HuMYB1 was localized exclusively to the nucleus and exhibited transcriptional inhibition capacities. Dual luciferase reporter assay demonstrated that HuMYB1 inhibited the expression of betalain-related genes: HuADH1, HuCYP76AD1-1 and HuDODA1. These results suggested that HuMYB1 is a potential repressor of betalain biosynthesis during pitaya fruit ripening. Our results provide the first genome-wide analyses of the R2R3-MYB subfamily involved in pitaya betalain biosynthesis and will facilitate functional analysis of this gene family in the future.


Introduction
Transcription factors (TFs) play important roles in plant growth and development by controlling gene expression, suppressing translation or activating transcription of target genes, or interacting with other proteins. Generally, TFs are composed of at least four discrete domains: DNA-binding domain, nuclear localization signal, transcriptionactivation domain and oligomerization site [1]. The MYB family is one of the largest TF families known in eukaryotes and is defined by the presence of one to four highly conserved MYB repeats. The MYB repeat contains approximately 52 amino acids with three regularly spaced tryptophan, encoding three α-helices. The third helix of each available genomic information. The aim of the present study is to identify and validate the R2R3-MYB involved in pitaya betalain biosynthesis and candidate R2R3-MYB genes can be used for genetic improvement of pitaya.

Identification and Sequence Analyses of Pitaya R2R3-MYB Subfamily
Plot hidden Markov model (HMM) profile of MYB DNA-binding domain (PF00249) was downloaded from Pfam (https://pfam.xfam.org/, accessed on 8 September 2019) [26] to identify MYB genes from pitaya genome [27]. To ensure the presence of the core MYB domains, the putative MYB sequences were further examined using Pfam and SMART (https://smart.embl-heidelberg.de/) (accessed on 27 March 2020). The sequences with two MYB domains were considered to be the final R2R3-MYB gene family members.

Chromosomal Locations and Synteny Analyses of Pitaya R2R3-MYB Members
The MapChart software was used to draw the location sites of HuMYBs according to their physical positions. The duplication pattern for each HuMYBs was analyzed by Multiple Collonearity Scan toolkit (MCscanX) software following the operation manual [29]. The synteny relationship of orthologous R2R3-MYB genes between pitaya, beet and Arabidopsis was shown by TBtools software [28].

Plant Materials
'Guanhuahong' pitaya (red peel and green scales with red pulp, H. monacanthus) from the orchard of Jinsuinong (Zhongluotan Village, Guangzhou City, Guangdong Province, China) were used as materials. Fruits of 'Guanhuahong' pitaya were collected on the 14th, 17th, 19th, 23rd, 25th, 27th and 32nd day after flowering (DAF). All samples were frozen in liquid nitrogen immediately and stored at −80°C until future analyses.

Measurement of Betalain Contents
Betalains were extracted and measured following our previously described method [30]. In brief, 0.5 g freeze-dried pitaya pulp powder were extracted with 5 mL 80% aqueous methanol (v/v) solution by sonication for 10 min and stirred for 20 min in dark at room temperature. After centrifuging at 5000× g for 15 min, the residues were subjected to a similar second extraction. The supernatants were measured through spectrophotometry (Infinite M200, Tecan Co., Ltd, Shanghai, China) at 478 nm for betaxanthins and 538 nm for betacyanins. All determinations were performed in three biological repetition with three technical replicates.

Gene Expression Analyses and Cloning
Transcriptome data of four fruit developmental stages (17th, 23rd, 25th and 32nd) (PRJNA704510) were used to draw the transcript abundance of 2R-MYB repressors by TBtools software. The primers were designed in NCBI (https://www.ncbi.nlm.nih.gov/ tools/primer-blast/, accessed on 19 May 2020). Total RNA was isolated using the EASYspin Plus Complex Plant RNA Kit (RN53) (Aidlab Biotechnology, Beijing, China) according to the manufacturer's protocol. Single-stranded cDNA was synthesized using the PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan). qRT-PCR was performed with an CFX384-Real-Time System (C1000 Touch Thermal Cycler, Bio-Rad, Irvine, CA, USA) using the RealUniversal Color PreMix (SYBR Green) (TIANGEN, Beijing, China) with specific primers (Table S1). The Actin (1) was used as internal control for gene expression analyses [31]. All experiments were repeated in technical replicates. The relative expression levels of each transcript were calculated using the comparative Ct method [32].
The full-length coding sequences of HuMYBs were cloned using I-5 TM 2× High-Fidelity Master Mix (MCLAB, San Francisco, CA, USA) with specific primers (Table S1). Alignment of the deduced protein sequences was performed by Clustal X and GeneDoc software.

Transcriptional Activation Analyses in Yeast Cells
Full-length of HuMYB1 were inserted into the pGBKT7 vector (primers are listed in Table S1). Then, the pGBKT7-HuMYB1, positive (pGBKT7-53 + pGADT7-T-antigen) and negative (pGBKT7) controls were separately transformed into yeast strain GoldY2H (Coolaber, Beijing, China) and grown on minimal medium plates without tryptophan (SD/-Trp) or without tryptophan, histidine and adenine (SD/-Trp-His-Ade). The transactivation activity of HuMYB1 protein was evaluated according to their growth status after 3 days in 30°C and confirmed by incubating with X-α-galactosidase (X-α-Gal) for 30 min.
The dual-luciferase transient expression system was conducted according to the previous study [34]. To assess the transactivation activity of HuMYB1, full-length coding sequences of HuMYB1 was inserted into the pGreenII BD-62-SK as effectors. The reporter vector was modified from the pGreenII 0800-LUC vector [33]. Then, effectors of pGreenII BD-62-SK-HuMYB1, pGreenII BD-62-SK, pGreenII BD-62-SK-VP16 (positive control) and reporter were separately transformed into Agrobacterium tumefaciens strain GV3101. The effector and reporter were infiltrated into N. benthamiana leaves with a mixture of 9:1, and LUC and REN activities were measured by dual-luciferase assay kit (Promega, Madison, WI, USA) after 72 h.
To assess the transcriptional effects of HuMYB1 on the promoters of ADH1, CYP76AD1-1 and DODA1, these promoters were inserted into the pGreenII 0800-LUC vector as reporters, while HuMYB1 were inserted into the pGreenII 62-SK vector as effectors. The reporter and effector were transformed into Agrobacterium tumefaciens strain GV3101 and then infiltrated into N. benthamiana leaves. LUC and REN activities were analyzed as above.

Identification of R2R3-MYB Genes from H. undatus
The pitaya MYB encoding genes were identified by the alignment of the pitaya genome and MYB domain HMM profile (PF00249) using BLASTP software. A total of 231 candidate deduced amino acid sequences containing MYB or MYB-like repeats were obtained. The MYB domains were subsequently analyzed by Pfam and SMART. In total, 75 1R-MYB proteins, 105 R2R3-MYB proteins, four R1R2R3-MYB proteins and one 4R-MYB protein were obtained from the H. undatus genome (Table S2). Besides this, the theoretical pI and molecular weight ranged from 4.26 kDa (HuMYB137) to 10.96 kDa (HuMYB171) and 8.91 kDa (HuMYB143) to 213.32 kDa (HuMYB154), respectively (Table S2).
To investigate the homologous domain sequence features, conservation and divergence of 105, 70 and 126 R2 and R3 repeats from H. undatus (betalain-producing plant), B. vulgaris (betalain-producing plant) and A. thaliana (anthocyanin-producing plant) were assessed using multiple alignment analyses ( Figure 1). The results showed that the basic regions of pitaya, beet and Arabidopsis MYB domains (about 104 amino acid residues) were rarely deleted or inserted, as previously reported in the other plants [35,36], while they were widely divergent in the length and amino acid outside of the MYB domains. However, the length of Chinese cabbage and rice were larger than that of Arabidopsis based on the space between neighboring tryptophan (Trp, W) residues [37]. In pitaya, six highly conserved triplet W residues were located in R2 repeat (position 6, 26 and 46), R3 repeat (position 78 and 97) ( Figure 1A,B) and the first Trp residue in R3 repeat (position 59) was replaced by a hydrophobic amino acid Phe (F). These conserved W (F) residues are indispensable in maintaining the helix-turn-helix structure of MYB domains [38]. Similar localization existed in both beet and Arabidopsis, suggesting the evolutionary conservation of MYBs among plants ( Figure 1C,F). Besides this, the 3' region of R2 domain in pitaya R2R3-MYBs contained a highly conserved LRPD motif. Compared to Arabidopsis, the positions at 12, 19, 25, 33 and 36 in R2 repeats and 69, 80 and 83 in R3 repeats showed difference in pitaya ( Figure 1). Difference was detected between pitaya and beet in the position 69 and 80 in R3 repeats, indicating the divergence of the MYB domain between pitaya and beet was less than that of pitaya and Arabidopsis.

Genomic Location of HuMYBs
Genome chromosomal location analyses showed that all the 105 R2R3-MYB genes were distributed throughout all 11 chromosomes (Chrs) in the pitaya genome ( Figure 2) and renamed as HuMYB1 to HuMYB105 according to their location in 11 Chrs. The largest number of MYB genes (fifteen) mapped on Chr 6, followed by thirteen on Chr 5. However, there were only three HuMYBs on Chr 9. In general, MYB genes are absent in the central sections of chromosomes including centromere and pericentromeric regions. Relatively high densities of HuMYB genes were observed at the end of Chr 3, 6 and 7 as well as the bottom of Chr 1, 4, 5, 6, 8 and 11 ( Figure 2). The same uneven distributions in all chromosomes were also reported in the other plants, such as A. comosus [13], Brassica rapa [37] and Eucalyptus grandis [39].

Genomic Location of HuMYBs
Genome chromosomal location analyses showed that all the 105 R2R3-MYB genes were distributed throughout all 11 chromosomes (Chrs) in the pitaya genome ( Figure 2) and renamed as HuMYB1 to HuMYB105 according to their location in 11 Chrs. The largest number of MYB genes (fifteen) mapped on Chr 6, followed by thirteen on Chr 5. However, there were only three HuMYBs on Chr 9. In general, MYB genes are absent in the central sections of chromosomes including centromere and pericentromeric regions. Relatively high densities of HuMYB genes were observed at the end of Chr 3, 6 and 7 as well as the bottom of Chr 1, 4, 5, 6, 8 and 11 ( Figure 2). The same uneven distributions in all chromosomes were also reported in the other plants, such as A. comosus [13], Brassica rapa [37] and Eucalyptus grandis [39].

Synteny Analyses of Pitaya R2R3-MYB Subfamily
Gene duplication is recognized to occur throughout plant evolution and plays an important role in expanding the large gene families in plants [40]. The duplication patterns and synteny relationship of the HuMYBs were analyzed using BLASTP and MCScanX methods. In total, 40 segmental duplication events with 61 HuMYBs happened in the pitaya genome ( Figure 3A; Table S3). HuMYBs were located within synteny blocks on almost all chromosomes except for Chr 9, indicating that HuMYBs underwent through expansion during genome evolution. If two or more MYB genes resided within 20 kb in the same chromosome, a tandem duplication event was defined [41]. In this study, three HuMYB tandem duplication pairs happened in Chr 5 (HuMYB48 and HuMYB49), Chr 6 (HuMYB67 and HuMYB68) and Chr 10 (HuMYB93 and HuMYB94). These results suggested that large-scale segmental duplication events were the major contributors to the expansion of the R2R3-MYB gene family in H. undatus.
To further investigate the potential evolutionary mechanism of R2R3-MYB subfamily, we constructed the comparative syntenic maps of H. undatus, B. vulgaris and A. thaliana in pairs ( Figure 3B). Finally, 46 collinear R2R3-MYB gene pairs between pitaya and beet ( Figure 3B(a); Table S4-1), 28 between pitaya and Arabidopsis ( Figure 3B(b); Table S4-2) and 10 between beet and Arabidopsis ( Figure 3B(c); Table S4-3) were identified. The number of orthologous events of HuMYB-BvMYB was much more than that of HuMYB-AtMYB and BvMYB-AtMYB, indicating that the evolutionary distance between pitaya and beet was closer than that of pitaya and Arabidopsis. In addition, the area of all synteny blocks between pitaya and beet was far greater than that of pitaya and Arabidopsis. These results suggested that the evolutionary distance between pitaya and beet is closer than that of pitaya and Arabidopsis.

Phylogenetic, Gene Structure and Motif Composition Analyses of R2R3-MYBs
To explore the putative function of the 105 HuMYBs, functional groups were performed using MYB genes from B. vulgaris and A. thaliana. A. thaliana was chosen because MYB genes have been extensively studied in this model plant while B. vulgaris was a vital plant material for betalain studies. Most MYB proteins shared similar functions cluster in the same phylogenetic clades, suggesting that most closely-related MYBs could recognize similar target genes and possess redundant, overlapping and/or cooperative functions [42].
By the ML method, an evolutionary tree of 105 HuMYBs, 70 BvMYBs and 126 AtMYBs was constructed ( Figure 4). In the classification of the R2R3-MYB genes, the subgroup (S) categories from A. thaliana were also labeled in the phylogenetic tree [7]. According to the phylogenetic tree topology, 105 HuMYB proteins were classified into nine groups (designated as Group 1 to Group 9). HuMYBs were distributed with AtMYBs and BvMYBs in most subgroups. However, species-specific AtMYBs were also detected in S10, S12 and S15 (clades without HuMYBs and BvMYBs).

Synteny Analyses of Pitaya R2R3-MYB Subfamily
Gene duplication is recognized to occur throughout plant evolution and plays an important role in expanding the large gene families in plants [40]. The duplication patterns and synteny relationship of the HuMYBs were analyzed using BLASTP and MCScanX methods. In total, 40 segmental duplication events with 61 HuMYBs happened in the pitaya genome ( Figure 3A; Table S3). HuMYBs were located within synteny blocks on almost all chromosomes except for Chr 9, indicating that HuMYBs underwent through ex-  Table S4-3) were identified. The number of orthologous events of HuMYB-BvMYB was much more than that of HuMYB-AtMYB and BvMYB-AtMYB, indicating that the evolutionary distance between pitaya and beet was closer than that of pitaya and Arabidopsis. In addition, the area of all synteny blocks between pitaya and beet was far greater than that of pitaya and Arabidopsis. These results suggested that the evolutionary distance between pitaya and beet is closer than that of pitaya and Arabidopsis.

Phylogenetic, Gene Structure and Motif Composition Analyses of R2R3-MYBs
To explore the putative function of the 105 HuMYBs, functional groups were performed using MYB genes from B. vulgaris and A. thaliana. A. thaliana was chosen because MYB genes have been extensively studied in this model plant while B. vulgaris was a vital plant material for betalain studies. Most MYB proteins shared similar functions cluster in the same phylogenetic clades, suggesting that most closely-related MYBs could recognize similar target genes and possess redundant, overlapping and/or cooperative functions [42]. By the ML method, an evolutionary tree of 105 HuMYBs, 70 BvMYBs and 126 AtMYBs was constructed ( Figure 4). In the classification of the R2R3-MYB genes, the subgroup (S) HpWRKY44 and BvMYB1 are involved in betalain biosynthesis but they share high homology with anthocyanin regulated TFs including AtWRKY44, AtMYB75, AtMYB90, AtMYB113 and AtMYB114 [17,43]. Thus, there is probably a link of transcription regulation between betalains and anthocyanins. Ten R2R3-MYB genes in Group 1 attracted our attention due to their gathering in S4 to S7 from Arabidopsis and Bv_ihfg ( Figure 4). Particularly, HuMYB12 was clustered with Bv_iogq and AtMYB111 in S7, which controlled flavonol biosynthesis [3,16]. HuMYB83 was close to Bv_ralf, Bv_jkkr and AtMYB90 in the S6 group, which is involved in betalain and anthocyanin biosynthesis [3,16,17]. HuMYB26 was near Bv_crae and AtMYB123 in S5, responsible for proanthocyanidins (PAs) biosynthesis in the seed coat of Arabidopsis [3]. In addition, five genes, i.e., HuMYB1, HuMYB21, HuMYB72, HuMYB78 and HuMYB101, were gathered in S4 which encoded transcriptional repressors [3]. HuMYB48 and HuMYB49 were grouped with Bv_ihfg containing negative flavonoid regulator FaMYB1 [16]. These results suggested that HuMYB12, HuMYB26 and HuMYB83 were candidate R2R3-MYB activators while HuMYB1, HuMYB21, HuMYB48, HuMYB49, HuMYB72, HuMYB78 and HuMYB101 were candidate transcriptional repressors.
Phylogenetic reconstruction of all 105 HuMYBs proteins were performed by the NJ method with 1000 bootstraps in MEGA7.0 ( Figure 5A). The 105 HuMYBs were divided into nine groups, which was consistent with the result of the phylogenetic tree constructed by ML method in Figure 4. To understand the distribution of conserved motifs, the 105 R2R3-MYB protein sequences were analyzed by MEME software with five conserved motifs (named as motif 1 to 5) ( Figure 5B and Figure S1). A total of 105 HuMYBs contained motifs 1, 2 and 3 which constituted R2 and R3 domains in the N-terminal, except for HuMYB42 in Group 6 and HuMYB9, HuMYB14, HuMYB38, HuMYB43, HuMYB75 and HuMYB89 in Group 9, which were found in the central region. Besides this, motif 4 and motif 5 were respectively conserved in Group 4 and Group 5, which is consistent with the phylogeny analyses of 105 HuMYBs proteins. categories from A. thaliana were also labeled in the phylogenetic tree [7]. According to the phylogenetic tree topology, 105 HuMYB proteins were classified into nine groups (designated as Group 1 to Group 9). HuMYBs were distributed with AtMYBs and BvMYBs in most subgroups. However, species-specific AtMYBs were also detected in S10, S12 and S15 (clades without HuMYBs and BvMYBs). HpWRKY44 and BvMYB1 are involved in betalain biosynthesis but they shared high homology with anthocyanin regulated TFs including AtWRKY44, AtMYB75, AtMYB90, AtMYB113 and AtMYB114 [17,43]. Thus, there is probably a link of transcription regulation between betalains and anthocyanins. Ten R2R3-MYB genes in Group 1 attracted our The exon-intron structures of 105 R2R3-MYB coding sequences were analyzed ( Figure 5C). Most R2R3-MYBs were clustered in the same group with the similar exon-intron structures, particularly with the same number of introns, such as Group 2 and Group 7. Except for six genes (HuMYB2, HuMYB4, HuMYB20, HuMYB46, HuMYB52 and HuMYB55) in Group 9 without any intron, most of HuMYBs were disrupted by introns. In addition, more than 73% HuMYBs had three exons and two introns, which was also found in other plants [13,14].
pressed on the 15 th DAF when the seeds were white ( Figure S3). Five genes were highly expressed on the 17 th DAF when the seeds became red ( Figure S3). Fourteen genes were predominantly expressed on the 19 th DAF when seeds were black ( Figure S3). Eleven genes were highly expressed on the 23 rd DAF when the betacyanin started to accumulate in pulps ( Figure S3). Twenty-two genes showed irregular trends during pitaya fruit ripening ( Figure S4). These results suggested that many HuMYBs may function as repressors during pitaya fruit ripening.

HuMYB1 Is a Nucleus-Localized Transcription Repressor
HuMYB1 probably plays an essential role in suppressing betalain biosynthesis in pitaya according to the above evolutionary, expression and sequence alignment analyses. To investigate the subcellular localizations of HuMYB1, the full-length coding sequence of HuMYB1 was fused with the GFP to construct 35S-HuMYB1-GFP vectors which were transiently transformed into N. benthamiana leaves. The nuclear signals were captured in protoplasts after overexpressing HuMYB1-GFP, whereas GFP alone showed cytoplasmatic and nuclear signals ( Figure 8A). vector was used as a calibrator (set as 1). BD-62SK-VP16 was used as a positive control. Asterisks represents highly significant differences at p value < 0.01 using two-tailed t-test, compared to BD-62SK. (E) Diagrams of the reporter and effector vectors. (F) HuMYB1 inhibited the transcription of betalain-related genes including HuADH1, HuCYP76AD1-1 and HuDODA1 by dual-luciferase transient expression assay in N. benthamiana leaves. The ratio of LUC/REN of the empty vector (62SK) plus promoter was used as a calibrator (set as 1). Asterisks represents highly significant differences at p value < 0.05 (one asterisks) and p value < 0.01 (two asterisks) using two-tailed t-test, compared to 62SK.  1). BD-62SK-VP16 was used as a positive control. Asterisks represents highly significant differences at p value < 0.01 using two-tailed t-test, compared to BD-62SK. (E) Diagrams of the reporter and effector vectors. (F) HuMYB1 inhibited the transcription of betalain-related genes including HuADH1, HuCYP76AD1-1 and HuDODA1 by dual-luciferase transient expression assay in N. benthamiana leaves. The ratio of LUC/REN of the empty vector (62SK) plus promoter was used as a calibrator (set as 1). Asterisks represents highly significant differences at p value < 0.05 (one asterisks) and p value < 0.01 (two asterisks) using two-tailed t-test, compared to 62SK.
Full-length coding regions of HuMYB1 were fused with the GAL4BD in the pGBKT7 vector to study the transcriptional activation abilities of HuMYB1 in yeast cells. As shown in Figure 8B, the transformed yeast cells of positive control (pGBKT753 + pGADT7-T) grew well in SD/-Trp-His-Ade and showed α-Gal activity, while yeast cells containing pGBKT7 (negative control) or pGBKT7-HuMYB1 did not, suggesting that HuMYB1 has no transactivation activities in yeast cells and probably functions as a transcriptional repressor in gene regulation.
The transcriptional repressor activities of HuMYB1 were confirmed in N. benthamiana leaves using the dual-luciferase reporter system ( Figure 8C). Compared to the expression of negative control (BD-62SK), the expression of the positive control (BD-62SK-VP16) resulted in a higher value of LUC/REN ratio while the BD-62SK-HuMYB1 showed a significantly lower value ( Figure 8D). It is consistent with the result from yeast cells. These results demonstrated that HuMYB1 was a nucleus-localized transcriptional repressor.

HuMYB1 Inhibited the Transcription of Three Betalain Biosynthesis-Related Genes
HuADH1, HuCYP76AD1-1 and HuDODA1 promoters were cloned from 'Guanhuahong' pitaya and the conserved cis-element motifs were predicted in Plant-CARE database. MYB binding sites were present in the promoters of HuADH1, HuCYP76AD1-1 and Hu-DODA1 (Table S6). Transient dual-luciferase assays in N. benthamiana leaves were performed to study whether HuMYB1 acted as the transcriptional repressor of HuADH1, HuCYP76AD1-1 and HuDODA1 ( Figure 8E). As shown in Figure 8F, compared to the empty control (pGreenII 62-SK), LUC/REN ratio was significantly reduced when HuMYB1 was co-transformed with HuADH1 or HuCYP76AD1-1 or HuDODA1 pro-LUC reporter. These results indicated that HuMYB1 participated in betalain biosynthesis by repressing the transcription of betalain biosynthesis-related HuADH1, HuCYP76AD1-1 and HuDODA1 genes ( Figure 8F).

Discussion
Although the genome-wide of MYB superfamily has been extensively studied in various plants, the identification and characterization of R2R3-MYB TFs based on the whole genome sequence of H. undatus were not reported yet. In our study, a total of 185 MYB genes (75 1R-MYB, 105 R2R3-MYB, four R1R2R3-MYB and one 4R-MYB) were identified by genome-wide search and distributed in all H. undatus Chrs (Figure 2). The R2R3-MYB proteins are the largest subfamily which is consistent with the findings in Arabidopsis [3], pineapple [13] and Salvia miltiorrhiza [36]. Compared with beet, the number of pitaya R2R3-MYB members (105) is fewer than Arabidopsis (126) [3] but more than beet (70) [16], indicating R2R3-MYB subfamilies in pitaya and Arabidopsis were expanded. Low tandem and high segmental duplications resulted in the expansion of the R2R3-MYB gene family in H. undatus, which also occurred in cotton [54] and pineapple [13]. Based on the comparison of DNA-binding domains (Figure 1), synteny ( Figure 3) and phylogenetic relationship (Figure 4) between pitaya, beet and Arabidopsis, closer genetic distance was detected between betalain-producing plants of pitaya and beet compared to pitaya and the anthocyanin-producing plant Arabidopsis [27].
Among 105 R2R3-MYB genes, expression patterns of 91 HuMYBs were analyzed in pulps of the seven fruit ripening stages of 'Guanhuahong' pitaya. Some genes showed preferential expressions on the 15th, 17th, 19th or 23rd DAF, suggesting that they play multiple regulatory roles during pitaya fruit ripening ( Figure S3). HuMYB3, HuMYB25 and HuMYB30 shared a similar expression pattern with betalain biosynthesis-related structure genes of HuADH1, HuCYP76AD1-1 and HuDODA1. However, 20.9% R2R3-MYB genes kept downward trends and were contrary to the expressions of HuADH1, HuCYP76AD1-1 and HuDODA1 during pitaya fruit ripening ( Figure 6). Additionally, HuMYB93 and HuMYB101 shared similar expression patterns with HucDOPA5GT2 involved in the formation of betanin ( Figure S3). ADH and ADT (arogenate dehydratase) can catalyze arogenate to form tyrosine and phenylalanine, respectively, and then form betalains and anthocyanins [55].
In the present study, HuMYB83 was clustered with BvMYB1 and AtMYB90. BvMYB1 is responsible for betalain biosynthesis [17], while AtMYB90 is involved in anthocyanin biosynthesis [3]. The transcription regulation of betalains and anthocyanins may exist a connection in their evolutionary relationship. Thus, HuMYB12, HuMYB26 and HuMYB83 may be involved in betalain biosynthesis since they were clustered with Arabidopsis subgroups 7, 5 and 6 involved in the regulation of flavonoid biosynthesis. Besides this, five AtMYB4-like repressors (HuMYB1, HuMYB21, HuMYB72, HuMYB78 and HuMYB101) and two FaMYB1-like repressors (HuMYB48 and HuMYB49) were identified according to the phylogenetic tree. In total, eight positive regulators and twenty-seven negative regulators were obtained by expression profiles and phylogenetic analyses (Figure 9). volved in betalain biosynthesis by targeting the CYP76AD1 and DODA promoter while HpWRKY44 has a functional role in the regulation of CYP76AD1 [17,43]. Thus, HuMYB TFs are probably involved in betalain biosynthesis by binding to the HuADH1, HuCYP76AD1-1 and HuDODA1 promoters in pitaya. Therefore, a hypothetical model of betalain biosynthesis in pitaya was proposed according to these candidate HuMYBs possibly involved in betacyanin biosynthesis by regulating the expression levels of HuADH1, HuCYP76AD1-1, HuDODA1 and/or GTs (Figure 9). Besides this, ABRE, E-box, G-box and W-box cis-elements were detected in the sequences of HuADH1, HuCYP76AD1-1 and Hu-DODA1 promoters, which were the binding sites of TFs such as AREB/ABF, bZIP, bHLH and WRKY [43,[56][57][58]. Therefore, these HuMYB TFs are possibly involved in betalain biosynthesis in pitaya. Moreover, the transcription activity of MYB proteins was influenced by the interaction with other TFs, such as bZIP, bHLH, NAC and WRKY TFs [59][60][61][62]. The activity of MYBs and bHLHs appears to be regulated by cytosolic WD40 repeat proteins by forming highly dynamic MBW complexes. WRKY can activate the transcription of anthocyanin-related genes by engaging with the MBW complex to form MBWW complex. Thus, further studies need to be carried out to determine whether the other TFs such as AREB/ABF, bZIP, bHLH and WRKY are involved in betalain biosynthesis, or HuMYB1 can coordinate with the other TFs to regulate betalain biosynthesis.  According to expression and sequence analyses, HuMYB1 showed a downward trend during fruit ripening ( Figure 7) and was closer to AtMYB4, which negatively regulated phenylpropanoid pathway [3]. Besides this, HuMYB1 acted as a nucleus-localized transcriptional repressor involved in betalain biosynthesis pathway by repressing the transcription of HuADH1, HuCYP76AD1-1 and HuDODA1 (Figure 8). BvMYB1 is involved in betalain biosynthesis by targeting the CYP76AD1 and DODA promoter while HpWRKY44 has a functional role in the regulation of CYP76AD1 [17,43]. Thus, HuMYB TFs are probably involved in betalain biosynthesis by binding to the HuADH1, HuCYP76AD1-1 and HuDODA1 promoters in pitaya. Therefore, a hypothetical model of betalain biosynthesis in pitaya was proposed according to these candidate HuMYBs possibly involved in betacyanin biosynthesis by regulating the expression levels of HuADH1, HuCYP76AD1-1, HuDODA1 and/or GTs (Figure 9). Besides this, ABRE, E-box, G-box and W-box cis-elements were detected in the sequences of HuADH1, HuCYP76AD1-1 and HuDODA1 promoters, which were the binding sites of TFs such as AREB/ABF, bZIP, bHLH and WRKY [43,[56][57][58]. Therefore, these HuMYB TFs are possibly involved in betalain biosynthesis in pitaya. Moreover, the transcription activity of MYB proteins was influenced by the interaction with other TFs, such as bZIP, bHLH, NAC and WRKY TFs [59][60][61][62]. The activity of MYBs and bHLHs appears to be regulated by cytosolic WD40 repeat proteins by forming highly dynamic MBW complexes. WRKY can activate the transcription of anthocyanin-related genes by engaging with the MBW complex to form MBWW complex. Thus, further studies need to be carried out to determine whether the other TFs such as AREB/ABF, bZIP, bHLH and WRKY are involved in betalain biosynthesis, or HuMYB1 can coordinate with the other TFs to regulate betalain biosynthesis.

Conclusions
In conclusion, this study is the first report on identification and characterization of R2R3-MYB gene family based on the genome of H. undatus. A total of 105 R2R3-MYB genes were obtained from the H. undatus genome. These R2R3-MYB genes were distributed in all 11 chromosomes with conserved R2 and R3 repeats and expanded its subfamily members through segmental duplication with a closer distance of beet than Arabidopsis. They were functionally classified into nine groups, consisting of conserved motif and exon-intron clustering. Phylogeny, sequence alignment and expression patterns of seven candidate R2R3-MYB repressors were further analyzed. HuMYB1 was grouped in the AtMYB1-like repressor group with bHLH binding motif and C1-C4 motifs. Besides this, HuMYB1 showed a downward expression trend opposite to betalain accumulation in pulps during fruit ripening of 'Guanhuahong' pitaya. HuMYB1 was a nucleus protein and could reduce the transcription activity of HuADH1, HuCYP76AD1-1 and HuDODA1 with transcription repression activities. The present study provides valuable information for a better understanding of MYB TFs involved in fruit ripening and betalain biosynthesis in pitaya.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are openly available in the NCBI PRJNA704510.