WOX genes were identified from the genome data representing the four Rosaceae species (P. bretschneideri, P. persica, F. vesca, and P. mume, respectively). Two different methods were used to identify WOX genes in the P. bretschneideri, P. persica, F. vesca, and P. mume genome: (1) BLASTP search using Arabidopsis and rice WOX protein sequences according to previous research methods (Cao et al., 2016a,d), and (2) the screening of Hidden Markov Model profile (PF00046) in four Rosaceae genome using DNAtools software with an e-value cut off of 0.001 (Gehring, 1992). All candidate WOX proteins were confirmed to have a complete WOX domain using both Pfam (Punta et al., 2011), SMART dadabases (Letunic et al., 2012) and InterProScan database (Zdobnov and Apweiler, 2001).

The multiple alignment of WOX proteins in five plant species (P. bretschneideri, P. persica, F. vesca, P. mume, and A. thaliana) was performed using CLUSTAL_X software (Thompson et al., 1997). Subsequently, we constructed NJ (neighbor-joining) tree using MEGA version 5.1 software (Tamura et al., 2011) with the following parameters: bootstrap (1000 replicate), pairwise deletion and Poisson correction. At the same time, we used ML (maximum-likelihood) and ME (Minimum-evolution) methods to generate the phylogenetic trees to validate the topologies.

The online program Gene Structure Display Server (Hu et al., 2014) was used to detect the exon–intron structure of cDNAs and genomic DNA sequences. Subsequently, the MEME (Multiple Em for Motif Elicitation, Version 4.11.1) program (Bailey et al., 2015) was used to obtain the motifs in all candidate WOX proteins, with the parameters: the maximum number of motifs at 20, and the optimum motif width between six and 200 residues. Furthermore, the Pfam database (Punta et al., 2011), SMART software (Letunic et al., 2012), and InterProScan database (Zdobnov and Apweiler, 2001) were used to annotate these structure motifs.

According to the comparisons of the specific regions containing WOX genes, we carried out microsynteny analysis across the four Rosaceae species. Similarly, the WOX genes of P. bretschneideri, P. persica, P. mume, and F. vesca were categorized based on their classification in the evolutionary tree. Subsequently, all WOX genes in P. bretschneideri, P. persica, F. vesca, and P. mume were set as anchor sites, according to their physical location. Then the flanking protein-coding genes of the WOX gene in one species were compared with those in other species. The criterion for dividing an interspecific synteny block is to locate three or more conserved homologous genes within 100 KB between genomes (BLASTP E-value < 10^-10) (Wang et al., 2012).

To future understand whether the WOX genes have undergone positive selection during evolution, maximum likelihood codon models (site models and branch-site models) in PAML software (Yang, 2007) were performed. Three pairs of models (M0 vs. M3, M1a vs. M2a, and M7 vs. M8) were utilized to detect positive selection sites. In the free site models, M0 (one ratio), M1a (neutral), M2a (selection), M3 (discrete), M7 (beta), and M8 (beta and ω) were evaluated by the likelihood ratio test (LRT). The LRT was used to judge which model was more suitable in the two models, and the amino acids sites with positive selection were obtained by the Bayesian method of PAML software (Yang, 2007).

Functional divergence analysis of amino acid sequence data was performed using Diverge 2.0 combined with constructed phylogenetic tree (Gu, 1999, 2006; Gu et al., 2013). The type I functional divergence led to a change of functional limitation, which was highly correlated with the evolution rate after gene duplication (Gu, 1999, 2006; Gu et al., 2013). Type II functional divergence did not result in a change in the functional limitation of the members after gene duplication, but the change of physical and chemical properties of amino acid residues (Gu, 1999, 2006; Gu et al., 2013).

To identify putative cis-elements in promoter regions of WOX genes, the PlantCARE database (Lescot et al., 2002) was used. 2000 bp genomic sequence upstream of the start codon (ATG) was used for cis-acting elements analysis.

The normalized data (Fragments Per Kilobase Exon model per Million mapped fragments, FPKM) during F. vesca development was reported by Darwish et al. (2013), and available from SGR GBrowse. A gene was thought to be expressed if the FPKM value was greater than or equal to 0 FPKM in at least one of the 14 tissues. Subsequently, the transcriptome data of FvWOXs was visualized using the R software¹.

For identification of WOXs gene families, the genome data of tested Rosaceae species was subjected to HMM and BLASTP searches. As a result it is revealed, presence of 9, 10, 10, and 14 WOX genes in P. mume, P. bretschneideri, P. persica, and F. vesca, respectively (Supplementary Table S1). These WOX genes were named according to method of Haecker et al. (1991). For this purpose, the phylogenetic tree was carried out based on multiple sequence alignments for the full-length WOX protein sequences of tested Rosaceae species and A. thaliana. The identified WOX genes of these tested four Rosaceae species were renamed according to the evolutionary relationship as shown in Supplementary Table S1. Subsequently, the distribution of these WOX genes on chromosomes was identified, based on genomic annotation information. As shown in Figure 1, it is discovered that the WOX genes were unevenly distributed among the chromosomes in each of tested four species. In the P. bretschneideri genome, four of WOX genes were distributed on chromosome 15, while remaining distributed on chromosomes 3, 5, 6, 12, and 17. In F. vesca, five and four of 14 WOX genes were distributed among chromosomes 3 and 5, respectively. In both P. persica and P. mume, three WOX genes distributed were found on one chromosome (no. 7 and no. 8, respectively), with others scattered across different chromosomes (Figure 1).

To investigate the possible evolutionary history of the WOX genes in the tested Rosaceae species, we carried out a joint phylogenetic analysis using three methods; ME, ML, and NJ. Based on previous report that WOX13 subfamily was an ancient member in the WOX gene family (Deveaux et al., 2008), the WOX13 subfamily was selected as an outgroup to root phylogenetic tree. All the tree topologies generated by the three methods (ME, ML, and NJ) were largely consistent with each other, with only minor changes in internal branches (Figure 2 and Supplementary Figures S1, S2). Therefore, only NJ phylogenetic tree was used in the following analysis. Previous studies on WOX genes have confirmed that motifs FYWFQNH, FYWFQNR, and YNWFQNR were representative markers for the WUS/modern clade, intermediate clade and ancient clade, respectively (Graaff et al., 2009; Ge et al., 2016). Confining the previous results (Deveaux et al., 2008; Xin et al., 2010), our evolutionary analysis exposed a total of 58 members of the WOX genes in P. bretschneideri, P. persica, F. vesca, P. mume along with A. thaliana. These 58 members of WOX genes were divided into three clades and nine subfamilies. The Modern clade contained a total of six subfamilies (WUS, WOX1, WOX2, WOX3, WOX4, and WOX5), and Intermediate clade included two subfamilies (WOX9 and WOX11), while Ancient clade just had a WOX13 subfamily, which were consistent with the evolutionary relationships of WOXs in other species (Xin et al., 2010; Hedman et al., 2013; Nardmann and Werr, 2013; Lian et al., 2014). Remarkably, we found that all subfamilies contained at least one WOX member from each of the four Rosaceae species (Figure 2). These results imply that rapid duplication of WOX genes occurred before these dicotyledonous species were diverged.

Phylogenetic analysis revealed that three pairs of paralogous genes were found among the WOX genes, which were consistent with the previous notion that most members of the WOX gene family are represented by pairs of orthologous genes. As shown in Figure 3, up to eight pairs of orthologous WOX genes were shared by P. persica and P. mume, whereas only two pairs of orthologous WOX genes found between P. bretschneideri and P. persica. However, no orthologous WOX genes were found between F. vesca and other species. These results were consistent with the evolutionary relationships among these four Rosaceae species (Dickinson et al., 2007; Cao et al., 2016b).

Previous studies have shown that gene structural diversity is an important resource for the evolution of multigene families (Liu et al., 2009; Cao et al., 2016c). To understand the structural diversity of the WOX genes in Rosaceae, gene structures of PbWOXs, PpWOXs, PmWOXs, and FvWOXs were deduced. It is revealed that these WOX genes contained different numbers of exons as shown in Figure 4A. For example, FvWOX11 only contained one exon, while FvWOX9A contained the largest number of exons (5). Moreover, 16, 15, and 3 of WOX genes contained two, three and four exons, respectively. These results suggested that the functional diversity of WOX genes may be in consequence due to exon loss or gain during the evolution of the WOX gene family. Subsequently, gene structures of the WOX paralogous and orthologous gene pairs were further analyzed. Among these genes, we found that the exon number of seven gene pairs had changed, including FvWOX9A/FvWOX9B, FvWOX13B/FvWOX13C, PmWOX3/PpWOX3, PbWOX4/PpWOX4, PmWUS/PpWUS, PmWOX11/PpWOX11, and PmWOX13B/PpWOX13B. By comparing among these seven gene pairs, it was found that one exon was lost in FvWOX9B, PmWOX3, PpWOX4, PpWUS, PmWOX11 and PmWOX13B, while one exon was obtained in FvWOX9A, PpWOX3, PbWOX4, PmWOXWUS, PpWOX11, and PpWOX13B. It may happen during the long evolutionary period. Previous studies have proposed that introns could be specifically inserted and remained in the plant genome during evolution (Rogozin et al., 2003; Carmel et al., 2007; Cao et al., 2016d). In our study, these phenomenon were observed, which might explain the functional differences and diversity of closely related WOX genes, such as PbWOX3 and FvWOX3, PbWOX13B, and FvWOX13A (Figure 4A).

Furthermore, it was observed that 20 of the conserved motifs were found in the 43 WOX proteins using MEME website (Supplementary Table S2). These motifs were annotated by using Pfam and SMART. Motif 1, present in all subfamilies, was identified to encode for a conserved homodomain. In addition to the homodomain, most of the WOX members within the same clade shared the similar motif compositions as shown in Figure 4B. These results reinforced the classification of WOX subfamilies. However, several motifs were unique to the proteins in some clades. For example, Motif 4 was unique to Ancient clade (clade I: WOX13 subfamily), Motif 5 to WUS clade (clades A–F) and Motif 15 to clade C (Supplementary Table S2 and Figure 4B). To some extent, these specific motifs may play an important role in the clade or subfamily, as well as contribution to the functional divergence of WOX genes.

Based on their amino acid sequences, the newly identified WOX gene family members were found to contain the conserved homeodomain by multiple sequence alignment of Cluxa2.0 with default parameters. The conserved homeodomain was selected for the visualized results by ESPript 3 (Gouet et al., 2005). The homeodomain structures of these four species were highly similar with each other. They contained a helix-loop-helix-turn-helix structure with either 65 or 66 amino acid residues. A total of 11 conserved sites (Q, L and Y in helix1; I, V, W, F, N, K, and R in helix3) of homeodomain reported previously (Gehring, 1992; Xin et al., 2010), were also conservative in the WOX proteins of the Rosaceae species (Figure 5). These findings suggested that these amino acid residues could play an important role in their functions. In addition to the previously reported conserved amino acid sites, other conserved amino acid sites have been identified in this study, such as P, L, and I in helix 2, Q and F in helix 3, as well as G in the turn region. Interestingly, an extra Y residue was found in the homeodomain of PbWUS, PmWUS, PpWUS, FvWUS and AtWUS, compared with other members of WOX gene family in P. bretschneideri, P. persica, P. mume and F. vesca and A. thaliana. Similar finding have been reported by Mayer et al. (1998) and Xin et al. (2010), that the homeodomains of A. thaliana WUS, O. sativa WUS, Z. mays WUS1, Z. mays WUS2, S. bicolor WUS and P. trichocarpa WUS were composed of 66 amino acid residues containing an extra Y residue by multiple sequence alignment, which indicates that this residue might play an important role on the function of WUS TF. Remarkably, in Arabidopsis, AtWOX5 (without Y residue between Helix1 and Loop) could replace AtWUS (containing Y residue) to maintain the dynamic balance of stem cells in the shoot apical meristem (Sarkar et al., 2007)

It was reported previously that the WUS protein contains three functional domains, including WUS-box, acidic region and EAR-like motif (Xin et al., 2010; Xiaoxu et al., 2016). These functional domains significantly contribute to its function as a TF (Xin et al., 2010; Xiaoxu et al., 2016). In present study, Motif 5 (WUS box: amino acids, TLLFP) was observed to be in the presence of all WOX proteins in WUX clade (clades A–F) (Supplementary Table S2 and Figure 4B). In clade C, the Motif 15 (EAR-like motif: amino acids, SLELSL) was found in all WOX proteins. However, no acidic region was identified in all WOX proteins (Supplementary Table S2 and Figure 4B). These results were consistent with previous findings that acidic region may be an important function domain only in Arabidopsis WUS gene (Xin et al., 2010; Xiaoxu et al., 2016).

Microsynteny has been surveyed in different species to understand the position of the homologous genes (orthology or paralogy) (Cannon et al., 2003; Yan et al., 2004; Cao et al., 2016a). In this study, microsynteny analysis was carried out for identification of homologous relationships within the WOX genes in P. bretschneideri, P. persica, F. vesca, and P. mume (Figure 3). Additionally, to measure the linkages and molecular history among WOX genes, a stepwise gene-by-gene reciprocal comparison was performed. In general, if the flanking genes in the chromosome region of the target gene contained three or more pairs of genes that are collinear, they could be considered as the conserved microsynteny (Lin et al., 2014; Cao et al., 2016a).

Primarily, the intraspecies microsynteny was investigated among four Rosaceae species. However, it was revealed that no collinear WOX genes were observed (Figure 3). These results suggested that independent duplication events were the main expansion pattern of WOX gene family members. Consequently, we analyzed the relationship of the WOX genes within each interspecies. The results exposed that the nine clades containing 38 WOX genes were found, among which 10 were from P. persica and F. vesca, 9 from P. bretschneideri and P. mume, respectively. Then several higher levels of microsynteny found in subfamilies A–G. Among these microsynteny some were remarkably inverted, duplicated such as PpWOX3/PmWOX3, PbWOX1/PpWOX1, and FvWOX13C/PpWOX13B (Figure 3). Usually, genome segments in the same group may evolve from a single sequence, which led to species differentiation (Tripoli et al., 2005; Jing et al., 2016). However, sequence fragments from the same group are considered to be homologous genes, and their genetic evolution resulted in species segregation (Tripoli et al., 2005; Jing et al., 2016). Remarkably, with the construction of the phylogenetic tree, the conservation of microsynteny in different families gradually emerged. Furthermore, some flanking genes were not conserved in each microsyntenic group. Therefore, it was speculated that these new genes were later than this duplication event. Interestingly, several lower levels of microsynteny were also found, such as PpWOX13B/FvWOX13A and PbWOX13A/PmWOX13B in clade I, PbWOX11/PpWOX11 and PpWOX11/PmWOX11 in clade H (Figure 3). These results strongly suggested that the ancient large-scale duplications could follow by gene rearrangement and loss.

To investigate whether WOX genes have undergone strong selection pressures in the evolution of WOX gene family, we used site and branch-site models in the CODEML program of PAML software to detect positive selection sites (Yang, 2007). However, no positive selection was detected among these genes (Supplementary Table S3). These finding imply that relaxed purifying selection might play a major role in the evolution of WOX genes. Nardmann and Werr (2013) have shown that the WOX gene expansion was resulted from the increased complexity of plant morphology (Nardmann and Werr, 2013). These results proposed that the novel members after gene expansion were retained with partly overlapping expression domains and functions. These results are similar with findings of Nardmann and Werr (2013) who reported that with the relaxed purifying, dosage effects will lead to a selective advantage (Nardmann and Werr, 2013).

For further investigation to comprehend significant differences in selection pressures among WUS clade, intermediate clade and ancient clade, the branch-site models were performed using PAML software. It was exposed that no significant positive selection observed in the different branches of WOX gene family (Supplementary Table S4). This result was in contrast to the previous report that some significant positive selection sites were fixed in WOX genes of peanut (Wang et al., 2015). Because significant positive selection usually exerts its effects only in few sites and in a short period of evolutionary processes, it is difficult to detect positive selection. Thus, the selected signal could be diluted by the purifying selection (Zhang, 2005). However, as the WOX coding regions are highly conserved among members of orthologous families, the absence of strong positive selection was expected in Rosaceae species.

Due to the fact that significant positive selection could only detect a limited number of adaptive selection events, we performed a functional divergence analysis according to method used by Cao et al. (2016a). The DIVERGE software was used to calculate functional divergence of type I or II between gene clades in WOX genes with posterior analysis. In general, type I functional divergence usually resulted in a specific amino acid selectivity change, i.e., evolutionary rate change. The type II functional divergence only led to the change of physical and chemical properties of amino acids, which were occurred after gene duplication. In present study, to avoid the emergence of false positives, the sites with a posterior probability Q K > 0.9 were set as the key amino acid sites arising the functional differences according to previous experimentation reported (Yin et al., 2013; Cao et al., 2016a). Our results showed that five key amino acid sites (144, 154, 161, 165, and 166) were identified as type I functional divergence between Ancient and Modern (Table 1), while just one key site (152) between Intermediate and Modern (Table 1). The chi-square test (x 2) found that the P-values of Ancient/Modern and Intermediate/Modern were less than 0.05, reaching a significant level. Interestingly, among these three clades, no specific type II functional divergence site (Q K > 0.9) was detected (Supplementary Table S5), suggesting that the physicochemical properties of amino acid sequences between these Rosaceae WOX genes were highly identical.

Two thousand bp sequences of upstream from start codon (ATG) among the putative WOX genes, were used for analysis of WOX promoters by searching, against the PlantCARE website. Consequently, we detected various types of cis-acting elements in the promoter region of 43 WOX genes (Supplementary Table S6). These results indicated that the same type of WOX might carry out different functions. MBS and ABRE elements were found to be distributed in promoter region of most WOX genes, implying that WOX genes were transcriptionally regulated upon salt stress and dehydration. Remarkably, we found that the cis-elements exhibit significant differences in the promoter regions of duplicated WOX genes. These results indicated that the duplicated WOX genes may exhibit different regulation features.

To explore the role of the WOX gene family in F. vesca development process, the expression of the FvWOX genes was explored. The results showed that their expression levels were divergent from each other, indicating that they may be functionally active among all tissues except FvWOX5B (Supplementary Figure S3), which was located in Pollen with no expression. At the same time, most of FvWOX genes exhibited developmental stage-specificity, such as higher expression of FvWOX13A, FvWOX3, FvWOX13B, and FvWOX1 in flowering, and FvWOX4, FvWOX5, and FvWOX5A in embryo (Supplementary Figure S3). Surprisingly, we found that FvWOX13A, FvWOX9A, and FvWOX1 were highly expressed among all tissues, indicating that these genes were persistent and very important during development process of F. vesca.