Genome-Wide Identiﬁcation and Expression Analysis of MYB Transcription Factors and Their Responses to Abiotic Stresses in Woodland Strawberry ( Fragaria vesca )

: Woodland strawberry ( Fragaria vesca ) is a diploid strawberry that is widely used as a model of cultivated octoploid strawberry ( Fragaria × ananassa ). It has also been used as a model for Rosaceae fruits, non-climacteric fruits, and stolons. The MYB superfamily is the largest transcription factor family in plants, and its members play important roles in plant growth and development. However, the complete MYB superfamily in woodland strawberry has not been studied. In this study, a total of 217 MYB genes were identiﬁed in woodland strawberry and classiﬁed into four groups: one 4R-MYB protein, ﬁve 3R-MYB proteins, 113 2R-MYB proteins, and 98 1R-MYB proteins. The phylogenetic relationship of each MYB subgroup was consistent in terms of intron/exon structure and conserved motif composition. The MYB genes in woodland strawberry underwent loss and expansion events during evolution. The transcriptome data revealed that most FveMYB genes are expressed in several organs, whereas 15 FveMYB genes exhibit organ-speciﬁc expression, including ﬁve genes (FveMYB101, -112, -44, and -8; FveMYB1R81) in roots, two genes (FveMYB62 and -77) in stolon tips, three genes (FveMYB99 and -35; FveMYB1R96) in open ﬂowers, and ﬁve genes (FveMYB76 and -100; FveMYB1R4, -5, and -86) in immature fruits. During fruit ripening of woodland strawberry, the expression levels of 84 FveMYB genes were decreased, of which ﬁve genes (FveMYB4, -22, -50, and -66; FveMYB1R57) decreased more than 10-fold, whereas those 18 FveMYB genes were increased, especially FveMYB10 and FveMYB74 increased more than 30-fold. In addition, the expression levels of 36, 68, 52, and 62 FveMYB genes were altered by gibberellic acid, abscisic acid, cold, and heat treatments, respectively, and among them, several genes exhibited similar expression patterns for multiple treatments, suggesting possible roles in the crosstalk of multiple signaling pathways. This study provides candidate genes for the study of stolon formation, fruit development and ripening, and abiotic stress responses. acid (ABA) treatments) were explored. This study provides basic information on the woodland strawberry MYB superfamily as well as important clues for studying the functions of possible MYB genes involved in stolon formation, fruit development and ripening, and responses to abiotic stresses in woodland strawberry.


Introduction
Higher plants have evolved complex regulatory mechanisms to adapt to a constantly changing environment. Transcription factors are important components of these regulatory mechanisms due to their ability to regulate cell development, differentiation, and growth by binding to specific DNA sites to regulate gene expression. According to the type of DNA binding domain, transcription factors can be classified into different gene families, and among them, the MYB transcription factors comprise an important family containing numerous members [1][2][3].
MYB transcription factors are characterized by one to four conserved MYB domains (Myb_DNA-binding domains). Each domain contains approximately 52 amino acids and consists of three alpha helices, of which, helices 2 and 3 form a helix-turn-helix structure that inserts into the main groove of the target DNA that is recognized by helix 3 [4]. According to the number of MYB domain repeats, the MYB transcription factors can be classified into four subfamilies. The 2R-MYB proteins, also called R2R3-MYBs, contain two MYB domain repeats (generally R2 and R3 repeats). The 3R-MYB or R1R2R3-MYB proteins contain three MYB domain repeats (generally R1, R2, and R3). The 4R-MYB proteins contain four MYB domain repeats, and 1R-MYB or MYB-related proteins contain a single or partial MYB domain repeat [5].
In higher plants, R2R3-MYBs are the most widespread MYB proteins [6]. R2R3-MYB proteins not only play important roles in phenylpropanoid metabolism [7], flavonoid biosynthesis [8], anthocyanin biosynthesis [9], indolic and aliphatic glucosinolate biosynthesis [10], epidermis and root formation [11], and embryogenesis [12], but also protect plants against environmental stresses such as high salinity, cold, high temperature, drought, phosphate starvation, nitrogen deficiency, and so on [13]. The MYB-related subfamily may be older than the R2R3-MYB subfamily in terms of their evolutionary process [2]. MYBrelated proteins may have different functions compared to R2R3-MYB and R1R2R3-MYB proteins due to their different manner of binding to target DNA [2]. MYB-related subfamily genes have attracted little attention, and few have been functionally studied; for example, some Circadian Clock Associated 1 (CCA1)-like genes are involved in maintaining the circadian rhythm [14][15][16], and some Caprice (CPC)-like genes are involved in controlling cell morphogenesis [17][18][19]. In animals, MYB proteins are encoded by a small gene family, usually 3R-MYBs [20]. As homologs, plant 3R-MYB proteins have been reported to participate in cell-cycle regulation and plant abiotic stress responses [21][22][23]. At present, the functions of 4R-MYB proteins remain unclear.
In Arabidopsis, MYB transcription factors have been investigated extensively. MYBs are widely involved in secondary metabolism, cell morphogenesis, signal transduction in plant growth, abiotic and biological stress responses, circadian rhythms, dorso-ventral relationships, and other physiological and biochemical activities [24]. Some MYB genes in other plant species are thought to be related to important agricultural features, such as reproductive development, wood formation, legume-specific nodules, and wine quality in grapes [25][26][27]. In addition, numerous MYB superfamily proteins in rice, cucumber, pear, peach, apple, and other plants have been defined and described as being involved in abiotic stress responses, hormones, anthers, and fruit development [27][28][29][30][31][32].
Cultivated strawberry is not only a popular fruit but also a typical non-climacteric fruit that has representative significance in the study of fruit development and ripening. As a diploid plant, the woodland strawberry genome is less complex than the cultivated strawberry genome, and it has become a model plant for studying strawberries, Rosaceae, and other non-climacteric fruits [33]. At present, there are limited reports on the woodland strawberry MYB superfamily, only the identification of R2R3-MYB subfamily members and a few functional studies [20,34,35]. Therefore, in this study, we comprehensively analyzed the entire MYB superfamily of woodland strawberry, including the less studied 1R-MYB, 3R-MYB, and 4R-MYB families. The phylogeny, gene structure, and chromosomal locations of these putative MYB genes were analyzed. In addition, the expression patterns of all MYB genes in various organs and tissues of woodland strawberry and seedlings under different treatments (cold and heat stress, and exogenous gibberellic acid (GA) and abscisic acid (ABA) treatments) were explored. This study provides basic information on the woodland strawberry MYB superfamily as well as important clues for studying the functions of possible MYB genes involved in stolon formation, fruit development and ripening, and responses to abiotic stresses in woodland strawberry.

Identification and Phylogenetic Analysis of MYB Proteins
The protein database and annotation information of woodland strawberry was downloaded from GDR (https://www.rosaceae.org/species/fragaria/fragaria_vesca (access Horticulturae 2021, 7, 97 3 of 21 on 13 May 2019)). First, the MYB DNA-binding domain HMM file (PF00249) from the Pfam database (https://pfam.xfam.org/(accessed on 13 May 2019)) was used as a query to look for possible MYB proteins in woodland strawberry. Secondly, each candidate MYB gene containing MYB DNA-binding domain was confirmed through the pfam website (https://pfam.xfam.org/(accessed on 13 May 2019)) with an E-value cut-off of 1 × e −4 . The longest protein was selected when a gene had alternative splicing. Finally, 217 MYB proteins in woodland strawberry were identified.
198 Arabidopsis MYB gene names were issued by Chen et al. [36]. The corresponding Arabidopsis protein sequences were downloaded from Phytozome (https://phytozome.jgi. doe.gov/pz/portal.html (accessed on 13 May 2019)). In order to analyze the phylogenetic relationship and evolutionary history of the MYB gene family, four phylogenetic trees were set up in this paper. The protein sequence alignment was performed through ClustalX2.0 software with default settings to obtain faste files. The phylogenetic trees were created by MEGA 6.0 using the ML (Maximum Likelihood) method, and the bootstrap test was repeated 1000 times.

Chromosomal Location, Gene Structure, Conserved Domain, and Motif Analysis of MYB Genes
Conserved domains of woodland strawberry MYB proteins were obtained through the pfam website using Batch search (https://pfam.xfam.org/search#tabview=tab1 (accessed on 13 May 2019)). The conserved motifs were analyzed through MEME (http://memesuite.org/tools/meme (accessed on 13 May 2019)), with the maximum number of motifs was 15, and other parameters default. The results of conserved domains and motifs were visualized through TBtools software. The chromosomal location, duplication events, and intron/exon structure of the woodland strawberry MYB genes were analyzed by utilizing TBtools software (Version No. 0.66763, South China Agricultural University, Guangzhou, China) [37].

Plant Materials and Treatments
The genome-sequenced diploid woodland strawberry "Hawaii 4" (Fragaria vesca) was used as plant material. The plants were grown in pots in a glass-enclosed greenhouse under native conditions. Immature fruits with seeds were collected at about 8-15 DPA, and roots, crowns, leaves, petioles, stolons, stolon tips, and fully open flowers were also sampled at the same time. All of the samples were rapidly frozen in liquid nitrogen, and then stored at −80 • C for RNA extraction. Roots, crowns, leaves, petioles, stolons, and immature fruits had three biological replicates, and stolon tips and fully open flowers had two biological replicates because of limited samples.
For cold, heat, GA 3 and ABA treatments, 'Hawaii 4' seedlings were grown in a climate chamber under a 16 h light (22 • C) with 15,000 Lux irradiance and 8 h dark (20 • C) for about two months. To avoid photoperiod effect on gene expression during treatments, the growth conditions were reset as 24 h light (22 • C), and the seedlings were re-adaptation for three days before treatments. Afterwards, seedlings were transferred to another chamber maintained at 40 • C for heat treatment and at 4 • C for cold treatment. For GA 3 and ABA treatment, the seedlings were sprayed with 50 mg/L GA 3 and 100 µM ABA solution, respectively. At 0, 2, 6, 12, 24, and 48 h after treatments, the seedlings were collected, rapidly frozen in liquid nitrogen, and then stored at −80 • C for RNA extraction. Three biological replicates were performed.

Heat Map Construction of FveMYB Gene Expressions
The reads per kilobase per million (RPKM) values of transcriptome data in Figure  6B,D of woodland strawberry were downloaded from Li et al. (2019) [38]. The FPKM values of FveMYB genes in Figure 6A,C, 8-11 were selected from our own transcriptome data (accession number: PRJNA723745). Then RPKM and FPKM values were transformed in log2 level, and heat maps were constructed by MeV4.8 software (Dana-Farber Cancer Institute, Boston, MA, USA).

Identification and Classification of MYB Proteins in Woodland Strawberry
The Myb_DNA-binding (PF00249) hidden Markov model profile was used as a query to search for potential MYB proteins in the woodland strawberry protein database. As a result, 217 proteins containing MYB domains were identified in woodland strawberry. According to the number of Myb_DNA-binding repeats, we categorized these proteins into four groups: one 4R-MYB, five 3R-MYBs, 113 2R-MYBs, and 98 1R-MYBs. Accordingly, these MYB proteins were named FveMYB4R1, FveMYB3R1-5, FveMYB1-113, and FveMYB1R1-98, respectively (Table S1).

Phylogenetic Analysis and Gene Structure of MYB Genes in Woodland Strawberry
According to a previous study in watermelon [39], the 1R-MYB subfamily exhibits greater differentiation and evolution capabilities than do the 2R-MYB and 3R-MYB subfamilies in plants. Thus, we used 2R-MYB and 3R-MYB proteins to construct one phylogenetic tree, while the 1R-MYB and 4R-MYB proteins were combined to construct another phylogenetic tree. Based on the bootstrap values and topology, the two phylogenetic trees are divided into 29 and 14 subgroups, respectively; these subgroups were named S1-29 and Sr1-14 ( Figures 1A and 2A).
To investigate the genetic structural diversity of FveMYBs, we analyzed the conserved motif structure and intron-exon organization. The motif compositions of most MYB proteins are similar in the 2R-MYB and 3R-MYB groups, but diverse in the 1R-MYB and 4R-MYB groups ( Figures 1B and 2B). In the 2R-MYB group, almost all proteins in subgroups S2-23 have Motifs 3, 5, 2, 4, and 1, except for FveMYB105 and FveMYB9, which lack Motif 5. In subgroup S1, the members have another conserved motif (Motif 11). S24 members have Motifs 3, 2, 4, and 1 as well as three specific motifs (Motif 15, 13, and 12). S26 and S27 members have Motifs 3, 9, 2, 4, and 1. S29 members have one Motif 8, one Motif 9, and two Motif 3s. S25 members, which belong to the 3R-MYB subfamily, have Motifs 3, 9, 2, 4, and 1 as well as one specific Motif 10 ( Figure 1B). In the 1R-MYB group, the motif composition is much more diverse among the 14 subgroups but is similar within each subgroup with one or two specific motifs. For example, Sr2 contains Motif 10, Sr4 contains Motif 4, Sr5 and Sr6 contain Motif 3, Sr10 contains Motif 9, Sr12 contains Motifs 14 and 13, Sr13 contains Motif 15, and so on ( Figure 2B). In addition, the motifs of most 2R-MYB proteins are usually located at the N terminus, whereas the motifs of 1R-MYBs are more scattered ( Figures 1B and 2B).
The intron-exon organization of 1R-MYB genes is much more variable than that of 2R-MYB genes ( Figures 1C and 2C);~66% of 2R-MYB genes have a conserved splicing pattern showing a structure of three exons separated by two introns ( Figure 1C), whereas~61% of 1R-MYBs have at least four introns, with a maximum number of 19 introns ( Figure 2C). Six 2R-MYB genes (FveMYB94, -110, -26, -25, -4, and -14) and two 1R-MYB genes (FveMYB1R79 and -90) have no introns ( Figures 1C and 2C). These findings imply that 2R-MYB genes are conservative, whereas 1R-MYB genes are diverse in their evolutionary processes. The five 3R-MYB genes contain at least five introns, with a maximum number of 11 introns. The only 4R-MYB gene (FveMYB4R1) contains eight introns. The degree of variable splicing in the genes of these two subfamilies is higher than that of the 2R-MYB subfamily, and closer to the 1R-MYB subfamily. The phylogenetic tree of five 3R-MYB proteins and 113 2R-MYB proteins in woodland strawberry. According to previous study in soybean [26], the phylogenetic tree was divided into 29 subgroups (S1-S29) based on bootstrap values (>50), which had only one R1R2R3 subgroup. Some proteins were not clustered together due to bootstrap values were below than 50. (B) The motif composition of 2R-MYB and 3R-MYB proteins from woodland strawberry. Fifteen motifs were marked with different color boxes. (C) The domain composition and intron-exon organization of 2R-MYB and 3R-MYB proteins from woodland strawberry. Red boxes represent Myb_DNA-binding domains. The exon structure was shown as yellow CDS (s) and green UTR (s). The space (grey lines) between the boxes were intron(s). According to previous study in soybean [26], the phylogenetic tree was divided into 29 subgroups (S1-S29) based on bootstrap values (>50), which had only one R1R2R3 subgroup. Some proteins were not clustered together due to bootstrap values were below than 50. (B) The motif composition of 2R-MYB and 3R-MYB proteins from woodland strawberry. Fifteen motifs were marked with different color boxes. (C) The domain composition and intron-exon organization of 2R-MYB and 3R-MYB proteins from woodland strawberry. Red boxes represent Myb_DNA-binding domains. The exon structure was shown as yellow CDS (s) and green UTR (s). The space (grey lines) between the boxes were intron(s).
Horticulturae 2021, 7, x FOR PEER REVIEW 6 of 22 To evaluate the phylogenetic relationships of FveMYB proteins in woodland strawberry and Arabidopsis, 2R-MYB and 3R-MYB proteins from woodland strawberry and Arabidopsis were used to construct one phylogenetic tree, while 1R-MYB and 4R-MYB proteins were used to build another phylogenetic tree. Except for 14 2R-FveMYB proteins, one 2R-AtMYB protein, and one 3R-FveMYB protein, the remaining 2R-MYB and 3R-MYB proteins could be grouped into 44 clades, designated C1-44, based on bootstrap values (>50). Clade C6 contained all 3R-MYB proteins, and the remaining 43 clades contained 224 of the 239 2R-MYB proteins ( Figure 3). Among the 43 clades, 30 clades contained 2R-MYB proteins from both woodland strawberry and Arabidopsis, indicating that the possible functions of these FveMYBs were similar to those of Arabidopsis AtMYBs ( Figure 3). However, clades C9, C38, C39, and C43 contained only AtMYBs, whereas C1, C20, C22, C23, C27, and C31 contained only FveMYBs, suggesting that these To evaluate the phylogenetic relationships of FveMYB proteins in woodland strawberry and Arabidopsis, 2R-MYB and 3R-MYB proteins from woodland strawberry and Arabidopsis were used to construct one phylogenetic tree, while 1R-MYB and 4R-MYB proteins were used to build another phylogenetic tree. Except for 14 2R-FveMYB proteins, one 2R-AtMYB protein, and one 3R-FveMYB protein, the remaining 2R-MYB and 3R-MYB proteins could be grouped into 44 clades, designated C1-44, based on bootstrap values (>50). Clade C6 contained all 3R-MYB proteins, and the remaining 43 clades contained 224 of the 239 2R-MYB proteins ( Figure 3). Among the 43 clades, 30 clades contained 2R-MYB proteins from both woodland strawberry and Arabidopsis, indicating that the possible functions of these FveMYBs were similar to those of Arabidopsis AtMYBs ( Figure 3). However, clades C9, C38, C39, and C43 contained only AtMYBs, whereas C1, C20, C22, C23, C27, and C31 contained only FveMYBs, suggesting that these proteins may play different roles in woodland strawberry and Arabidopsis. In addition, clade C6 contained five 3R-MYB proteins from Arabidopsis, but only four of the five 3R-MYB proteins from woodland strawberry belonged to this clade; the remaining 3R-MYB protein (FveMYB3R5) did not belong to any clade, suggesting a different role for FveMYB3R5 ( Figure 3). proteins may play different roles in woodland strawberry and Arabidopsis. In addition, clade C6 contained five 3R-MYB proteins from Arabidopsis, but only four of the five 3R-MYB proteins from woodland strawberry belonged to this clade; the remaining 3R-MYB protein (FveMYB3R5) did not belong to any clade, suggesting a different role for FveMYB3R5 ( Figure 3). Arabidopsis. According to previous study in soybean [26], the phylogenetic tree was divided into 44 phylogenetic clades (C1-C44) based on bootstrap values (>50). Red dots represent MYB proteins from woodland strawberry. Grey dots represent MYB proteins from Arabidopsis. Blue stars represent MYB proteins from woodland strawberry not fit into any clades (bootstrap support values <50). S1-S29 are group numbers of strawberry 2R-MYB and 3R-MYB proteins named in Figure 1. The function of some members from Arabidopsis was annotated outside the circle.

Chromosomal Distribution and Gene Duplication of MYB Genes in Woodland Strawberry
Using woodland strawberry gene annotation, the chromosomal distribution of the 217 FveMYB genes was studied. The results revealed that the 2R-MYB and 1R-MYB genes are unevenly distributed over all seven chromosomes ( Figure 5)  Arabidopsis. According to previous study in soybean [26], the phylogenetic tree was divided into 23 phylogenetic clades (Cr1-Cr23) based on bootstrap values (>50). Red dots represent MYB proteins from woodland strawberry. Grey dots represent MYB proteins from Arabidopsis. Blue stars represent MYB proteins from woodland strawberry not fit into any clades (bootstrap support values <50). Sr1-Sr14 are group numbers of strawberry 1R-MYB and 4R-MYB proteins named in Figure 2. The function of some members from Arabidopsis was annotated outside the circle.

Chromosomal Distribution and Gene Duplication of MYB Genes in Woodland Strawberry
Using woodland strawberry gene annotation, the chromosomal distribution of the 217 FveMYB genes was studied. The results revealed that the 2R-MYB and 1R-MYB genes are unevenly distributed over all seven chromosomes ( Figure 5) (Fvb2, 4, 5, 6, and 7). The only strawberry 4R-MYB gene is located in Fvb5 ( Figure 5). In summary, strawberry Fvb5 and Fvb6 contain the largest numbers of FveMYB genes. tionary history of the MYB superfamily in woodland strawberry, the duplication type of FveMYB genes was studied. The results revealed five 2R-MYB gene clusters (Fvb2: FveMYB19-22; Fvb5: FveMYB54 and -55, FveMYB61 and -62, and FveMYB66 and -67; Fvb6: FveMYB91 and -92) and one 1R-MYB gene cluster (Fvb1: FveMYB1R4 and -R5) in the seven woodland strawberry chromosomes ( Figure 5). However, no segment duplication of the FveMYB genes was found through genome investigation.

Expression Profiles of FveMYB Genes in Various Organs of Woodland Strawberry
To investigate the expression patterns of FveMYB genes in different organs of woodland strawberry, we analyzed transcription data from the root, crown, stolon, stolon tip, leaf, petiole, open flower, and immature fruit of woodland strawberry variety "Hawaii 4". Genes with fragments per kilobyte per million (FPKM) values >1 were considered to be expressed [40,41]. Twenty-one 2R-MYB genes had FPKM values <1 in all eight organs, suggesting that these genes may not have roles in these organs. Twenty-nine 2R-MYB genes had FPKM values >1 in all eight organs, suggesting widespread function of these genes in organ development. Sixty-three 2R-MYB genes were expressed in one to seven organs; notably, 10 of them were specifically expressed in certain woodland strawberry organs ( Figure 6A; Table S2  Gene expansion plays an important role in species evolution. To examine the evolutionary history of the MYB superfamily in woodland strawberry, the duplication type of FveMYB genes was studied. The results revealed five 2R-MYB gene clusters (Fvb2: FveMYB19-22; Fvb5: FveMYB54 and -55, FveMYB61 and -62, and FveMYB66 and -67; Fvb6: FveMYB91 and -92) and one 1R-MYB gene cluster (Fvb1: FveMYB1R4 and -R5) in the seven woodland strawberry chromosomes ( Figure 5). However, no segment duplication of the FveMYB genes was found through genome investigation.

Expression Profiles of FveMYB Genes in Various Organs of Woodland Strawberry
To investigate the expression patterns of FveMYB genes in different organs of woodland strawberry, we analyzed transcription data from the root, crown, stolon, stolon tip, leaf, petiole, open flower, and immature fruit of woodland strawberry variety "Hawaii 4". Genes with fragments per kilobyte per million (FPKM) values >1 were considered to be expressed [40,41]. Twenty-one 2R-MYB genes had FPKM values <1 in all eight organs, suggesting that these genes may not have roles in these organs. Twenty-nine 2R-MYB genes had FPKM values >1 in all eight organs, suggesting widespread function of these genes in organ development. Sixty-three 2R-MYB genes were expressed in one to seven organs; notably, 10 of them were specifically expressed in certain woodland strawberry organs ( Figure 6A; Table S2). For example, FveMYB101, -112, -44, and -8 were only expressed in roots; FveMYB62 and -77 were expressed in stolon tips; FveMYB99 and -35 were expressed in open flowers; and FveMYB76 and -100 were mainly expressed in immature fruits. These results suggest that these genes function in the development of specific organs. In addition, FveMYB6 and -30 were expressed in the six vegetative organs but were hardly expressed in the two reproductive organs. By contrast, FveMYB55 and -90 were mainly expressed in the two reproductive organs but were almost undetectable in the six vegetative organs ( Figure 6A; Table S2).  Table S2. RPKM values in (B,D) were located in Supplementary S3. The heat map showed log2 level.

Expression Analysis of MYB Genes in Woodland Strawberry under ABA and GA Treatments
The plant hormone ABA is a regulator of strawberry fruit maturation and abiotic stress responses [42,43], whereas GAs have been reported to be involved in strawberry stolon formation and early fruit development [42,44]. Therefore, we analyzed transcrip- Transcription data for the 1R-MYB genes showed that 62% (61/98) of the genes were strongly expressed in all eight organs. Eleven 1R-MYB genes had FPKM values <1 in all eight organs, and the remaining 26 1R-MYB genes were expressed in one to seven organs ( Figure 6C; Table S2). Five 1R-MYB genes exhibited organ-specific expression. For example, FveMYB1R96 was mainly expressed in open flowers, whereas FveMYB1R81 was only expressed in roots. Three 1R-MYB genes (FveMYB1R4, -5, and -86) were mainly expressed in immature fruits. FveMYB1R73 and -78 were preferentially expressed in vegetative organs. In addition, three 3R-MYB genes (FveMYB3R2, -4, and -5) had FPKM values <1 in all eight organs ( Figure 6A; Table S2); FveMYB3R1 was expressed in all eight organs, with the highest expression in the stolon tip, and FveMYB3R3 was expressed in multiple organs, but the FPKM values were low. The only 4R-MYB gene was expressed in all eight organs, indicating that FveMYB4R1 may play a role in the development of multiple organs ( Figure 6C; Table S2). In conclusion, FveMYB genes may be associated with the growth and development of various organs in woodland strawberry.

Expression Analysis of MYB Genes in Woodland Strawberry under ABA and GA Treatments
The plant hormone ABA is a regulator of strawberry fruit maturation and abiotic stress responses [42,43], whereas GAs have been reported to be involved in strawberry stolon formation and early fruit development [42,44]. Therefore, we analyzed transcription data from woodland strawberry seedlings treated with exogenous ABA and GA3 to explore the response of FveMYBs to ABA and GA treatments.
The expression levels of 39 2R-MYB genes and 29 1R-MYB genes were altered after exogenous ABA treatment (Figure 7A,C; Table S4). Among them, 20 2R-MYB genes (  Table S4). In addition, the expression levels of two 2R-MYB genes (FveMYB4 and -53) and two 1R-MYB genes (FveMYB1R57 and -68) continuously increased after ABA treatment ( Figure 7A,C; Table S4). By contrast, the expression levels of two 1R-MYB genes (FveMYB1R6 and -92) continuously decreased after ABA treatment. In addition, one 1R-MYB gene (FveMYB1R20) had a FPKM value >5.06 before treatment, but this value decreased to <1 after ABA treatment ( Figure 7C; Table S4). In summary, the results suggest that many FveMYB genes are regulated by ABA and thus, may participate in ABA-related biological processes in woodland strawberry.  Supplementary Table S4, and the heat map showed log2 level.

Expression Patterns of FveMYB Genes under Heat and Cold Treatments
As a perennial plant, strawberry is subjected to the heat of summer and the cold of winter; thus, we analyzed the response of woodland strawberry MYB genes to cold and heat stresses. The transcription data revealed that a large number of FveMYB genes in woodland strawberry responded to cold and heat stresses. The detailed results follow.

Expression Patterns of FveMYB Genes under Heat and Cold Treatments
As a perennial plant, strawberry is subjected to the heat of summer and the cold of winter; thus, we analyzed the response of woodland strawberry MYB genes to cold and heat stresses. The transcription data revealed that a large number of FveMYB genes in woodland strawberry responded to cold and heat stresses. The detailed results follow.

Identification and Comparison of MYB Genes in Woodland Strawberry and Arabidopsis
The MYB superfamily is a large transcription factor family in plants. At present, numerous MYB genes with various numbers of MYB domain repeats have been identified in different plants, such as 680 members in Brassica napus, 256 members in peach, 198 members in Arabidopsis, etc. [32,36,45]. By contrast, fewer MYB genes have been identified in some species, such as 162 members in watermelon and 122 members in Brachypodium distachyon [39,46]. The number of MYB genes may be linked to the genome size and

Identification and Comparison of MYB Genes in Woodland Strawberry and Arabidopsis
The MYB superfamily is a large transcription factor family in plants. At present, numerous MYB genes with various numbers of MYB domain repeats have been identified in different plants, such as 680 members in Brassica napus, 256 members in peach, 198 members in Arabidopsis, etc. [32,36,45]. By contrast, fewer MYB genes have been identified in some species, such as 162 members in watermelon and 122 members in Brachypodium distachyon [39,46]. The number of MYB genes may be linked to the genome size and evolution of the species. In strawberry, although the R2R3-MYB subfamily had been previously identified, the other MYB subfamilies and their expression patterns have not been comprehensively analyzed [47]. In this study, we identified 217 FveMYB genes in woodland strawberry (Table S1). Similar to Arabidopsis and rice [28], 113 (52%) and 98 (45%) FveMYBs were identified as 2R-MYB and 1R-MYB genes, respectively (Table S1). In addition, we found six gene clusters (five 2R-MYB gene clusters and one 1R-MYB gene cluster) and no segment duplications in the woodland strawberry genome ( Figure 5). This is similar to the duplication type of watermelon MYB genes [39].
The 3R-MYB proteins are universal in vertebrates but rare in plants. In Arabidopsis and rice, five and four 3R-MYB proteins have been identified, respectively [28]. Similarly, we also found five 3R-MYB proteins in woodland strawberry. Interestingly, only four strawberry 3R-MYB proteins (FveMYB3R1-4) and five Arabidopsis 3R-MYB proteins (AtMYB3R1-5) are classified within the same clade (C6) (Figure 3). Based on motif structure, the conserved motif of FveMYB3R5 is significantly different from that of FveMYB3R1-4. For example, FveMYB3R5 additionally has Motifs 2, 4, and 1, whereas FveMYB3R1-4 have Motifs 10 and 9, which FveMYB3R5 lacks (Figure 1). At the same time, the sequence alignment showed that the three repeats of FveMYB3R1-4 and AtMYB3R1-5 are very similar, whereas FveMYB3R5 only has one repeat ( Figure S1). Therefore, we speculate that FveMYB3R5 has a different function from those of other woodland strawberry 3R-MYB proteins. In addition, one woodland strawberry 4R-MYB protein (FveMYB4R1) was identified that clustered with Arabidopsis 4R-MYB proteins, suggesting a similar function in these two plants.

Possible Roles of FveMYB Genes and Their Crosstalk with the ABA Signaling Pathway in Response to Abiotic Stresses
High temperatures in summer and low temperatures in winter greatly impact the vegetative and reproductive growth of cultivated strawberry. To identify the FveMYB genes that may be linked to cold and heat tolerance, we analyzed the expression changes in FveMYB genes under cold and heat treatments. In plants, many MYB genes are associated with the response to abiotic stresses, and the Arabidopsis MYB genes have been the most extensively investigated. For example, AtMYB21, -33, -35, -60, -65, -80, -99, and -110 participate in drought stress; AtMYB20, AtMYB73, and HPPBF-1 (AT5G13820) are involved in salt stress; PHR1 (AT1G12370), PHR2 (At1g79430), PHL1 (AT5G29000), and CCA1 (AT2G46830) participate in phosphate (Pi) starvation; RVE1 (AT5G17300), MYBC1 (AT2G40970), CIR1 (AT5G37260), AtMYB15, and AtMYB44 respond to freezing stress or cold; the expression of AtMYB3, AtMYB6, and MYBL2 (AT1G71030) is regulated by high temperatures, and so on [13]. In our study, 73% (159/217) and 70% (151/217) of woodland strawberry MYB genes were induced by cold and heat stress, respectively. Among them, FveMYB1R61 was upregulated at 2, 6, 12, and 24 h after cold stress ( Figure 8C; Table S5). This is consistent with the performance of its homolog CIR1 (AT5G37260) under freezing stress [48], indicating that FveMYB1R61 and CIR1 have high functional conservation. Rowan et al. found that two Arabidopsis MYB genes (AtMYB3 and AtMYB6) were upregulated by elevated temperatures, whereas another (AtMYB4) was downregulated [49]. In this study, these three AtMYB genes clustered together with five FveMYB genes (FveMYB13, -52, -66, -67, and -73) in clade C28, and the five FveMYB genes all exhibited downregulated expression under high temperature stress (Figures 3 and 8B; Table S5). This indicates that the genes in clade C28 may have undergone severe functional differentiation during the evolution process. This phenomenon also exists in clade Cr21. For example, MYBL2 (AT1G71030) is upregulated by high temperature, whereas FveMYB1R57 on the same branch is downregulated by high temperature (Figures 3 and 8D; Table S5).
To enhance their adaptability to various abiotic stresses, plants generally upregulate ABA. In plants, the MYB transcription factor superfamily is a primary regulator of ABAbiosynthesis genes [43]. Previous reports indicated that some MYB genes mediate ABA signals in response to abiotic stresses. For example, AtMYB96 is simultaneously induced by drought stress and ABA, and has been demonstrated to mediate ABA signaling in response to drought stress in Arabidopsis [50]. Overexpression of AtMYB15 enhances the drought and freezing tolerance of Arabidopsis and is responsive to ABA [51,52]. AtMYB2 is a transcriptional activator in the ABA signaling pathway, and it also responds to abiotic stresses such as drought, salt, and Pi starvation [13,53]. In this study, several genes exhibited similar or opposite expression patterns under cold or heat stress and ABA treatment. For example, the expression level of FveMYB1R47 increased 1.5-and 3.2-fold at 2 h under cold stress and ABA treatment, respectively, and then decreased (Figures 7C and 8C). Similarly, FveMYB1R20 was expressed in the control, but was scarcely expressed under heat stress and exogenous ABA treatment (Figures 7C and 8D). FveMYB14 was upregulated 9.2-and 10.3-fold at 2 h under heat stress and exogenous ABA treatment, respectively, and then gradually downregulated (Figures 7A and 8B). In addition, the expression of FveMYB59 peaked after 6 h of heat stress and then decreased. In the ABA treatment, FveMYB59 expression was lowest at 6 h and then increased (Figures 7A and 8B). We speculate that these genes may be associated with the crosstalk between the cold or heat stress response and the ABA signaling pathway.

Possible Roles of FveMYB Genes in Plant Organs and Their Responses to Plant Hormones
At present, the functions of some MYB genes in plant vegetative growth have been characterized. In Arabidopsis, slight drought stress can promote the expression of AtMYB60 to regulate root growth and increase water absorption capacity, whereas severe drought stress can inhibit the expression of AtMYB60, leading to stomatal closure and root growth inhibition [54]. AtMYB59 is abundantly expressed in roots and is involved with the cell cycle and root growth [55]. In this study, FveMYB38, a homolog of AtMYB59, had the highest expression in roots ( Figure 6A; Table S2), indicating that FveMYB38 may be associated with root growth in woodland strawberry. Interestingly, three FveMYB genes (FveMYB11, -63, and -97) that clustered with AtMYB60 in clade C37 were hardly expressed in roots but could be induced by ABA, and FveMYB63 and -97 were transiently upregulated by ABA, similar to AtMYB60 in Arabidopsis ( Figure 7A; Table S4). Therefore, we speculate that FveMYB63 and -97 may have lower expression levels in roots under normal conditions but can be induced under abiotic stress conditions to promote root growth. In addition, five FveMYB genes (FveMYB101, -112, -44, and -8; FveMYB1R81) were specifically expressed in roots and thus, may also participate in root formation and development, but this needs further study ( Figure 6A,C; Table S2).
Stolon propagation is the main asexual propagation method for cultivating strawberries. The stolon is a unique organ of cultivated strawberries and has not been found in model plants such as Arabidopsis, rice, and tomato. According to previous reports, GA can regulate the initiation and elongation of stolons. The formation of stolons is promoted by exogenous GA and inhibited by exogenous GA inhibitors [44]. This study found that five FveMYB genes (FveMYB1, -27, -53, -69, and -102) were primarily expressed in the stolon or stolon tip ( Figure 6A; Table S2), and four of these genes (FveMYB1, -27, -53, and -69) were transiently downregulated by GA ( Figure 7B; Table S4). We assume that these genes are involved in the GA signaling pathway to regulate stolon or stolon tip development. Wenjun Sun et al. reported that StMYB19 had the highest expression in potato stolons, and also had higher expressions in other seven organs [56]. In this study, three FveMYBs genes (FveMYB5, -25 and -14) were clustered with StMYB19 ( Figure S2). Among them, FveMYB5, -25 and -14 were not the highest expression in stolons, but they were expressed in various organs of strawberry ( Figure 6A), indicating that FveMYB5, -25 and -14 play a role in various organs of strawberry including stolons. FveMYB1, -27, -69, and -102 had no orthologous genes in potato, while FveMYB53 was clustered with potato StMYB143 ( Figure S2). FveMYB53 was expressed in most organs, with the highest expression in strawberry stolons ( Figure 6A). StMYB143 was also expressed in all potato tissues, but did not have the highest expression in potato stolons [56]. Therefore, the MYB genes, which regulate the growth and development of stolons, may have functional differentiation in strawberry and potato. In addition, the FPKM values of FveMYB62 and FveMYB77 were only >1 in the stolon tip and were almost undetectable in the other organs ( Figure 6A; Table S2). These genes may play a role in stolon elongation.
The roles of many MYB genes in the development of floral organs have also been reported. For example, AtMYB21 and AtMYB24 are flower-specific transcription factors [57,58]. Clustered together with AtMYB21 and AtMYB24 in clade C13, FveMYB99 was abundantly transcribed in open flowers but almost unexpressed in the other seven organs (Figures 3 and 6A; Table S2). In detail, FveMYB99 was abundantly expressed in mature anthers (Anther_12), mature carpels (Carpel_12), and unpollinated styles (0 d after pollination) (Table S6). This indicates that FveMYB99 is functionally conserved in both strawberry and Arabidopsis. AtMYB80 (formerly AtMYB103) can regulate the development of pollen and is functionally conserved in various crops such as wheat, rice, canola, and cotton [59][60][61]. As a homolog of AtMYB80, transcripts of FveMYB41 increased significantly in Anther_9, then decreased rapidly, and were almost nonexistent in Anther_12 (Table  S6). This indicates that FveMYB41 may be primarily involved in pollen development at the Anther_9 stage of woodland strawberry. AtMYB17 is involved in early inflorescence development via the regulation of the expression of LFY and AGL15 [62]. In this study, FveMYB106, a homolog of AtMYB17, was expressed abundantly in young carpels (Carpel_7-8) and young anthers (Anther_7-8). It was then gradually downregulated during the development of these floral organs (Figure 3; Table S6). This indicates that FveMYB106 is also involved in the development of premature floral organs. In addition, AtMYB33 and AtMYB65 are involved in the development of anthers [63], and their homolog, FveMYB35, was only expressed in open flowers, indicating that FveMYB35 may be involved in the formation of floral organs and primarily regulate the development of anthers. Other FveMYB genes (such as FveMYB1R96 and FveMYB1R69) expressed in flowers may also play roles in the development of woodland strawberry floral organs, although their orthologous genes in Arabidopsis have not been studied.

MYB Transcription Factors That Color Fruit
At present, many MYB genes have been reported to regulate the biosynthesis of anthocyanins in fruits, such as MrMYB1 of Chinese bayberry [64]; MdMYB10 and MdMYB1 of apple [63,65]; LcMYB1 of lychee [66]; MYB7 of peach [32]; MYB114, PbMYB10b, and PbMYB9 of pear [67,68]; SlMYB75 and SlMYBATV of tomato [69,70], etc. In Arabidopsis, AtMYB75, -90, -113, and -114 have also been reported to regulate anthocyanin biosynthesis [71,72]. In this study, two FveMYB genes (FveMYB10 and FveMYB75) clustered together with these four Arabidopsis genes in clade C21 ( Figure 3). Interestingly, from the immature to ripening stage of strawberry fruit, only FveMYB10 was upregulated, whereas FveMYB75 was not detected ( Figure 6B; Table S3). This indicates that FveMYB10 may regulate anthocyanin biosynthesis during the ripening of woodland strawberry fruit, whereas FveMYB75 may have a different function. This function of FveMYB10 has been reported by Chen and Castillejo C et al. [73,74], corroborating our results. During the fruit ripening of octoploid strawberry (Fragaria × ananassa), FaMYB10 promotes the biosynthesis of anthocyanins through the independent regulation of light and ABA [75]. Similarly, the expression of FveMYB10 was increased after ABA treatment ( Figure 6A; Table S4), but this result was only obtained by ABA treatment on woodland strawberry seedlings, ABA effect on FveMYB10 expression in woodland strawberry fruits needs further studies. In addition, transcription data from various organs showed that FveMYB76 was mainly transcribed in fruits and was almost unexpressed in other organs ( Figure 6A; Table S2). Further analysis revealed that FveMYB76 was only expressed in the ghost and wall tissues of the achene, and it was gradually upregulated with development of the achene (Table S6). Similarly, FveMYB100 was abundantly expressed in fruits, but hardly expressed or detected in other organs ( Figure 6A; Table S2). Detailed fruit tissue data showed that FveMYB100 was also mainly expressed in achene tissues (embryo, ghost, and wall), especially in Ghost_4 (Table S6). However, it is puzzling that FveMYB76 and FveMYB100 were expressed at very low levels in the immature and ripening fruit of the two woodland strawberry varieties ( Figure 6B; Table S3); therefore, further research is needed. Five MYB genes (FveMYB4, -22, -50, and -66; FveMYB1R57) stood out because they were highly expressed in early developing fruit but were significantly downregulated from the immature to the ripening stage ( Figure 6B,D; Table S3), indicating their roles in the early development of woodland strawberry fruit. These FveMYB genes may regulate the development and ripening of woodland strawberry fruit.

Conclusions
In this study, we identified 217 MYB proteins that were divided into four subgroups in woodland strawberry. The phylogenetic relationship of each MYB subgroup was consistent with the intron/exon structure and conserved motif composition. The RNA-seq analysis showed that most of the FveMYB genes expressed in multiple organs, only 11 FveMYB genes showed organ-specific expression. Eighty-four genes showed decreased expression, and eighteen genes showed increased expression during fruit ripening. In addition, several genes could be regulated by GA3, ABA, cold and heat treatments. This study provides some candidate MYB genes for the study of stolon formation, fruit development and ripening, and abiotic stress.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/horticulturae7050097/s1, Figure S1: Multiple alignment of 10 3R-MYB protein sequences in Arabidopsis and woodland strawberry, Figure S2: The phylogenetic tree of 2R-MYB proteins from woodland strawberry and potato, Table S1: Rename of 217 MYB genes in woodland strawberry,