Genome-wide identification and expression analysis of the R2R3-MYB gene family in tobacco (Nicotiana tabacum L.)

The R2R3-MYB transcription factor is one of the largest gene families in plants and involved in the regulation of plant development, hormone signal transduction, biotic and abiotic stresses. Tobacco is one of the most important model plants. Therefore, it will be of great significance to investigate the R2R3-MYB gene family and their expression patterns under abiotic stress and senescence in tobacco. A total of 174 R2R3-MYB genes were identified from tobacco (Nicotiana tabacum L.) genome and were divided into 24 subgroups based on phylogenetic analysis. Gene structure (exon/intron) and protein motifs were especially conserved among the NtR2R3-MYB genes, especially members within the same subgroup. The NtR2R3-MYB genes were distributed on 24 tobacco chromosomes. Analysis of gene duplication events obtained 3 pairs of tandem duplication genes and 62 pairs of segmental duplication genes, suggesting that segmental duplications is the major pattern for R2R3-MYB gene family expansion in tobacco. Cis-regulatory elements of the NtR2R3-MYB promoters were involved in cellular development, phytohormones, environmental stress and photoresponsive. Expression profile analysis showed that NtR2R3-MYB genes were widely expressed in different maturity tobacco leaves, and however, the expression patterns of different members appeared to be diverse. The qRT-PCR analysis of 15 NtR2R3-MYBs confirmed their differential expression under different abiotic stresses (cold, salt and drought), and notably, NtMYB46 was significantly up-regulated under three treatments. In summary, a genome-wide identification, evolutionary and expression analysis of R2R3-MYB gene family in tobacco were conducted. Our results provided a solid foundation for further biological functional study of NtR2R3-MYB genes in tobacco.

third helices can be folded to form a helix-turn-helix (HTH) structure, which is involved in DNA binding [2]. Three regularly spaced tryptophan residues (or other hydrophobic residues) form a hydrophobic core in the three-dimensional structure center of HTH, which is important for maintaining the configuration of HTH [3]. Sometimes, the tryptophan residues are replaced by other amino acids, such as aromatic amino acids and hydrophobic amino acids, especially in R3 domains. The C-terminal of MYB transcription factor usually contains a transcriptional activation region rich in acidic amino acids, which is responsible for multiple protein regulatory activities. Based on the number and types of MYB repeats, the MYB family is subdivided into four major groups, namely 4R-MYB, 3R-MYB/R1R2R3-MYB, 1R-MYB/MYB-related, 2R-MYB/R2R3-MYB [4].
R2R3-MYB transcription factor is predominantly present in plants, which contains R2 and R3 domains. It was reported that the R2R3-MYB genes probably evolved from the loss of the R1 repeat in the R1R2R3-MYB gene or from the duplication of the R1 repeat in 1R-MYB gene [5,6]. R2R3-MYB gene family is widely involved in the regulation of various biological processes, including plant growth and development, hormone signaling, primary and secondary metabolism [7][8][9]. For example, AtMYB106, AtMYB16 and AtMYB17 are involved in the regulation of trichome branching, petal epidermal cell morphogenesis and early inflorescence development, respectively [10][11][12]. The biosynthesis of flavonoids is directly or indirectly controlled by three genes (P, C1, Pl) encoding R2R3-MYB domain in maize [1,13,14]. MYB7 gene in Actinidia deliciosa positively regulated the biosynthesis of carotenoids and chlorophyll [15]. Knockdown the expression of MYB305 of the ornamental tobacco which contains a conserved R2R3-MYB DNA binding domain resulted in the decrease expression of related genes in nectarins and flavonoid biosynthetic [16]. In addition, most R2R3-MYB genes have also been suggested to regulate plant responses to biotic and abiotic stress conditions. For instance, overexpression of AtMYB75 in Arabidopsis can increase secondary metabolites (anthocyanins and flavonols), which can protect against pests [17]. OsMYB6 gene of rice as a stress-responsive factor which plays the role as a positive regulator in response to drought and salt stress resistance [18]. Increasing the expression of NtMYB4a may promote anthocyanin accumulation, and thereby increase antioxidant capability and tolerance to low temperature of tobacco plant [19]. Tobacco NtMYB12 overexpression adapts to low Pi stress environment via regulating the contents of flavonol and phosphorus [20].
Tobacco is one of the most important model plants [21]. Like other plants, tobacco is often subjected to stress such as low temperature, salt, strong light, drought and other stresses, which leads to the reduction of production. Leaf senescence is an important production process with a positive and orderly process accompanied by the changes of leaf color, cell structure, biochemical metabolism and gene expression level, along with a series of degradations, which is related to secondary metabolites, such as flavonoids, carotenoids, chlorophyll, etc. [22,23]. The regulatory role of R2R3-MYB genes in abiotic stress response, plant growth and development has been assessed in several studies [12,18,24]. The investigation of the R2R3-MYB gene family and their expression patterns under various stresses including abiotic stress and senescence in tobacco is of great significance for the study of plant physiology and development.
In this study, a comprehensive investigation of R2R3-MYB gene family, including gene structures, chromosomal localization, phylogenetic relationship, motif composition, duplication events and cis-element compositions was performed using the current tobacco genome sequence data. Moreover, the gene expression profile in different senescence stages of tobacco leaves and expression pattern under various abiotic stresses such as cold, salt, and drought treatments were analyzed. The objectives of this study were to systematically analyze the sequence structures of tobacco R2R3-MYB gene family and explore the evolutionary relationship of R2R3-MYB gene family in plant, and thereby, reveal the expression regulation of the R2R3-MYB gene family members under various stresses or adversity condition of senescence. The information derived from this study lay a foundation for further functional investigation on the R2R3-MYB gene family in tobacco.

Characterization and distribution of R2R3-MYB genes in tobacco genome
In this study, a total of 174 NtR2R3-MYB genes were identified in tobacco and were renamed from NtMYB1 to NtMYB174 (Table 1). The information of these NtR2R3-MYB genes and their corresponding proteins are showed in Table 1 and Additional file 1: Table S1, namely, gene ID, location, number of exons, protein length (aa), molecular weight (MW), theoretical isoelectric point (pI) and subcellular location. The protein lengths varied greatly from 192aa (NtMYB37, NtMYB39) to 505aa (NtMYB96). The molecular weights ranged from 22,003.55 Da (NtMYB37) to 55,710.67 Da (NtMYB96), and the theoretical isoelectric point (pI) ranged from 4.80 (NtMYB56) to 9.77 (NtMYB3). Subcellular localization prediction indicated that the majority NtR2R3-MYB members were located in the nucleus. In addition, the chromosome positions of A total of 107 NtR2R3-MYB genes were unevenly distributed on 24 chromosomes of tobacco, while 67 genes were mapped to unattributed scaffolds (Fig. 1). Chromosome 4 contained the biggest number of NtR2R3-MYBs (13 genes), while chromosome 17 contained 10 NtR2R3-MYBs, and chromosome 22 had 8 NtR2R3-MYBs. In contrast, chromosome 11 only contained 1 NtR2R3-MYB gene. In our study, 3 pairs of tandem duplication genes on chromosome 4 (NtMYB72/112, NtMYB72/138, NtMYB112/138) and 62 pairs of segmental duplication genes were identified in tobacco R2R3-MYB gene family (Fig. 2, Additional file 2: Table S2).

Phylogenetics and gene structure of the NtR2R3-MYBs
To investigate the evolutionary relationships among NtR2R3-MYB genes, a phylogenetic tree was constructed based on the amino acid sequences of 174 NtR2R3-MYB genes by using the MEGA-X software (Fig. 3A). The NtR2R3-MYBs were classified into 24 subgroups (A to X) with at least 50% bootstrap supported on phylogenetic trees (Fig. 3A). However, 5 NtR2R3-MYBs (NtMYB139, NtMYB9, NtMYB61, NtMYB103 and NtMYB27) could not be assigned to any of the 24 subgroups due to the low bootstrap values (< 50%). Among these subgroups, subgroup A (27 members) and X (21 members) were the two largest groups and these two subgroups represented more than 27% of the total NtR2R3-MYB members. In contrast, subfamilies B, C, F, J, K, N and T only contained two members.
Gene structure (Fig. 3B) analysis of NtR2R3-MYBs showed that the number of introns was varied from 0 to 11. The highest number (11) of introns was possessed by NtMYB173 and NtMYB174. Notably, a great number of R2R3-MYB genes (119 members, 68.39%) had a conserved gene structure with two introns and three exons, while 10 NtR2R3-MYBs completely lacked the introns. Similar exon-intron structural patterns were found among members within the same subgroup, especially the exon number and exon length were relatively conservative (Fig. 3).

Domain and motif analysis of the NtR2R3-MYBs
A total of 20 conserved motifs were predicted for 174 NtR2R3-MYB proteins using the online MEME program (Fig. 4). The lengths and conserved sequence of each motif is listed in Additional file 3: Table S3. The motif composition and distribution were found relatively conservative among members within the same subgroup (Fig. 4). The motif1, motif2 and motif3 located in the N-terminal of the majority NtR2R3-MYB protein sequences. There were conserved tryptophan residues in these three motifs, which were related to R2 and R3 domains. The type and number of motifs were similar in same subgroup, which suggested that motif pattern might be related to the function of MYB protein. Different subgroups usually possessed specific motifs, most of which located in the C-terminal. For example: motif7 and motif9 were specific to P; motif10 and motif13 only appear in A; motif8, motif14, and motif18 was presented in G, L, and Q alone, respectively.
To further explore the conservative domain of the NtR2R3-MYB proteins, the multiple alignment of the 174 NtR2R3-MYB protein sequences was performed based on DNAMAN software, and the R2 and R3 sequence logos of MYB were generated by WebLogo (Fig. 5). As a result, NtR2R3-MYB members possess the typical characteristics of MYB conserved domains, and the R2 and R3 repeat of NtR2R3-MYB contain about 52 amino acid residues. The R2 repeat contains three highly conserved tryptophan residues (W), forming a hydrophobic core zin HTH structure. The first tryptophan residues (W) of R3 repeat often was replaced by a phenylalanine (F), isoleucine (I) or leucine (L) residues, whereas the second and third tryptophan residues were highly conserved, and this result was consistent with A. thaliana [4]. We also observed that some amino acid residues showed highly conservative, such as G-2, E-8,  repeat. These highly conserved amino acid residues may be associated with conserved tryptophan residues to maintain the helix-turn-helix (HTH) structure of MYB transcription factor.

Promoter cis-elements analysis of NtR2R3-MYBs
Promoter cis-elements play critical roles in the initiation of gene expression. A total of 37 cis-regulatory were identified in the promoter region of NtR2R3-MYB genes, which could be classified into four categories elements, including cellular development, phytohormones, environmental stress and photoresponsive elements (Fig. 6, Additional file 4: Table S4). There were seven cis-acting elements related to cell development, such as CAT-box, MSA-like, GCN4_motif, CCAAT-box, MBSI, HD-Zip 1, and RY-element. Twelve phytohormone-responsive elements were identified, namely, CGTCA-motif, TGACGmotif, ABRE, P-box, TGA-element, TCA-element, AuxRR-core, TATC-box, GARE-motif, AuxRE, A-box and O2-site. These cis-elements are involved in JA/ MeJA, abscisic acid, gibberellin, auxin, salicylic acid responsiveness and zein metabolism regulation. Meanwhile, the ABRE responsiveness elements were the most common in the NtR2R3-MYB gene promoters. In addition, ten light responsive elements were calculated, including GT1-motif, G-Box, Box 4, MRE, ATC-motif, Sp1, ATCT-motif, ACE, 3-AF1 binding site, and AAACmotif. Almost all NtR2R3-MYB genes contained at least one phytohormone-responsive element in their promoter regions. There were 8 cis-regulatory elements that associated with response to external or environmental stresses were also present. This category includes lowtemperature responsive element (LTR), anaerobic induction elements (ARE, GC-motif ), drought-inducibility element (MBS), defense and stress responsive element (TC-rich repeats), circadian control element (circadian), wound-responsive element (WUN-motif ), as well as ATrich element. G-Box, Box 4, and ARE elements appear in most promoters of NtR2R3-MYB genes. The expression of these genes might be regulated by phytohormones, diverse light-responsiveness cis-elements, defense signaling transduction, and abiotic stresses during tobacco growth. To investigate the phylogenetic relationships of the R2R3-MYB gene family, phylogenetic tree was generated based on the 174 tobacco R2R3-MYB protein sequences and 126 Arabidopsis R2R3-MYB protein sequences by using MEGA-X with maximum likelihood (ML) method (Fig. 7). According to the bootstrap values (> 50%) of the phylogenetic tree, the R2R3-MYB family were clustered into 38 subfamily. Among them, the R2R3-MYB members of tobacco were distributed in 34 subgroups (named N1-N34). There was no NtR2R3-MYB member distributed on the subfamily of S3, S6, S12 and S15, while three subfamilies (N10, N20, N21) only contain the R2R3-MYB members from tobacco. In addition, the NtR2R3-MYB members were mainly distributed in N9 (9), N11(11), N22(13), N25(11), N27(9) and N34 (10), and the number of R2R3-MYB members of tobacco were about twice to triple than that in Arabidopsis in these subfamilies. These results indicated that there were some common ancestors of R2R3-MYB genes between tobacco and Arabidopsis, and specific expansion and divergence also occurred after their separation during the evolution process.

Expression changes of NtR2R3-MYB genes in five senescence stages of tobacco leaves
To analyze the expression pattern of NtR2R3-MYBs in tobacco leaves, the FPKM values of NtR2R3-MYB genes at five senescence stages of tobacco leaves were obtained from the transcriptome data (Additional file 5: Table S5), and NtR2R3-MYB genes with no expression or low expression level (FPKM < 0.5) were excluded. Finally, the expression profiles of NtR2R3-MYB genes of 78 NtR2R3-MYB genes were generated (Fig. 8). The results showed that the members of NtR2R3-MYB genes exhibited differential expression in tobacco leaves at different senescence stages (Fig. 8) and these 78 NtR2R3-MYB genes were classified into three groups (Fig. 8 I to III). A total of 9 NtR2R3-MYB members were including in group I, and these genes had high expression level at the M5 stage. In contrast, these genes showed relative low expression level at other four senescence stages (M1-M4). In group II, most of NtR2R3-MYB genes showed high expression level at the M1 stage and decreased regularly with the increase of maturity. In terms of group III, the expression level of NtR2R3-MYB genes showed increase first and then decreased with the increasing of the senescence degrees. The results indicated that there may be  Table S4 functional diversity among NtR2R3-MYB members during tobacco growth and development.

Expression of NtR2R3-MYB genes in response to abiotic stress
To analyze the expression pattern of NtR2R3-MYBs in response to abiotic stress, gene expression was investigated by using the Genevestigator tools based on transcriptome data. NtR2R3-MYB genes with no expression or low expression level (FPKM < 0.5) were excluded (Additional file 6: Table S6). Finally, the expression profiles of 69 NtR2R3-MYB genes were generated. The result showed that many genes showed significant up-regulated or down-regulated compared with the control group under cold and salt stress conditions (Fig. 9), and these genes including NtMYB34, NtMYB38, NtMYB42, NtMYB44, NtMYB46, NtMYB63, NtMYB67, NtMYB73, NtMYB79, NtMYB82 and NtMYB104 were clustered together with S1 and S7 subfamilies of Arabidopsis. In addition, it has been reported that the members of Arabidopsis R2R3-MYB gene family in S1 and S7 subgroups were related to various stress responses [25]. It is known that homologous genes in the same subgroup may have similar biological functions. To further explore the possible function of the tobacco R2R3-MYB genes, 15 tobacco R2R3-MYB genes that clustered with Arabidopsis R2R3-MYB gene members in S1 and S7 subgroups of the phylogenetic tree (Fig. 7) were selected for qRT-PCR analysis under cold, drought and salt stresses (Fig. 10). Compared with the control, the expression of four NtR2R3-MYB genes (NtMYB36, NtMYB45, NtMYB46 and NtMYB110) showed significantly up-regulated. Interestingly, the expression of NtMYB46 showed significantly up-regulated in response to all the stresses. In addition, the expression patterns of 15 genes in Fig. 9 and Fig. 10 were not completely consistent, and this phenomenon may be due to trial or sampling differences. The result Fig. 9 Expression of 69 NtR2R3-MYB genes in response to cold and slat treatments. RNA-seq data was obtained by using Genevestigator software. Expression values from RNA-seq data were log10-transformed implied the functional dissimilation among the tobacco R2R3-MYB genes.

Discussion
R2R3-MYB gene family members are widely distributed in eukaryotes [26]. With the development of genome sequencing, the whole-genome analysis of R2R3-MYB gene family has been identified in numerous species, including 126 in Arabidopsis thaliana, 89 in watermelon, 110 in rice, 157 in maize, 244 in soybean, and so on [7,14,[27][28][29]. In this study, 174 R2R3-MYB genes of tobacco were identified, the number was greater than those identified in Arabidopsis, maize, rice, watermelon, but less than that in soybean. As an allotetraploid, the genome size of Nicotiana tabacum is 4.5Gb, while that of Arabidopsis, rice, maize, and soybean is 125 Mb, 430 Mb, 2300 Mb, and 1.025Gb, respectively [30][31][32][33][34]. In this case, it seems that there is no direct correlation between the number of R2R3-MYB genes and genome size in these plants. It was reported that polyploidization and gene region-specific duplication (tandem duplication and segmental duplications) were important mechanisms for the generation and expansion of gene families in plant [35]. In our study, phylogenetic analysis found that the majority subfamilies contained the R2R3-MYB members both from tobacco and Arabidopsis, and however, a few subfamilies possessed only the R2R3-MYB members either from tobacco or Arabidopsis, suggesting that they might be derived from a common ancestor, and moreover, the R2R3-MYB gene family could also be undergone species specific differentiation after their separation. A total of 3 pairs of tandem duplication genes and 62 pairs of segmental duplication genes were identified in tobacco R2R3-MYB gene family, implying that the segmental duplication events were the main source for the expansion of R2R3-MYB gene family in tobacco, and this result was possible due to the allotetraploid of tobacco. The evolution of gene family largely depends on the organization of gene structure. The varied length of nucleotide sequence among 174 NtR2R3-MYBs indicated the complexity in the Nicotiana tabacum L. genome. The molecular weight and isoelectric point values of NtR2R3-MYB proteins were also different among family members, suggesting their functional divergence. In addition, NtR2R3-MYB proteins contained 20 conserved motifs with different compositions, and their members were clustered in the same subfamily containing similar type and number of motifs, demonstrating the conservation and diversification of R2R3-MYB gene family of tobacco. It has been reported that the ancestor MYB gene has no intron, but the intron insertion event occurred in the MYB domain region under very occasional condition, and this intron pattern was kept conserved during the long evolution process [5,36]. In our study, the majority NtR2R3-MYB genes possess typically splicing pattern of three exons and two introns, which exists in the conserved R2 and R3 repeats, and this result was in consistent with the previous reports in other plants [14,37]. In addition, the NtR2R3-MYB proteins are comprised of the highly conserved MYB domain (R2 repeat and R3 repeat), R2 repeat contained a conserved LRPD motif at the C-terminal and three highly conserved tryptophan residues (W), whereas R3 repeat exist diversity at the first tryptophan residues (W), which could be substituted by other residues such as phenylalanine (F), isoleucine (I) and leucine (L). Meanwhile, substitution of amino acid occurred frequently in some sites of the MYB domain, and these regions may play important roles for the evolution and functional differentiation of tobacco R2R3-MYB protein. Similar results have been reported in other species, including Arabidopsis, Zea mays (Fig. 5) [14,25].
The cis-elements of promoter play a key role in initiating gene expression. Genes with different cis-regulatory elements in the promoter sequences of genes may result in different expression patterns in pepper [38]. In our study, a total of 37 cis-elements related sequences were identified in the promoter region of NtR2R3-MYB genes, and among them, 7 were related to cell development, 12 engaged in phytohormone-responsive, 10 for light responsive elements, and 8 involved in stresses cisregulatory elements, suggesting the different function of regulatory elements of NtR2R3-MYB genes. These highly diverse cis-regulatory elements in the promoter region of NtR2R3-MYB genes may also reflect the functional divergence at the transcriptional level.
R2R3-MYB genes have been proved to be related to biotic and abiotic stress [8]. For example, TaMYB344overexpressing of wheat enhanced the tolerance of transgenic tobacco to drought, heat and high salt stress [39]. OsMYB2 was induced by cold, salt, and dehydration stress in rice [40]. In present study, the profiling data of gene expression was used to dissect the functional roles of NtR2R3-MYB genes, and different expression patterns were identified among NtR2R3-MYB genes in response to various abiotic stresses (Fig. 9). This result indicated the functional differentiation of NtR2R3-MYB genes. In addition, the expression of 15 NtR2R3-MYB genes were analyzed by qRT-PCR in response to three abiotic stress conditions, including cold, salinity, and drought. A total of 13 NtR2R3-MYB genes could be significantly regulated by at least two treatments except NtMYB44 and NtMYB104, implied that NtR2R3-MYBs may be involved in the cross-talk of different signaling pathways under stress. Generally, genes with similar structure will be clustered in the same subfamily and these genes may have similar biological functions. It was reported that two R2R3-MYB genes of Arabidopsis (AtMYB60 and AtMYB94) clustered in subfamily S1 (Fig. 7) were involved in the physiological regulation under salt and drought treatments [41,42], and therefore, the orthologous clustered in S1 may have similar function. In this study, two NtR2R3-MYB genes (NtMYB38 and NtMYB46) which clustered with these two R2R3-MYB genes of Arabidopsis (AtMYB60 and AtMYB94) in S1 subfamily showed similar response patterns under salt and drought stresses, suggesting that NtMYB38 and NtMYB46 may be involved in the response under salt and drought stresses, and further experiments showed be conducted to validation these functions.
Tobacco leaf senescence is often regarded as a kind of adversity, and R2R3-MYB genes were proved to be involved in the senescence accompany with secondary metabolites, including flavonoids, carotenoids, chlorophyll, and so on [23,43]. In this study, diverse expression patterns were found among the NtR2R3-MYB genes inferring that the functional differentiation of family gene members should also be coexisting (Fig. 8). For example, NtMYB3/4 and NtMYB3/6 appeared to belong to segmental duplication genes pairs (Fig. 2, Additional file 2: Table S2), and our result showed that NtMYB4 and NtMYB6 expressed in all the five senescence stages, with especial high expression level in M5 stage, whereas NtMYB3 had no expression in all investigated senescence stages (M1-M5) (Fig. 8). It has been reported that there are three outcomes in the evolution of duplicate genes theoretically, including nonfunctionalization, neofunctionalization and subfunctionalization [44]. Therefore, it was inferred that NtMYB3 degenerated and lost its original function during the long evolution process. Notably, the expression level of majority members in group II including such as NtMYB146, NtMYB108, NtMYB90, NtMYB159, NtMYB120, NtMYB17 (Fig. 8) were decreased precisely corresponding consistent to the increasing of the senescence degrees, and these gene may closely relate to leaf senescence and can be further developed as the measure of leaf senescence or maturity. Whether or not the differential expression of NtR2R3-MYB genes leads to the changes of secondary metabolites such as flavonoids, carotenoids or chlorophyll still needs to be investigated. In short, our results provide using information for their further functional exploration.

Conclusions
In this study, a total of 174 R2R3-MYB genes were identified in tobacco (Nicotiana tabacum L.) genome and these genes were divided into 24 subfamilies. The NtR2R3-MYB genes were distributed randomly on 24 tobacco chromosomes. A total of 3 pairs of NtR2R3-MYB genes were founded to be originated from tandem duplication and 62 pairs NtR2R3-MYB genes were originated from segmental duplication. Cis-regulatory elements of the NtR2R3-MYB promoters were involved in cellular development, phytohormones, environmental stress and photoresponsive. The members of NtR2R3-MYB genes showed differential expression pattern in different maturity tobacco leaves, and differential response were also found under different abiotic stresses (cold, salt and drought) for 15 NtR2R3-MYBs. Our results provided valuable information for further functional study of NtR2R3-MYB genes in tobacco.

Chromosome localization and gene duplication
The physical position and chromosomal distribution information of NtR2R3-MYB genes were obtained by using the MapInspect software (http:// mapin spect. softw are. infor mer. com/) [53]. The possible segmental duplication and tandem duplication events were defined based on the method reported by Wang et al. (2010) [54]. Both chromosomal localization and duplication events of the NtR2R3-MYB genes were graphic displayed using the TBtools software [55].

Multi-sequence alignment and phylogenetic classification
To explore the evolutionary relationship of R2R3-MYB gene family, the full protein sequences of R2R3-MYB from Arabidopsis and tobacco were used for the phylogenetic tree construction. Multiple sequence alignment was performed using ClustalW program in MEGA-X software [56], and the phylogenetic tree was constructed using the maximum likelihood (ML) with 1000 bootstraps.

NtR2R3-MYB genes expression analysis
The tobacco variety of Cuibi 1 (CB-1) was used in this study. To analyze the NtR2R3-MYB genes expression profile at different senescence stages, five senescence stages (M1, M2, M3, M4 and M5) of middle leaves (8th to 10th) judged by the appearance characteristics were collected for tests [57]. The FPKM (Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) value of the NtR2R3-MYB genes at five senescence stages of tobacco leaves were extracted from our recent RNA-Seq data [58]. The expression profile of NtR2R3-MYB genes at different senescence stages were measured by their FPKM value, and the heat map was generated using the Heatmap function of R gplots package [59]. Additionally, expression analysis of NtR2R3-MYBs under salt (SRP193166) [60] and cold stress (SRP097876) [61] of tobacco were performed using Genevestigator software. Genes with low expression level (FPKM < 0.5) were filtered.

Plant treatments and quantitative real-time PCR analysis
To further decipher the expression pattern of NtR2R3-MYB genes in response to various abiotic, tobacco seeds were sown in sterilized mixed soil (vermiculite: humus = 1:1) under the condition of 22 °C and 16 h light/8 h dark photoperiod for 60 days [62]. The plantlets of 60 days were transplanted into a tray with a nutrient solution for 3 days in growth chamber, and then were exposed to the abiotic treatments, including cold (4 °C), drought (10% polyethylene glycol) and salt (200 mM NaCl), respectively. Untreated plantlets were used as control (CK). The samples for gene expression analysis were collected 6 h after treatment, and three biological replicates per treatment and 3 leaves for each sample from different plantlet were gathered and these samples were immediately stored at − 80 °C prior to RNA extraction.
Total RNA was extracted using the Hipure Plant RNA Mini Kit (Magen Biotech, Shanghai, China) and the cDNA was synthesized using the SMART kit (Takara) according to the manufacturer's protocol. The qRT-PCR primers of NtR2R3-MYB genes were designed by online