Identification of the non-redundant MYB genes in the pear genome¹ was performed using DNATOOLS software. First, the consensus protein sequences (PF00249) of the MYB hidden Markov model (HMM) were downloaded from Pfam². Then, this HMM profile was used as a query to identify all MYB-containing sequences in pear by searching against the pear genome with an E-value of <1e⁻³. Subsequently, all candidate PbMYBs were verified using Pfam and SMART³ to confirm that they contained the core domains. Based on the sequence alignments generated by ClustalX software, all potentially redundant MYB sequences were discarded.

Analysis of exons and introns was carried out using Gene Structure Display Server 2.0⁴ (Guo et al., 2007) by comparing the coding sequences (CDS) with their corresponding gene sequences. The isoelectric point (PI) and protein molecular weight (kDa) of each PbMYB were obtained using ExPASy proteomics server⁵. MapInspect ⁶ was used to determine the locations of the PbMYB genes on the pear chromosomes. The MYB proteins were analyzed using MEME (Multiple Expectation Maximization for Motif Elicitation ⁷; Bailey and Elkan, 1995) to confirm the presence of the conserved motifs with the following parameters: an optimum motif width of no less than 6 and no more than 200 and a maximum number of motifs of 20. The conserved motifs were annotated by Pfam⁸ and SMART³.

The sequences of the R2 and R3 domains of 105 PbR2R3-MYBs were aligned using ClustalW⁹ to analyze their sequence features. Multiple alignment files for these domains were submitted to WebLogo¹⁰ (Crooks et al., 2004) using the default settings to acquire sequence logos.

To analyze the promoter regions of the PbMYB genes, 2,000 bp regions upstream of these genes were examined based on the positions of the genes provided in pear GigaDB database¹. PLACE database¹¹ (Higo et al., 1999) was employed to examine putative cis-acting regulatory DNA elements in the promoter regions of these genes.

A total of 133 MYB protein sequences in A. thaliana were obtained from TAIR¹². In addition, 18 protein sequences of well-known plant R2R3-MYB genes were collected from GenBank¹³. Based on the findings of previous studies (Riechmann et al., 2000; Jiang et al., 2004b), AtCDC5 is an R2R3-MYB; thus, we determined whether its orthologs exist in the pear genome in our phylogenetic analysis. We also included the 4RMYB and 3R-MYB proteins in this analysis. A phylogenetic dendrogram was constructed from the ClustalW-aligned MYB proteins by the neighbor-joining (NJ) method using MEGA5.2 software¹⁴ (Tamura et al., 2011) with the default parameter values.

In this study, to identify collinearity, the whole-genome sequences of pear were downloaded to our local server, and then Multiple Collinearity Scan toolkit (MCscanx; Wang et al., 2012) was used to detect the microsynteny relationships between each pair of chromosomes. The resulting microsynteny chains were evaluated using ColinearScan software with an E-value of <1e⁻¹⁰.

The PbMYBs protein sequences were aligned to the NCBI non-redundant (nr) protein database by BlastP using Blast2GO software¹⁵ (Conesa et al., 2005) with the default parameters. After obtaining a GO annotation for each PbMYB, we plotted the GO classifications using WEGO online tool¹⁶ (Ye et al., 2006).

The position of each PbMYB was marked using a Perl script. Orthologs of the MYB genes in pear, strawberry (Fragaria vesca), yang mei (Prunus mume), and peach (Prunus persica) were identified using OrthoMCL¹⁷ (Li et al., 2003). These MYB orthologs were confirmed by performing a reciprocal BLAST search. Then, we used Circos¹⁸ (Krzywinski et al., 2009) to plot the relationships of the orthologous pairs among the four species.

DnaSP v5.0 software (Librado and Rozas, 2009) was used to calculate Ks and Ka substitution rates. Then, the selection modes of PbMYB paralogs were determined by analyzing the Ka/Ks ratios. Additionally, sliding window analysis was performed to determine the Ka/Ks ratios of all of the encoding sites of the PbMYB paralogs.

To examine PbMYB gene expression, pear (P. bretschneideri cv. Dangshan Su) fruit samples were collected on 19 April (15 days), 14 May (39 days), 21 May (47 days), 30 May (55 days), 6 June (63 days), 22 June (79 days), 15 July (102 days), and 29 August (145 days) in 2015 after flowering (DAF). Three or more fruits were collected at each stage from 40-years-old pear trees grown on a farm in Dangshan, Anhui, China. Total RNA was extracted from the collected samples using Trizol reagent (Invitrogen). Then, the DNase-treated RNA was reverse transcribed using M-MLV reverse transcriptase (Invitrogen). Primers (Supplementary Table S1) were designed for real-time quantitative PCR (qRT-PCR) using Primer Express 3.0 software (Applied Biosystems). The tubulin gene (Wu et al., 2012) was used as an internal reference, and the primers for this gene were synthesized by Sangon Biotech, Co., Ltd. (Shanghai, country-region China). qRT-PCR was performed using an ABI7500 instrument to examine the gene expression in cDNA samples from the cross-pollinated varieties at different developmental stages. Each reaction was performed in triplicate. The relative expression levels of the PbMYB genes were calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). The reaction mixtures contained the following in a total volume of 20 μl: 10 μl of SYBR Premix Ex Taq II (2x), 1 μl of template cDNA, 0.5 μl of forward and reverse primers, and water. PCR amplification was carried out under the following conditions: 50°C for 2 min and 95°C for 30 s, followed by 40 cycles of 95°C for 15 s, 60°C for 20 s, and 72°C for 20 s.

The full-length CDS of PbMYB25 and PbMYB52 were cloned from pear and constructed into pCambia1304 and pCambia1301 vectors (Clontech, Beijing, country-region China), which both contain a CaMV 35S promoter and GFP gene, resulting in fusion PbMYB-GFP genes. The primers used for gene cloning and vector construction are shown in Supplementary Table S1. The constructed pPbMYB-GFP vectors were electroporated into Agrobacterium tumefaciens EHA105 using a Gene Pulser Xcell (BIO-RAD, country-region USA). The suspensions were infiltrated into the leaves of Nicotiana tabacum using the injection method. The expressed PbMYB-GFP was observed by confocal laser scanning microscopy (CarlZeiss LSM710, Germany).

To identify the MYB family genes in the pear genome, the MYB HMM profile (Pfam: 00249) was used as query in a BlastP search against the pear genome database. This search resulted in the identification of 138 candidate PbMYB proteins in pear, and their sequences were examined to further verify whether they contained MYB-DNA binding domains using Pfam database⁸. Four PbMYB candidates were discarded because they lacked these binding domains. Notably, PbMYB38 (Chr15: 9477433–9479211) and PbMYB39 (Chr15: 9296590–9298368) encode proteins that are the same length and have the same PI and molecular weight (Supplementary Table S2). Further analysis showed that PbMYB80 (17385133–17386109) and PbMYB83 (17513104–17514080) on chromosome 10, PbMYB103 (8444013–8447010), PbMYB104 (8385864–8388861) and PbMYB105 (7592807–7596377) on chromosome 3, and PbMYB131 (24626034–24629526) and PbMYB133 (19686471–19689983) on chromosome 17 were also identical. Thus, the five overlapping genes were removed. As a result, a total of 129 non-redundant MYB genes, including 22 R1-MYBs, 105 R2R3-MYBs, and two R1R2R3-MYBs, some of which might be pseudogenes, were identified in pear based on the presence of one, two or three MYB-DNA binding domain repeats (Supplementary Table S2).

We performed phylogenetic reconstruction of all of the PbMYB genes using the NJ method (Figure 1). A phylogenetic tree was created, in which the PbMYB genes were categorized into 31 major groups (C1-C31) with supported bootstrap values. The results of subsequent analysis of their exon–intron structures were consistent with those of phylogenetic analysis. Most genes clustered in the same group exhibited similar exon–intron structures, particularly with regard to the number of introns, such as C1, C2, C3, etc. (Figure 1). However, a few exceptional cases were observed. For example, the members of C19 contained different numbers of introns. In addition, their intron lengths were variable, ranging from several 10s of 1000s to approximately 25,000 nucleotides. These results demonstrate the presence of highly conserved structures within the PbMYB subfamily and high sequence diversity among the different PbMYB groups.

A composite evolutionary tree was generated to evaluate the evolutionary relationships of the MYB genes between pear and Arabidopsis. These members could be divided into 41 groups with bootstrap support of > 50%. Each subfamily is outlined in a different color in Figure 2. Most of the PbMYB genes were clustered with their counterparts in Arabidopsis with high bootstrap support in the phylogenetic tree (Figure 2). In previous studies (Stracke et al., 2001; Dubos et al., 2010), 126 AtMYBs from Arabidopsis were grouped into 25 subfamilies, and the defined clades are labeled in our composite evolutionary tree (Figure 2). Most large subfamilies in our study (C1, C2, C3, C26, etc.) are supported by the previous studies with high bootstrap support. Nevertheless, some small subfamilies (C23, C24, C15 etc.) were not retrieved from the NJ trees (Stracke et al., 2001; Dubos et al., 2010). For example, the subfamily C24 only included PbMYBs and not AtMYBs, indicating that some changes occurred among MYBs of different species during the evolutionary process.

To demonstrate conservation at particular positions, WebLogo was used to generate sequence logos. As shown in Figure 3, the conserved amino acids among the PbMYB domains are very similar to those of Z. mays, A. thaliana, Cucumis sativus, Vitis vinifera, Populus trichocarpa, Glycine max, and O. sativa (Du et al., 2012b). These results indicate that many conserved amino acids are present in the R2 and R3 repeats of PbMYBs, particularly the characteristic Trp (Figures 3A,B). Three conserved Trp residues were identified in the R2 repeat. However, only the second and third Trp were conserved in the R3 repeat, and the first Trp was commonly substituted with Phe or Ile. Subsequently, 20 conserved motifs were identified in the pear MYB proteins using MEME web server. The MYB DNA-binding domains were represented by motifs 1, 2, 4, and 8. As shown in Supplementary Figure S1 and Supplementary Table S3, the majority of the PbMYB proteins contained several of motifs 1, 2, 4, and 8. For example, subfamily C31 only contained motif 8, while subfamilies C28 and C29 contained motifs 1 and 4. Most of the closely related genes had the same motif compositions, indicating that there are functional similarities between MYB proteins within the same subfamily. Furthermore, we detected some subfamily specific motifs, which might be required for subfamily specific functions, such as motif 20 in subfamily C21 and motif 8 in C31. In addition, some motifs, such as motifs 3 and 6, were present in nearly every subfamily, reflecting their importance to the functions of PbMYB proteins.

Based on the starting positions of the MYB genes within the chromosomes, 122 PbMYB genes were found to be unevenly distributed among 16 pear chromosomes (Supplementary Figure S2). Another seven genes (PbMYB26, PbMYB49, PbMYB51, PbMYB66, PbMYB102, PbMYB109, and PbMYB116) could not be conclusively mapped to any chromosome. Chromosomes 9 and 14 contained the largest numbers of PbMYB genes (14), followed by chromosome 15 (11), and the lowest numbers of these genes were found on chromosomes 7 and 15 (2). High densities of PbMYBs were observed in several specific regions, including the tops of chromosomes 10 and 15, the tops and bottoms of chromosomes 6 and 9, and the bottoms of chromosomes 8, 11, 14, and 17; in contrast, PbMYBs were not detected in the middle parts of the chromosomes.

We also investigated segmental duplication events. A total of 70 collinear gene pairs (Figure 4; Supplementary Table S4) were found in the pear genome, which might have resulted from whole-genome duplication in pear millions of years ago (Li et al., 2014). Among them, PbMYB26 on scaffold 1015.0 was collinear with PbMYB24, and PbMYB51 on scaffold 570.0 was collinear with PbMYB55, etc. The collinearity relationships of the MYB genes in the pear chromosomes are presented in Figure 4 and Supplementary Table S4. Furthermore, two duplicated pairs of genes, PbMYB50/PbMYB59 and PbMYB113/PbMYB118, were located on the same chromosome, and the other 68 pairs of genes were located on different chromosomes or scaffolds. Their locations within the chromosomes indicate that tandem (2) or segmental duplication (68) has occurred in these regions (Supplementary Table S4).

To explore the selective constrains among the duplicated MYB genes, we calculated Ks, Ka and the Ka/Ks ratio for the 70 duplicated gene pairs (Supplementary Figure S3; Supplementary Table S4). Generally, a Ka/Ks ratio of <1 indicates that the gene pairs are under negative or purifying selection. Further, a Ka/Ks of >1 indicates positive selection, and a Ka/Ks of =1 suggests neutral evolution. In this study, the Ka/Ks ratios for all PbMYB gene pairs were <1 (Supplementary Figure S3; Supplementary Table S4), suggesting that the MYB gene pairs have mainly evolved under purifying selection in pear, with limited functional divergence after the duplication events. Compared with the other replication pairs, the Ka/Ks ratios of the three pairs PbMYB6/PbMYB11, PbMYB22/PbMYB131, and PbMYB23/PbMYB87 were larger with values of 0.82, 0.90, and 0.95, respectively. Sliding window analysis was performed to evaluate the Ka/Ks ratios of the CDS at different sites. The results showed that the Ka/Ks ratios of some coding sites were greater than 1, indicating strong positive selection (Supplementary Figure S4; Supplementary Table S4). As shown in Supplementary Figure S4, the R2 domain has undergone purifying selection, whereas the R3 domain has been mainly affected by positive selection. These results are consistent with the results of analysis of the PbMYB domain characteristics.

The promoters of the predicted PbMYB genes were analyzed, including a 2000 bp region upstream of the transcription start site (ATG). Five regulatory elements, including ABRE (Narusaka et al., 2003; Yamaguchi-Shinozaki and Shinozaki, 2005), LTRE (Yamaguchi-Shinozaki and Shinozaki, 2005), DRE (Narusaka et al., 2003), BOXP (Logemann et al., 1995) and BOXL (Maeda et al., 2005), were detected by performing a search of these promoter sequences against PLACE database¹⁹. Surprisingly, 91% of the PbMYB genes were found to contain putative ABRE elements in their promoter regions (Supplementary Table S5), indicating that ABA can affect the expression levels of PbMYB genes. In addition, some PbMYB genes contained BOXP and BOXL components, suggesting that they may be involved in the regulation of lignin biosynthesis. A comparison of the distributions of the five regulatory elements (ABRE, LTRE, DRE, BOXP, and BOXL) in the promoter regions among the PbMYBs revealed significant differences in the cis-elements of the promoter sequences of the 70 duplicated MYB genes (Figure 4; Supplementary Table S4), implying that the duplicated genes may have different regulatory features.

The functions of the putative PbMYB proteins in pear were predicted by GO annotation analysis. Based on amino acid similarity, 129 PbMYB proteins were categorized into 23 functional categories (Supplementary Figure S5; Supplementary Table S6) of the three main ontologies, namely biological process, cellular component, and molecular function. Among these categories, cell, cell part, binding, biological regulation, cellular process, metabolic process, and pigmentation were predominant. Several genes were assigned to the “death” and “anatomical structure formation” categories. Analysis of the molecular function annotations revealed that most of the PbMYB proteins function in DNA binding, followed by chromatin binding. Further, analysis of the cellular component annotations revealed that these proteins are predominantly localized to the nucleus.

The orthologous MYB genes were identified to explore the evolutionary history of the MYB family in Rosaceae plants. We identified 128 orthologous gene pairs between pear and peach and 114 orthologs between pear and yang mei but only 91 orthologs between pear and strawberry (Figure 5; Supplementary Table S7). These findings likely reflect the closer relationship of pear with peach and yang mei (Li et al., 2014). Remarkably, some collinear gene pairs identified between pear and peach were not detected between pear and yang mei, such as PbMYB22/ppa011751m, PbMYB26/ppa007594m and PbMYB43/ppa018123m, indicating that these orthologous pairs formed after yang mei diverged from the common ancestor of pear and peach. Some collinear gene pairs identified between pear and peach and pear and yang mei were not identified between pear and strawberry. For example, ppa017136m and Pm013279 are orthologous to PbMYB16, while ppa007438m and Pm009892 are orthologous to PbMYB17, which indicates that these orthologous pairs formed after strawberry diverged from the common ancestor of pear, peach and yang mei. In addition, two or more MYB genes from strawberry, peach and yang mei matched one pear MYB gene (Figure 5). For example, Pm027045, Pm013279 and Pm005570 are orthologous to PbMYB77, ppa007438m, and ppa016708m are orthologous to PbMYB17, and mrna25149 and mrna03712 are orthologous to PbMYB6, suggesting that these genes are probably paralogous gene pairs.

Increasing evidence indicates that MYB TFs play key roles in the regulation of plant metabolism. For example, AtMYB32, AtMYB4, AmMYB330, EgMYB2, PttMYB21a, etc. are involved in the regulation of lignin synthesis in plants. To identify additional MYB genes that participate in the regulation of lignin synthesis in pear, we generated a composite evolutionary tree that included many MYBs in other plants that take part in the regulation of lignin synthesis and all PbMYBs using the NJ method (Supplementary Figure S6). According to the phylogenetic relationships, we identified 20 PbMYBs (chosen from lignin-related gene subgroups, Supplementary Figure S6) that might participate in lignin synthesis in pear fruit. The qRT-PCR results showed that these genes exhibited diverse expression patterns at 15, 39, 47, 55, 63, 79, 102, and 145 DAF (Figure 6). The expression of 12 genes (PbMYB18, PbMYB25, PbMYB29, PbMYB32, PbMYB49, PbMYB55, PbMYB57, PbMYB90, PbMYB93, PbMYB95, PbMYB103, PbMYB119, and PbMYB127) was significantly increased at 15 days after blooming. Therefore, they were inferred to be involved in lignin synthesis as upstream genes. Expression of 7 MYB genes (PbMYB12, PbMYB25, PbMYB51, PbMYB52, PbMYB80, PbMYB96, and PbMYB89) was significantly increased at 47, 55, or 63 DAF, with similar expression patterns as the key genes involved in regulation of the lignin synthesis pathway (Xie et al., 2013). These results suggest that these TFs might be specific activators of the lignin synthesis pathway. PbMYB94 expression was significantly increased at 102 days after blooming, indicating that this gene is involved in not only regulating the lignin pathway but also mediating the synthesis of other materials, such as flavonoids.

We further analyzed the expression of five pairs of segmental duplicated genes at eight time points after flowering. The results revealed that two of these genes had similar expression patterns at several time points, while others exhibited variable expression patterns. For example, the expression of PbMYB55 was significantly increased at 15 DAF, while that of PbMYB51 was significantly increased at 55 days after blooming.

Nuclear localization of TFs is critical for their regulatory functions. Many factors influence the process of nuclear entry, such as environment stress, the cell cycle, and the developmental stage. Although, 40–60 kD proteins can pass through nuclear pores by diffusion, the movement of a protein through a nuclear pore is a proactive process. The entry process requires one or more nuclear localization signals. To examine the subcellular localization of PbMYB25 and PbMYB52, pPbMYB-GFP expression vectors were constructed and transformed into N. tabacum. As shown in Figure 7, green fluorescence signals from the expressed fusion PbMYB25-GFP and PbMYB52-GFP genes were specifically distributed within the nuclei. In contrast, green fluorescence from the expressed GFP alone was observed throughout the whole cell, demonstrating a constitutive expression pattern.