Genome-Wide Analysis of Specific PfR2R3-MYB Genes Related to Paulownia Witches’ Broom

Paulownia witches’ broom (PaWB), caused by phytoplasmas, is the most devastating infectious disease of Paulownia. R2R3-MYB transcription factors (TF) have been reported to be involved in the plant’s response to infections caused by these pathogens, but a comprehensive study of the R2R3-MYB genes in Paulownia has not been reported. In this study, we identified 138 R2R3-MYB genes distributed on 20 chromosomes of Paulownia fortunei. These genes were classified into 27 subfamilies based on their gene structures and phylogenetic relationships, which indicated that they have various evolutionary relationships and have undergone rich segmental replication events. We determined the expression patterns of the 138 R2R3-MYB genes of P. fortunei by analyzing the RNA sequencing data and found that PfR2R3-MYB15 was significantly up-regulated in P. fortunei in response to phytoplasma infections. PfR2R3-MYB15 was cloned and overexpressed in Populus trichocarpa. The results show that its overexpression induced branching symptoms. Subsequently, the subcellular localization results showed that PfR2R3-MYB15 was located in the nucleus. Yeast two-hybrid and bimolecular fluorescence complementation experiments showed that PfR2R3-MYB15 interacted with PfTAB2. The analysis of the PfR2R3-MYB15 gene showed that it not only played an important role in plant branching, but also might participate in the biosynthesis of photosystem elements. Our results will provide a foundation for future studies of the R2R3-MYB TF family in Paulownia and other plants.


Plant Materials
The tissue-cultured seedlings used in this study were from the Institute of Paulownia, Henan Agricultural University, Zhengzhou, Henan Province, China. Healthy P. fortunei seedlings (PF), phytoplasma-infected P. fortunei seedlings (PFI), and wild-type P. trichocarpa (WT) were cultured as described by Zhai et al. (2010) [22]. PF and PFI were collected at 30 days, immediately frozen in liquid nitrogen, and stored at −80 • C for further study.

Expression Pattern Analysis of PfR2R3-MYB Genes Responsive to Phytoplasma and qRT-PCR Verification
In previous studies, our research group found that the phytoplasma content decreased in PFI treated with 100 mg·L −1 of rifampin (Rif) with the extension in treatment time, and the plant morphologically showed reduction in branching and gradually reached a healthy state [28]. To further understand the role of PfR2R3-MYB genes in the occurrence of PaWB, we downloaded the RNA-seq data (https://www.ncbi.nlm.nih.gov/sra/?term=Paulownia, accessed on 27 June 2021), including PF and PFI, as well as PFI treated with 100 mg·L −1 Rif for 5, 10 and 20 days (PFI100-5, PFI100-10, and PFI100-20) (SRR11787883, SRR11787894, SRR11787905, and SRR11787912, and SRR11787921). Clean reads were aligned to the P. fortunei genome using BWA software. Gene expression levels were quantified by RSEM (RNASeq by expectation maximization) software package. The fragments per kb per million reads (FPKM) values were used to calculate the expression of genes. The heat map was generated using heatmap on TBtools. To validate the sequencing data, qRT-PCRs were performed to detect the expression patterns of eight differentially expressed genes (DEGs). We designed mRNA and universal reverse primer using Primer5 (Supplementary Materials: Table S1). The qRT-PCRs were performed on a CFX 96 Touch RT-PCR detection system (Bio-Rad, Hercules, CA, USA) with SYBR Green Real-time PCR Master Mix (Toyobo, Osaka, Japan) as described previously [29]. Expression levels were calculated by the 2 −∆∆CT method [30]. Actin (Forward primer: AATGGAATCTGCTGGAAT, Reverse primer: ACTGAGGACAATGTTACC) was used as the internal reference gene. Three biological replicates were used.

Vector Construction and Agrobacterium Tumefaciens Transformation
The PfR2R3-MYB15 gene primers were as follows: forward primer was 5 -ATGGGAA GAGCACCTTGCTG-3 and reverse primer was 5 -CAAATATCAACAGAAAGTTCTGAGT TTC-3 . The PfR2R3-MYB15 gene was driven by the cauliflower mosaic virus (CaMV) 35S promoter and was transformed into Agrobacterium tumefaciens (A. tumefaciens) by the freezethaw method [31]. WT was cultured on 1/2 woody plant medium (WPM) for 30 days, then stems were infected with A. tumefaciens GV3101 harboring the PfR2R3-MYB15 gene using cultures with an OD at 600 nm of 0.6; the details of this method referred to Song et al. (2006).

Detection of Transgenic Seedlings
PCR was used to identify the transgenic seedlings. The primer sequences designed by Primer Premier 5.0 were as follows: forward primer was 5 -GGATCCATGGGAAGAGCACC-3 and reverse primer was 5 -GAGCTCCTATTTGTACAATT-3 . The PCR cycling conditions were 95 • C for 3 min followed by 35 cycles of 95 • C for 30 s, 52 • C for 30 s, and 72 • C for 90 s, and melting at 72 • C for 10 min. QRT-PCR was used to validate OE-PfR2R3-MYB15. The method is the same as that in Section 2.6. Ribosomal PtActin RNA was used as the reference gene. The primer sequences were designed by Primer Premier 5.0 (Supplementary Materials: Table S1).

Data Availability Statement
All the raw read data will be provided during review.

Genome-Wide Identification of PfR2R3-MYB in P. fortunei
We identified 245 MYB genes in the P. fortunei genome; 101 were R1-MYBs, 6 were R1R2R3-MYBs and 138 were R2R3-MYBs (Supplementary Materials: Table S4). On the basis of their chromosome localization, we named the 138 PfR2R3-MYB genes from PfR2R3-MYB1 to PfR2R3-MYB138 (Supplementary Materials: Table S5). The length of  Table S5). The PfR2R3-MYB proteins were predicted to be located in the cell nucleus, which is consistent with the location of R2R3-MYB proteins in Arabidopsis [2] and poplar [8].

Localization, Gene Replication and Collinearity Analysis
The 138 PfR2R3-MYB genes were assigned to the 20 chromosomes of P. fortunei (Figure 1a), and their distribution on the chromosomes was scattered and widespread, which was consistent with the distribution of the PtR2R3-MYB genes on the poplar chromosomes [8]. Sixteen PfR2R3-MYBs were on chromosome 6, eleven each on chromosomes 11 and 13, six each on chromosomes 4, 8, 12, 14, 16, and 20, five each on chromosomes 10, 15, 17, and 19, four each on chromosomes 5 and 8, and three on chromosome 1.  Figure S1). The results showed that the PfR2R3-MYBs had common ancestors and were highly conserved during the evolutionary process.

Analysis of Chromosomal Location and Synteny of PfR2R3-MYB Genes
To clarify the classification of PfR2R3-MYB proteins, we constructed an NJ rootless Gene duplication is important for gene family formation and extension in plants [3]. In this study, we performed a collinearity analysis of the PfR2R3-MYB genes and identified 92 PfR2R3-MYB repeat pairs (Figure 1b). The PfR2R3-MYB gene duplication events were fragment duplications in all the P. fortunei chromosomes, implying that 92 fragment duplications played important roles in the evolution of the PfR2R3-MYB gene family (Supplementary Materials: Table S6). To assess the limiting factors of the evolution of the PfR2R3-MYB gene family, we calculated the Ka/Ks ratios of the 92 repeat pairs and found that they ranged from 0.071 (PfR2R3-MYB76/PfR2R3-MYB86) to 0.546 (PfR2R3-MYB58/PfR2R3-MYB103) (Supplementary Materials: Table S6). The ratio was much less than 1, which indicated that the PfR2R3-MYB gene family had undergone gene purification and selection in the evolutionary process.
To further investigate the evolution of the R2R3-MYB gene family in P. fortunei and A. thaliana, we constructed a collinearity map of 138 PfR2R3-MYBs and 126 AtR2R3-MYBs, which represented 182 PfR2R3-MYB orthologous gene pairs (Supplementary Materials: Figure S1). The results showed that the PfR2R3-MYBs had common ancestors and were highly conserved during the evolutionary process.

Analysis of Chromosomal Location and Synteny of PfR2R3-MYB Genes
To clarify the classification of PfR2R3-MYB proteins, we constructed an NJ rootless phylogenetic tree using 138 PfR2R3-MYB and 126 AtR2R3-MYB protein sequences from P. fortunei and A. thaliana ( Figure 2). The 138 PfR2R3-MYB proteins were divided into 27 subfamilies (named I-XXVII), most of which were consistent with those of A. thaliana in Figure 2. However, none of the PfR2R3-MYB proteins were found in subfamily XVII, suggesting that some R2R3-MYB genes may have been lost during the long-term evolution of the P. fortunei genome. Subfamilies I and XXIII each had twelve PfR2R3-MYBs, accounting for 8.63% of the total P. fortunei PfR2R3-MYB family, whereas subfamily XXVII had only one PfR2R3-MYB. Based on the evolutionary relationship of R2R3-MYB proteins in P. fortunei and poplar, we found that the 138R2R3-MYB and 192 MYB proteins mostly clustered in a branch, showing a sister-group relationship (Supplementary Figure S1).
The intron and exon positions of the PfR2R3-MYB genes were visualized using the GSDS software. Significant differences and diversity were observed in the PfR2R3-MYB genes' structure ( gesting that some R2R3-MYB genes may have been lost during the long-term evolution of the P. fortunei genome. Subfamilies I and XXIII each had twelve PfR2R3-MYBs, accounting for 8.63% of the total P. fortunei PfR2R3-MYB family, whereas subfamily XXVII had only one PfR2R3-MYB. Based on the evolutionary relationship of R2R3-MYB proteins in P. fortunei and poplar, we found that the 138R2R3-MYB and 192 MYB proteins mostly clustered in a branch, showing a sister-group relationship (Supplementary Figure S1).

Prediction of Cis-Regulatory Elements in PfR2R3-MYB Gene Family
To better understand the potential functions of the PfR2R3-MYBs, we detected cis-regulatory elements in the 2-kb upstream sequences of the PfR2R3-MYB genes using the PLACE database. We found that most of the cis-acting elements were light-responsive elements, whereas others were defense-and stress-responsive as well as hormone-responsive elements. Some common cis-acting elements were also found in the promoter regions of the PfR2R3-MYB genes (Figure 4, Supplementary Materials: Figure S3). The regulation of gene expression is mainly dependent on the presence of cis-regulatory elements in the gene promoter regions [32]. The presence of light-responsive, stress-responsive and hormoneresponsive elements in the upstream sequences leads us to speculate that they may have important roles in the PfR2R3-MYB gene family.

Prediction of Cis-Regulatory Elements in PfR2R3-MYB Gene Family
To better understand the potential functions of the PfR2R3-MYBs, we detected cisregulatory elements in the 2-kb upstream sequences of the PfR2R3-MYB genes using the PLACE database. We found that most of the cis-acting elements were light-responsive elements, whereas others were defense-and stress-responsive as well as hormone-responsive elements. Some common cis-acting elements were also found in the promoter regions of the PfR2R3-MYB genes (Figure 4, Supplementary Materials: Figure S3). The regulation of gene expression is mainly dependent on the presence of cis-regulatory elements in the gene promoter regions [32]. The presence of light-responsive, stress-responsive and hormone-responsive elements in the upstream sequences leads us to speculate that they may have important roles in the PfR2R3-MYB gene family.

PfR2R3-MYBs Involved in the Phytoplasma Interaction and Verification of qRT-PCR
We analyzed the expression patterns of the 138 PfR2R3-MYBs in phytoplasma-infected P. fortunei via RNA-Seq. The results showed that 119 PfR2R3-MYBs responded to phytoplasma infestation ( Figure 5a) and 26 (15 down-regulated and 11 up-regulated) of them were significantly differentially expressed in PF versus PFI (Table S7), indicating that these genes may actively respond to the PaWB phytoplasma. In previous studies, our research group found that the phytoplasma content decreased in PFI treated with 100 mg·L −1 of Rif with an extension of the treatment time, and the plant morphologically showed a reduction in branching and gradually reached a healthy state [28]. The 26 PfR2R3-MYBs actively responded to the PaWB phytoplasma in Rif 100-5/PFI, Rif100-10/PFI and Rif100-20/PFI (Figure 5a).
To validate the accuracy of the RNA-Seq results, we randomly selected eight PfR2R3-MYB differentially expressed genes in PFI vs. PF for the qRT-PCR analysis. The trends of their expression were similar, which indicated that the RNA-Seq data were reliable (Figure 5b-i).

PfR2R3-MYBs Involved in the Phytoplasma Interaction and Verification of qRT-PCR
We analyzed the expression patterns of the 138 PfR2R3-MYBs in phytoplasma-infected P. fortunei via RNA-Seq. The results showed that 119 PfR2R3-MYBs responded to phytoplasma infestation ( Figure 5a) and 26 (15 down-regulated and 11 up-regulated) of them were significantly differentially expressed in PF versus PFI (Table S7), indicating that these genes may actively respond to the PaWB phytoplasma. In previous studies, our research group found that the phytoplasma content decreased in PFI treated with 100 mg·L −1 of Rif with an extension of the treatment time, and the plant morphologically showed a reduction in branching and gradually reached a healthy state [28]. The 26 PfR2R3-MYBs actively responded to the PaWB phytoplasma in Rif 100-5/PFI, Rif100-10/PFI and Rif100-20/PFI (Figure 5a).
To validate the accuracy of the RNA-Seq results, we randomly selected eight PfR2R3-MYB differentially expressed genes in PFI vs. PF for the qRT-PCR analysis. The trends of their expression were similar, which indicated that the RNA-Seq data were reliable (Figure 5b-i).
To further study the functions of PfR2R3-MYB15, the positive seedlings were transplanted into pots. When the plants were cultured up to 90 days, the phenotypic characteristics showed that OE-PfR2R3-MYB15 had obvious branches compared with WT ( Figure  6). The results show that PfR2R3-MYB15 plays an important role in plant branching symptoms, which lays a foundation for future studies of the R2R3-MYB family in other plants.
To further study the functions of PfR2R3-MYB15, the positive seedlings were transplanted into pots. When the plants were cultured up to 90 days, the phenotypic characteristics showed that OE-PfR2R3-MYB15 had obvious branches compared with WT ( Figure  6). The results show that PfR2R3-MYB15 plays an important role in plant branching symptoms, which lays a foundation for future studies of the R2R3-MYB family in other plants. To further study the functions of PfR2R3-MYB15, the positive seedlings were transplanted into pots. When the plants were cultured up to 90 days, the phenotypic characteristics showed that OE-PfR2R3-MYB15 had obvious branches compared with WT ( Figure 6). The results show that PfR2R3-MYB15 plays an important role in plant branching symptoms, which lays a foundation for future studies of the R2R3-MYB family in other plants.

PfR2R3-MYB15 Located in the Nucleus
To verify the predicted nucleus location of the PfR2R3-MYB15 protein, 35S::PfR2R3-MYB15-GFP was transfected into N. benthamiana cells. The fluorescence signal of MYB15-GFP was detected in the nucleus, whereas GFP (the positive control) was in the full cell.
The results confirmed that PfR2R3-MYB15 was located in the nucleus (Figure 7), which was consistent with the Plant-mPLoc prediction.

PfR2R3-MYB15 Located in the Nucleus
To verify the predicted nucleus location of the PfR2R3-MYB15 protein, 35S::PfR2R3-MYB15-GFP was transfected into N. benthamiana cells. The fluorescence signal of MYB15-GFP was detected in the nucleus, whereas GFP (the positive control) was in the full cell. The results confirmed that PfR2R3-MYB15 was located in the nucleus (Figure 7), which was consistent with the Plant-mPLoc prediction.

The Discovery of Related Proteins with PfR2R3-MYB15 as Bait
PfR2R3-MYB15 self-activating activity was detected by Y2H screening. The results showed that 30 mM 3-AT can inhibit the autoactivation of PfR2R3-MYB15 (Figure 8a,c). To explore the molecular mechanism of PfR2R3-MYB15, a Y2H screen was performed to identify the potential binding partners of the PfR2R3-MYB15 protein. A total of 24 proteins from the cDNA library of PFI were found to interact with PfR2R3-MYB15 (Supplementary Materials: Tables S7). Among them, it has been verified by Y2H screening and BiFC that PFR2R3-MYB15 interacted with PfTAB2 (Paulownia_LG3G001014) (Figure 8c). ATAB2 was involved in the biogenesis of photosystem I (PSI) and II (PSII) [33]. We speculate that PfR2R3-MYB15 may be involved in the biosynthesis of the photosystem elements, but this needs to be further verified in future studies.

The Discovery of Related Proteins with PfR2R3-MYB15 as Bait
PfR2R3-MYB15 self-activating activity was detected by Y2H screening. The results showed that 30 mM 3-AT can inhibit the autoactivation of PfR2R3-MYB15 (Figure 8a,c). To explore the molecular mechanism of PfR2R3-MYB15, a Y2H screen was performed to identify the potential binding partners of the PfR2R3-MYB15 protein. A total of 24 proteins from the cDNA library of PFI were found to interact with PfR2R3-MYB15 (Supplementary Materials: Table S7). Among them, it has been verified by Y2H screening and BiFC that PFR2R3-MYB15 interacted with PfTAB2 (Paulownia_LG3G001014) (Figure 8c). ATAB2 was involved in the biogenesis of photosystem I (PSI) and II (PSII) [33]. We speculate that PfR2R3-MYB15 may be involved in the biosynthesis of the photosystem elements, but this needs to be further verified in future studies.

Discussion
R2R3-MYBs play important roles in plant growth and development, morphogenesis and resistance to stress [11][12][13]. The R2R3-MYB family has been widely studied in plants.
In A. thaliana, 126 R2R3-MYBs in 25 subfamilies have been identified and their functions have been predicted by phylogenetic and RNA-Seq analyses [2,34]. In maize, 157 R2R3-MYB TFs with diverse expression patterns and different functions have been reported [35]. In Solanum lycopersicum, 121 R2R3-MYB TFs were found, and the RNA interference of pGSMyb12 was shown to decrease the lateral meristems [36]. Despite their important roles in other plants, this study is the first to report R2R3-MYBs in the P. fortunei genome and analyze their physicochemical properties, structure and evolution. Furthermore, we

Discussion
R2R3-MYBs play important roles in plant growth and development, morphogenesis and resistance to stress [11][12][13]. The R2R3-MYB family has been widely studied in plants. In A. thaliana, 126 R2R3-MYBs in 25 subfamilies have been identified and their functions have been predicted by phylogenetic and RNA-Seq analyses [2,34]. In maize, 157 R2R3-MYB TFs with diverse expression patterns and different functions have been reported [35]. In Solanum lycopersicum, 121 R2R3-MYB TFs were found, and the RNA interference of pGSMyb12 was shown to decrease the lateral meristems [36]. Despite their important roles in other plants, this study is the first to report R2R3-MYBs in the P. fortunei genome and analyze their physicochemical properties, structure and evolution. Furthermore, we analyzed the expression patterns of the R2R3-MYB genes in PF and PFI and focused on the functional validation of the significantly overexpressed gene PfR2R3-MYB15. Our results provide new perspectives on the roles of PfR2R3-MYB TFs in plants' responses to PaWB infestations.
We identified 138 PfR2R3-MYB genes in the P. fortunei genome (Figure 1), which is more than those in A. thaliana [5], but less than those in Populus [8]. We speculated that the number of gene family members may be related to the genome size (Figures 1 and 2). Phylogenetic trees based on the R2R3-MYB protein sequences of P. fortunei and A. thaliana showed that the 138 PfR2R3-MYB genes clustered into 27 subfamilies (named I-XXVII), most of which contained A. thaliana R2R3-MYB TFs. Some exceptions were the XXI subfamily, which did not have any AtR2R3-MYBs, and the XVII subfamily, which did not have any PfR2R3-MYBs, implying that the R2R3-MYB gene family had undergone great differentiation in the evolution of these species. The most closely related genes likely have similar functions. For example, in A. thaliana, the overexpression of AtRAX2 resulted in an increased number of branches compared with the wild type, whereas the number of branches decreased in the AtRAX2 mutant. In the XXIII subfamily, PfR2R3-MYB15 was most closely related to At2g36890 (RAX2), suggesting that PfR2R3-MYB15 may have similar functions to AtRAX2.
Plant genomes are known to undergo whole-genome replication, segmental replication and tandem repetition events which contribute to the expansion of distinct gene families. Gene duplication involves three mechanisms: tandem replication, segmental duplication and transposition events [37]. Fragment and tandem duplication are the main evolutionary mechanisms for gene family expansion and evolution in plants [37]. In this study, the chromosome distribution analysis showed that the PfR2R3-MYBs were unevenly distributed on the 20 chromosomes. The collinearity analysis indicated that 92 fragment duplication events had occurred. These results suggested that fragment duplication played an important role in the evolution of the PfR2R3-MYB gene family. We also found that the Ka/Ks ratios of the P. fortunei segmental repeat pairs were much less than one, implying that the PfR2R3-MYB family experienced strong purification selection during its evolution. The collinear relationship of the R2R3-MYB genes in A. thaliana and P. fortunei showed that 105 PfR2R3-MYBs and 99 AtR2R3-MYBs formed 182 lineal homologous gene pairs (Figure 1c). AtRAX1 (MYB37) proteins interact with CUC2 to further regulate axillary meristems' expression and control lateral bud production [38]. RAX1, RAX2 (AtMYB38) and RAX3 (AtMYB84) have obvious effects on branch development in A. thaliana. Plants overexpressing RAX2 showed lateral bud symptoms, whereas RAX2 knockout mutants had fewer lateral buds than the wild type. Furthermore, the knockout line RAX1/RAX2/RAX3 was shown to have significantly reduced branches [19]. Together, our results provide new ideas for studying the evolutionary relationships of the PfR2R3-MYB gene family in P. fortunei.
Numerous studies have shown that members of the PfR2R3-MYB family play important roles in plant branching, so we analyzed the expression patterns of the PfR2R3-MYB genes in PF and PFI by RNA-Seq. We detected 119 PfR2R3-MYB genes that responded to PaWB phytoplasma infestationS, and 26 (15 down-regulated and 11 up-regulated) of them were significantly differentially expressed in PF versus PFI (FC ≥ 2.0, false discovery rate < 0.01) (Figure 5a). We focused on PfR2R3-MYB15, which was significantly up-regulated in PF versus PFI. To confirm its function, PfR2R3-MYB15 was cloned and transfected into P. trichocarpa, and its overexpression induced lateral bud symptoms (Figure 6b). Y2H and BiFC experiments demonstrated that PfR2R3-MYB15 interacted with PfTAB2 (Paulownia_LG3G001014). Barneche et al. (2006) found PfTAB2 was the localization of the chloroplast. It may be due to the localization of PfTAB2 that PfR2R3-MYB15 interacted with PfTAB2 in the cell membrane and nucleus via BiFC. Previous studies have shown that ATAB2 is involved in the biogenesis of photosystems I (PSI) and II (PSII) among the PfRAX2-related proteins that are involved in photosynthesis [33]. In this study, we speculate that PfR2R3-MYB15 might participate in the biosynthesis of the photosystem by influencing PfTAB2. These results provide a valuable resource to explore the potential functions of PfR2R3-MYB genes and provide a foundation for future studies of the R2R3-MYB TF family in Paulownia and other plants.

Conclusions
In this study, we identified 138 P. fortunei R2R3-MYB TFs, which were distributed on 20 chromosomes and were divided into 27 groups. Their gene structure, phylogenetic relationship, cis-regulatory elements, evolutionary relationships and rich segmental replication events were analyzed. By combining RNA-seq data, we profiled the expression patterns of the 138 PfR2R3-MYB members in phytoplasma-infected seedlings of P. fortunei. We found that PfR2R3-MYB15 is significantly up-regulated in infected P. fortunei. The overexpression of PfR2R3-MYB15 induced branching symptoms in P. trichocarpa. The yeast two-hybrid screening and bimolecular fluorescence complementation showed that PfR2R3-MYB15 interacted with PfTAB2. Our results will provide a foundation for future studies of the R2R3-MYB TF family in Paulownia.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/genes14010007/s1, Table S1: The primers used for qRT-PCR in PF and PFI. Table S2: Primers of 35S::PfR2R3-MYB15-GFP and 35S::GFP were designed. Table S3: Primers of Paulownia_LG3G001014 were designed on Primer Premier 5.0. Table S4: List of 192 identified P. fortunei R2R3-MYB transcription factors. Table S5: Identification of R2R3-MYB gene family members of Paulownia and analysis of their physicochemical properties. Table S6: The KA/KS analysis of PfR2R3-MYB genes. Table S7: 26 PfR2R3-MYB genes respond to phytoplasma infestations in PF versus PFI. Figure S1: Synteny analysis of PfR2R3-MYB genes between A. thaliana and P. fortunei. The blue line represents the gene pair of PfR2R3-MYB. Figure S2: The predictive motif information of R2R3-MYB transcription factors in P. fortunei. Figure S3: Cis-regulatory elements analysis of the promoter region of PfR2R3-MYB genes.  Data Availability Statement: All the sequencing data from P. fortunei transcriptome used in this study are available from the SRA-Archive (http://www.ncbi.nlm.nih.gov/sra accessed on 27 June 2021) of NCBI under the study accession numbers SRR11787883, SRR11787894, SRR11787905 and SRR11787912-SRR11787921.