Genome-wide identification of cassava R2R3 MYB family genes related to abscission zone separation after environmental-stress-induced abscission

Cassava plants (Manihot esculenta Crantz) resist environmental stresses by shedding leaves in leaf pulvinus abscission zones (AZs), thus leading to adaptation to new environmental conditions. Little is known about the roles of cassava R2R3 MYB factors in regulating AZ separation. Herein, 166 cassava R2R3 MYB genes were identified. Evolutionary analysis indicated that the 166 R2R3 MYB genes could be divided into 11 subfamilies. Transcriptome analysis indicated that 26 R2R3 MYB genes were expressed in AZs across six time points during both ethylene- and water-deficit stress-induced leaf abscission. Comparative expression profile analysis of similar SOTA (Self Organizing Tree Algorithm) clusters demonstrated that 10 R2R3 MYB genes had similar expression patterns at six time points in response to both treatments. GO (Gene Ontology) annotation confirmed that all 10 R2R3 MYB genes participated in the responses to stress and ethylene and auxin stimuli. Analysis of the putative 10 R2R3 MYB promoter regions showed that those genes primarily contained ethylene- and stress-related cis-elements. The expression profiles of the genes acting downstream of the selected MYBs were confirmed to be involved in cassava abscission zone separation. All these results indicated that R2R3 MYB plays an important regulatory role in AZ separation.

BLASTP searches of the complete genomes in Manihot esculenta (annotation v.6.1) were performed by using Arabidopsis MYB protein sequences as queries 3 . A total of 216 putative amino acid sequences in the cassava genome containing MYB or MYB-like repeats were identified form the cassava genome database. Subsequently, the redundant sequences of candidate MYBs were discarded from the data set according to the similarity of the sequences. A simple modular architecture research tool was used to analyze all the putative MYBs to confirm the MYB domain 3 . The remaining cassava MYBs possessing incomplete ORFs were also excluded from further analysis. All cassava MYB proteins were manually inspected to ensure that the putative genes contained 2 or 3 MYB repeats. In total, 166 genes were defined as R2R3 MYB proteins in the cassava genome (Additional file 1). The gene set represented approximately (166/33,033) 0.5025% of the annotated genes in the cassava genome (33,033 genes), a proportion greater than that of Arabidopsis genes (0.4596%).
To confirm that the R2 and R3 MYB repeats are highly conserved across all 166 R2R3 MYB proteins identified form the cassava genome, multiple alignment analysis was conducted to identify the homologous R2R3 MYB domain in all 166 MYB proteins by investigating the frequency of the most prevalent amino acids within the MYB domain 3 . In general, the basic regions of the R2 and R3 repeats both had 50 basic residues (Fig. 1). The R2 and R3 repeats of cassava were highly conserved in sequences; 10 out of 50 in R2 and 24 out of 50 amino-acids in R3 were 100% conserved in all MeMYB proteins (Fig. 1). As compared with those in other plant species, the R2 and R3 MYB repeats in cassava R2R3 MYB family contain characteristic amino acids. For example, highly conserved Trp residues were distributed in 166 cassava R2R3 MYB that have been demonstrated to play a role in the sequence-specific binding domain 3 (Fig. 1).

Phylogenetic reconstruction of the R2R3 MYB superfamily between cassava and Arabidopsis.
To study the phylogenetic relationships among the cassava R2R3 MYB super-family members in different plant species, all identified cassava R2R3 MYBs and all R2R3 MYBs from Arabidopsis were subjected to multiple sequence alignment 3 , and a phylogenetic tree was generated by using the R2R3 MYB amino acid sequences from cassava and Arabidopsis. The resulting phylogenetic tree contained 11 groups, termed groups A to K (Fig. 2). Most of the groups contained R2R3 MYB from both the cassava and the Arabidopsis genomes, and all these members from both Arabidopsis and cassava were grouped together, thus indicating that these genes may have the same function 3 . The C group had the most members, including 64 members from both cassava and the Arabidopsis genome. Specifically, 40 members were identified form cassava, and 24 were identified from Arabidopsis. The K group contained the least number of members. This group included only 1 member from the Arabidopsis genome. The D and J groups contained sister gene pairs identified in both Arabidopsis and cassava on the basis of the phylogenetic tree. The two sister gene pairs, AtMYB121 and Manes.02G058900 and AtMYB103 and Manes.02G017300 were grouped in the D group, whereas AtMYB91 and Manes.09G092900 were grouped in the J group.

Phylogenetic analysis and analysis of conserved gene structures and protein motifs of the R2R3 MYB gene family in cassavas.
To study the phylogenetic relationships among the cassava R2R3 MYB superfamily members, a phylogenetic tree was generated with 166 R2R3 MYB amino acid sequences from cassava (Fig. 3). On the basis of the clades with at least 50% bootstrap support, the 166 typical members of the cassava R2R3 MYB gene family were subdivided into 17 subgroups and designated as S1 to S17 (Fig. 3). Moreover, the tree topology after maximum likelihood (ML) analysis was essentially the same as that of the former unrooted phylogenetic tree 3 , thus indicating that these phylogenetic trees were in good agreement. The low bootstrap support for the internal nodes of these trees was consistent with phylogenetic analysis of MYBs in other organisms. A total of 62 sister pairs of putative paralogous genes were identified among the 166 R2R3 MYB genes, and 26 of them had high bootstrap support (≥ 98%).
In order to identify the conserved motifs in the MYB protein sequences, the Multiple Em for Motif Elicitation (MEME; version 4.9.0) program tool was used to analyze all the full-length protein sequences of 166 MeMYBs. Ten motifs were identified in the MeMYBs, and the motif lengths identified by MEME were between 8 and 50 amino acids (Fig. 3). The number of the conserved motifs in each MYB varied between 1 and 9. Most MeMYBs had motifs 1, 2, 3, and 6, and most members in the same subgroup shared more than one identical motif. The majority of the close members in the phylogenetic tree exhibited similar motif compositions, thus suggesting that the members may have the similar function in the same subgroup 3 .

Transcriptome identification of R2R3 MYB genes expressed in cassava leaf pulvinus-petiole AZs during both ethylene-and water-deficit stress-induced leaf abscissions. To analyze the R2R3
MYB gene expression profiles during ethylene and water-deficit stress-induced leaf abscission, cassava genome microarray (NimbleGen) analyses were performed. The results indicated that 41 R2R3 MYB genes were differentially expressed after ethylene-induced leaf abscission, and 38 R2R3 MYB genes were differentially expressed after water-deficit stress-induced leaf abscission (Additional table 2). In total, 26 R2R3 MYB genes were differentially expressed after both ethylene-and water-deficit stress-induced leaf abscissions. In contrast, 15 R2R3 MYB genes were differentially expressed exclusively after ethylene-induced leaf abscission, and 12 were differentially expressed after water-deficit stress-induced leaf abscission (Additional table 2). Of the R2R3 MYB genes induced by the water-deficit stress treatment, six SOTA clusters (WS1-WS6) were separated into four groups of primary expression patterns ( Fig. 4 and Additional table 3). The first group, cluster WS1, was up-regulated at the early and middle stages during abscission compared with the control, whereas the second group, clusters WS2 and WS3, was up-regulated at the later experimental time points. The third group, clusters WS4 and WS6, was down-regulated at the later experimental time points T4, T5 and T6. The fourth group, cluster WS5, was down-regulated at the early experimental time points ( Fig. 4 and Additional table 3).
The six ethylene-induced R2R3 MYB gene SOTA clusters (ES1-ES6) were divided into four main expression patterns ( Fig. 4 and additional table 4). The first group, cluster ES1, was up-regulated at the later experimental time points. The second group, cluster ES2, was up-regulated throughout the abscission process, with the highest levels of expression at the early and middle experimental points. The third group, clusters ES3 and ES4, was down-regulated later at T4, T5 and T6. The fourth group, clusters ES5 and ES6, was down-regulated at the early and middle experimental time points compared with the control.
Comparison of R2R3 MYB expression profiles between ethylene-and water-deficit stress-induced leaf abscissions indicated that R2R3 MYB subfamily genes are widely expressed in the cassava abscission zone. To identify the R2R3 MYB genes that participated in both ethylene-and water-deficit stress-induced leaf abscissions, the R2R3 MYB expression profiles were compared in both treatments by using SOTA clustering. Because the expression patterns of R2R3 MYB genes were nearly identical between ethylene-and water-deficit stress-induced leaf abscissions, the similar expression patterns at each time point were compared in both ethylene and water-deficit stress treatments.
The R2R3 MYB genes that were up-regulated during the early and middle experimental periods in response to both treatments (Fig. 4) were first examined. During water-deficit stress treatment-induced leaf abscission, 3 R2R3 MYB genes in the WS1 cluster exhibited increased expression during the early stages of water-deficit stress-induced leaf abscission, and GO annotations indicated that all 3 genes participate in the responses to an ethylene stimulus (GO:0009723) and an auxin stimulus (GO:0009733) (Fig. 5). During ethylene treatment-induced leaf abscission, 8 R2R3 MYB genes in the ES2 cluster exhibited increased expression during the early and middle stages of leaf abscission. GO annotation indicated that 5 of the genes participate in the response to stress (GO: 0009651), 4 genes participate in the response to an ethylene stimulus (GO:0009723), and 4 genes participate in  (166) and Arabidopsis (126) were aligned using ClustalW, and the phylogenetic tree was constructed with MEGA 5.0 using the neighbor-joining method based on the p-distance model with 1000 bootstrap replicates. Each subfamily is represented by a specific color. The phylogenetic tree was constructed with MEGA 5.0 using the neighbor-joining (NJ) method with 1,000 bootstrap replicates based on a multiple alignment of 166 amino acid sequences of R2R3 MYB genes from cassava. Bootstrap values greater than 50% are presented on the nodes. The 17 major subfamilies are indicated, with S1 to S17 marked with colorful backgrounds. Protein motif schematic diagram of the conserved motifs in the R2R3 MYB proteins in cassava, which were elucidated using MEME. Each motif is represented by a number in the colored box. The black lines represent the non-conserved sequences. the response to an auxin stimulus (GO:0009733) (Fig. 5). Two genes, Manes.15G040700 and Manes.05G177900, were expressed during the early and middle stages during leaf abscission in response to both ethylene and water-deficit stress, and both participate in the pathways in response to an ethylene stimulus (GO:0009723) and to an auxin stimulus (GO:0009733). Manes.15G040700 encodes MYB15, a member of the R2R3 factor gene family. Manes.15G040700 exhibits a high expression ratio (compared with the T1 time point) at the T4 time point in both ethylene and water-deficit stress treatments. Manes.05G177900 encodes MYB73, a member of the R2R3 factor gene family. A high expression ratio (compared with the T1 time point) appeared at the T4 time point in ethylene treatment and the T3 time point in water-deficit stress.
Later in leaf abscission (T4, T5 and T6), WS2, WS3 (water-deficit stress treatment) and ES1, ES2 (ethylene treatment) exhibited similar expression patterns (Fig. 4). In response to water-deficit treatment, 6 R2R3 MYB genes exhibited increased expression. GO annotation indicated that 4 participate in response to an abscisic acid stimulus (GO:0009737), and 3 participate in response to stress (GO:0009651) (Fig. 5). In response to ethylene treatment, 17 R2R3 MYB genes exhibited increased expression. GO annotation indicated that 10 of the genes participate in response to stress (GO:0009651), 8 genes participate in response to an ethylene stimulus (GO:0009723), and 5 of genes participate in response to an auxin stimulus (GO:0009733) (Fig. 5). Comparative analysis indicated that 8 genes, i.e., Manes.05G114400, Manes.01G118700, Manes.15G081900, Manes.02G194800, Manes.06G066600, Manes.05G012100, Manes.14G104200, and Manes.03G117500, were expressed during both ethylene-and water-deficit stress-induced leaf abscission. Manes.05G114400 encodes MYB62, a member of the R2R3 MYB transcription family. Manes.05G114400 exhibited a high expression ratio (compared with the T1 time point) at the T6 time point in both ethylene and water-deficit stress treatments. Manes.01G118700 encodes MYB43, and a high expression ratio (compared with the T1 time point) appeared at the T5 time point in water-deficit stress treatment. However, this gene was down-regulated in ethylene-induced abscission. Manes.15G081900 encodes MYB4, a member of the R2R3 MYB transcription family, which is involved in wounding and the osmotic stress response. This gene exhibited a high expression ratio (compared with the T1 time point) at the T6 time point in both ethylene and water-deficit stress treatments. Manes.02G194800 encodes MYB3, and a high expression ratio (compared with the T1 time point) appeared at the T6 time point in both ethylene and water-deficit stress treatments. Manes.06G066600, which encodes MYB62, a member of the R2R3 MYB transcription family involved in regulation of phosphate starvation responses, exhibited a high expression  ratio (compared with the T1 time point) at the T6 time point in both ethylene and water-deficit stress treatments. Manes.05G012100 encodes MYB78, a member of the R2R3 MYB transcription family. This gene exhibited a high expression ratio (compared with the T1 time point) at the T5 time point in ethylene treatment and the T6 time point in water-deficit stress treatment. Manes.4G104200 encodes MYB62, which had a high expression ratio (compared with the T1 time point) at the T6 time point in both ethylene and water-deficit stress treatments. Manes.03G117500 encodes MYB4, which exhibited a high expression ratio (compared with the T1 time point) at the T4 time point in both ethylene and water-deficit stress treatments.
Promoter motif prediction for R2R3 MYB subfamily genes expressed during both ethylene-and water-deficit stress-induced leaf abscission. As described above, 10 genes were expressed in response to both ethylene and water-deficit stress treatments. To further understand the potential functions of the genes in regulating abscission zone development, the promoters of the 10 R2R3 MYB subfamily genes were analyzed to identify cis-elements. For this analysis, 2-kbp sequences of putative promoter regions were examined for potential cis-regulatory elements that are responsive to water-deficit stress and ethylene. Three drought stress response cis-elements, S000415, S000414 and S000176, as well as two ethylene response cis-elements, S000037 and S000457, were frequently identified within the promoter regions of the genes. The results suggested that all the genes contained more than 20 drought stress elements and more than 5 ethylene response elements (Additional  table 5).
Quantitative real-time RT-PCR analysis of expression profiles for 10 selected R2R3 MYB genes in cassava leaf pulvinus-petiole AZs during both ethylene-and water-deficit stress-induced leaf abscission. To confirm the reliability of the 10 selected R2R3 MYB genes that were up-regulated in the transcriptome data, FDR-corrected P-values < 0.005 were used as the significance criterion in at least one of the time points in both ethylene-and water-deficit stress-induced leaf abscission, to evaluate the expression patterns in transcriptome. As shown in Table 1, most of the FDR-corrected P-values of 10 selected R2R3 MYB genes in the transcriptome that were expressed at six time points in both ethylene and water-deficit stress treatments were less than 0.005. Moreover, some of the FDR-corrected P-values in at least one of the time points in both treatments were less than 1e-6. Compared with the FDR-corrected P-values in water-deficit stress-induced abscission, more values of the FDR-corrected P-values less than 1e-6 appeared in ethylene stress-induced abscission. Together, the data suggested that the 10 selected R2R3 MYB genes that were up-regulated in the transcriptome in both ethylene-and water-deficit stress-induced leaf abscissions were reliable.
To further confirm the expression patterns of 10 R2R3 MYB subfamily genes in leaf abscissions, their expression levels were identified by quantitative real-time RT-PCR. The relative expression of 10 R2R3 MYB genes was significantly different in both ethylene-and water-deficit stress-induced leaf abscissions (Fig. 6) 3 . The study is the first to identify and characterize the R2R3 MYB genes in cassava. Herein, 166 cassava R2R3 MYB genes were classified into 11 subfamilies (Additional file 1 and Fig. 2). The number of cassava R2R3 MYB subfamily members (166) far exceeded that in the Arabidopsis genome (126), thus indicating that the abundance of R2R3 MYB genes in cassava have expanded potentially via gene duplications. Gene duplication is considered to be the primary driving force underlying new gene functions 25,31 . In this study, numerous R2R3 MYB genes were observed to participate in both ethylene-and water-deficit stress-induced abscissions, thus suggesting that the R2R3 MYB subfamily members have important functions in cassava plant development.

R2R3 MYB subfamily genes suggested the function of regulating the progression of cassava leaf abscission by participating in various biological pathways. In our previous work, we have
reported that ROS regulate the separation of abscission zones under water-deficit stress through transcriptomic, physiological, cellular, molecular, metabolic, and transgenic methods. Moreover, ROS are increased by the accumulation of proline and polyamine in cassava abscission zones under water-deficit stress 25 . The levels of proline, polyamine, ROS, ethylene and auxin have been proven to have the function of regulating the process of abscission zone development 25 . Here, molecular methods were used to identify the key R2R3 MYB genes that regulate cassava AZ separation in response to ethylene and water-deficit stress. To identify the R2R3 MYB genes expressed during cassava leaf abscission in response to both treatments, the R2R3 MYB gene expression profiles of pulvinus-petiole AZs were used for SOTA clustering. In total, 41 and 38 R2R3 MYB genes were expressed during ethylene-and water-deficit stress-induced leaf abscission, respectively (Additional table 2). Among the R2R3 MYB genes, 26 were expressed during both ethylene-and water-deficit stress-induced leaf abscission (Additional table 2), thus suggesting that similar genes regulate leaf abscission in response to these stresses. Comparative analysis of R2R3 MYB gene expression profiles during ethylene-and water-deficit stress-induced leaf abscission resulted in the selection of 9 R2R3 MYB genes distributed among the different time points, thus suggesting that R2R3 MYB subfamily genes play important roles in cassava abscission zone development.
During early and middle leaf abscission, high levels of MeMYB15 and MeMYB73 expression were detected (Fig. 4). MYB15 decreases oxidative damage and increases proline in transgenic plants under salinity and dehydration conditions 9,18 . Under the stresses of salinity, dehydration and heat, many stress-responsive genes, such as LEA5, ERD10D, PLC3, LTP1, HSF2, ADC, P5CS, SOD and CAT, have been reported to act down stream of MYB15 9,18 . Proline and ROS regulate the development of cassava abscission 25 , and ADC and P5CS have been reported to up-regulate during the separation of cassava abscission zones 25 . Ethylene signaling plays essential roles in mediating plant responses to biotic and abiotic stresses. MYB73 also regulates antioxidant enzyme activity 24 . In addition, wheat Myb73 overexpression leads to the activation of the ABA-independent pathway genes RAB18 and CBF3, as well as the ABA-dependent genes ABF3 and RD29B, thereby improving salinity stress tolerance in Arabidopsis 24 . In addition, AtMYB15 expression is increased in response to the constitutive expression of ACO1 in Arabidopsis 18 , thus suggesting that ethylene signaling can regulate the function of MYB15. In addition, promoter motif analysis has also confirmed that MYB15 responds to the ethylene pathway 32 . All these results suggest that MeMYB15 regulates the separation of cassava abscission zones. High levels of MeMYB15 and MeMYB73 expression were detected at the early and middle leaf abscissions, thus suggesting that MYB genes participate in the regulation of abscission zone separation through ROS pathways and proline production.
Later in abscission, 7 R2R3 MYB genes, including 1 MeMYB78, 1 MeMYB3, 2 MeMYB4 and 3 MeMYB62 genes, were highly expressed between T2 and T6 in response to both ethylene and water-deficit stress treatment (Fig. 4). The genes that exhibited this expression pattern primarily contributed to cassava leaf abscission. MYB78 regulates reactive oxygen species (ROS)-dependent cell death 8   and GST are most prominently increased, with an average of 3-to greater than 10-fold increased expression in transgenic plants compared with controls 8,23 . In our previous work, we have found that the ROS levels decrease after co-overexpression of the ROS-scavenging proteins SOD and CAT1 in cassava. ROS also regulates cassava abscission zone separation 25 . In addition, MeMYB78 regulates ROS-dependent cell death, thus suggesting that it participates in the process of PCD in abscission zone. MeMYB3 regulates ROS via antioxidant biosynthesis pathways 33 . MYB3 also contains a conserved MYB domain and an ethylene responsive element binding factor-associated amphiphilic repression (EAR) repression domain 17,33 , thus suggesting that this gene is involved in the response to ethylene. MYB4 has been reported to respond to salicylic acid, ethylene, abscisic acid and methyl jasmonate hormones and also to enhance cellular antioxidant capacity through radical scavenging mechanisms and increased activities of phenylpropanoid and isoprenoid metabolic processes involving various abscisic acid (ABA), jasmonic acid (JA), salicylic acid (SA), ethylene and reactive oxygen species (ROS) responsive genes 15,16 . MYB62 regulates responses to environmental stresses via changes in GA metabolism and signaling 21 . GA metabolism and signaling have been reported to induce programmed cell death (PCD) in wheat aleurone cells 34 . The expression of ASYMMETRIC LEAVES1 (AS1), a R2R3 MYB transcription factor from Arabidopsis, is also regulated by the gibberellin (GA) pathway. AS1 regulates abscission zone placement in Arabidopsis flowers via restricting expression of the KNOX genes and BREVIPEDICELLUS (BP) in the sepals. Moreover, abscission of the medial sepals is delayed in as1 flowers 35,36 . Cassava MYB62 and Arabidopsis AS1 have been found to have high sequence similarity through blast analysis. KNOX genes regulate organ abscission via the IDA peptide signaling pathway through the HAESA (HAE) and HAESA-LIKE2 receptor-like kinases 35,36 . We have also confirmed that GA regulate cassava abscission zone separation 25 . MeMYB62 potentially participates in regulating cassava abscission zone separation.
The expression profiles of the genes acting downstream of the selected MYBs indicate involvement in cassava abscission zone separation. The above analyses identified many genes acting downstream of the selected MYBs; therefore, we analyzed the involvement of the putative downstream genes in the process of cassava AZ. CATs acted downstream of both MYB78 and MYB15, which regulate reactive oxygen species (ROS)-dependent cell death in AZ separation 8,25 . Two CAT genes were confirmed to express in abscissions in both ethylene and water-deficit stress treatments (Fig. 7). These two genes were expressed at high levels at late stages in both ethylene-and water-deficit stress-induced abscission. The highest expression ratio of CAT1 and CAT2 appeared at the T6 time point at 2.3-fold and 2.7-fold increased expression, respectively, compared with the T1 time point (Fig. 7). KNOXs act downstream of R2R3 MYB 35,36 , and negatively regulate the development of AZs. Two KNOXs were confirmed to be downregulated in abscissions in both ethylene and water-deficit stress treatments (Fig. 7). The two genes were expressed at low levels at late stages in both ethylene-and water-deficit stress-induced abscissions. The lowest expression ratio of KNOX1 and KNOX2 appeared at the T6 time point at 0.42-fold and 0.5-fold reduced expression, respectively, compared with the T1 time point (Fig. 7). ADC, P5CS and ACO1 act downstream of MYB73 and MYB15 9,18 and regulate cassava AZ abscission, as noted in our previous work 25 . MYB3 and MYB4 respond to ethylene, and ethylene has also been confirmed to increase in the process of abscission in cassava AZs 25 . All these genes that regulate AZ separation acted downstream of the selected MYBs and exhibited distinct expression in both ethylene-and water-deficit stress-induced abscissions in cassava. These results indicated that the selected MYBs regulate cassava AZ separation via these functional genes. In conclusion, this study sought to identify the R2R3 MYB genes that regulate cassava abscission zone development. More R2R3 MYB subfamilies were identified in the cassava genome than in the Arabidopsis genome. Comparative analysis of transcriptome expression profiles in response to ethylene and water-deficit stress indicated that most R2R3 MYB genes were expressed in response to both treatments. These results suggested that the same R2R3 MYB genes may function in leaf abscission in response to different treatments. Comparisons of R2R3 MYB genes expressed at the same time points during abscission in response to both treatments indicated that 9 R2R3 MYB subfamily genes were highly expressed during leaf abscission. Further analysis of promoter cis-elements confirmed that the R2R3 MYB subfamily responds to ethylene and regulates cassava abscission zone development. These data should provide a foundation for future research on the functional characterization of R2R3 MYB genes and signal transduction pathways.

Materials and Methods
Identification of cassava R2R3 MYB family genes. Cassava genome sequences suggested to contain a MYB domain were isolated and identified as candidate R2R3 MYB genes by using the Arabidopsis genome as a refs 3,37. R2R3 MYB genes were isolated from the cassava genome (JGI database, version 6.1) using annotations and BLAST analysis with an E-value cutoff set as 0.00001 3 . The Arabidopsis R2R3 MYB genes were identified using the Arabidopsis Information Resource (TAIR) 3,37 . Each putative cassava R2R3 MYB gene was searched against the TAIR database by using BLAST to ensure that no additional related genes were selected 3,37 . Phylogenetic analysis. The phylogenetic tree was constructed using neighbor-joining methods.
Phylogenetic analysis was performed using MEGA software, version 5. The resulting tree was tested for reliability using bootstrapping with 1,000 replicates and amino acid p-distance parameters 3,38 . Gene motif detection in cassava. In order to investigate the conserved motifs other than MYB repeats in the MYB protein sequences, we used the program of Multiple Em for Motif Elicitation (MEME; version 4.9.0) tool to analyze all of the full-length protein sequences of 166 MeMYBs 3 . The following parameter settings were used: distribution of motifs, zero or one per sequence; maximum number of motifs to find, 10; minimum width of motif, 6; maximum width of motif, 300 (to identify long R2R3 domains). Other options used the default values. Only motifs with an e-value of < 1e-20 were retained for further analysis 3 .
Plant materials and treatments. SC5 plants were grown as previously described 25,30 . In detail, cassava plants were planted in plastic pots at 28 °C under a 16-h light photoperiod (130 μ mol·m −2 ·s −1 ) for 3 months in a greenhouse 25 . Three plants were planted in one pot, and three pots served as a biological replicate. Three-month-old cassava plants with uniform growth statuses were chosen for ethylene and water-deficit stress treatments. Water-deficit stress and ethylene treatments were evaluated using the chlorophyll fluorescence parameter Fv/Fm 25 . For ethylene treatment, leaves were sprayed with 100 μ M ethylene. Water-deficit stress treatment plants were planted in a pot without water. Fv/Fm values were used to select six time points for AZ sample collection. Samples (approximately 1-2 mm) were cut from each pulvinus-petiole, including the AZs. AZs were collected from the middle of the cassava plants in each pots, and three pots served as a biological replicate. The procedure was repeated 3 times for each time point. AZs that were approximately 1-2 mm in size were frozen in liquid nitrogen for RNA extraction. RNA was purified using the plant RNA reagent from Invitrogen (Carlsbad, CA) to produce a pure and high-quality RNA preparation as indicated by spectroscopic and gel electrophoresis analysis 25 .
The cassava genome microarray: design, hybridization, and data analysis. To compare the gene expression profiles during ethylene and water-deficit stress-induced cassava leaf abscission, leaf chlorophyll fluorescence (Fv/Fm) values were measured and represented as six time points (T1-T6) as described previously 25 . To verify the reliability and accuracy of differential gene expression profiling, AZs samples were collected with the same Fv/Fm values for both ethylene-and water-deficit stress-induced leaf abscission. Then used time course cassava genome microarray analysis to examine the AZ gene expression during ethylene-and water-deficit stress-induced leaf abscission using the T1 time point as a control. A time series of cassava microarray analyses based on the principle of the "loop design" was performed as previously described 39 . For either ethylene or drought treatment, 18 distinct AZ samples (three biological replicates at each of 6 time points) to be compared, the experimental design included 36 two-color microarray slides for both treatments, allowing three technical replicates of each sample to be observed. The cassava microarray was constructed as previously described 25,39 . In detail, Two public databases were used for cassava microarray construction: the great majority of the ESTs originated from JGI database (http://www.phytozome.net/cassava.php) and the minority based on sequences from NCBI with E < 1e-5. Custom-designed 60-mer nimblegen DNA microarrays were synthesized by maskless in situ photolithographic synthesis 25 . The fluorescent dye (Cy3-dCTP)-labeled cassava cDNA was produced as previously described using CapitalBio cRNA Amplification and Labeling Kit (CapitalBio). After completion of double-stranded cDNA (dsDNA) synthesis, the dsDNA products were purified using a PCR NucleoSpin Extract II Kit (MN). The resulting cRNA was labeled according to Nimblegen recommendations. The procedures of Array hybridization, washing, scanning and data analysis were performed at CapitalBio Corporation (Beijing, China) according to the NimbleGen's Expression user's guide. The expression data of probes were normalized using quantile normalization and expression data of genes were generated using the Robust Multichip Average (RMA) algorithm 25 . Time course analysis. To identify differences in R2R3 MYB gene expression, statistical analysis was used to screen the microarray data for genes that were differentially expressed in the treatment samples (T2-T6) compared with the T1 control. For comparative analysis, differences in gene expression between samples (T2-T6) and reference (T1) were identified using significant analysis of microarray software (SAM, version 3.02) 25,39 . Changes in gene expression exceeding a threshold of < 0.5 or > 2.0-fold change with a Wilcoxon Rank-Sum test significance level of 0.01 (P < 0.01) and a false discovery rate (FDR) threshold of < 1% in the SAM output were considered to be differentially expressed. Time-dependent differentially expressed genes were classified with self-organizing tree algorithm clustering (SOTA) using MeV 4.0 software 25,40 . GO analysis. GO annotation of gene clusters was performed using BiNGO according to Maere et al. 25,41 .
Significant GO categories were identified using a hypergeometric test with a significance threshold of 0.01 after a Benjamini-Hochberg FDR correction 42 . GO categories were classified by hierarchical clustering using MeV 4.0 software.
Analysis of cassava R2R3 MYB subfamily putative promoter region cis-elements. The putative promoter regions were analyzed for cis-elements according to Wu et al. 43 . In detail, the 2-kbp putative promoter regions upstream of each R2R3 MYB subfamily coding DNA sequence were examined for cis-elements by using the PLACE website (http://www.dna.affrc.go.jp/PLACE/).

Real-time RT-PCR validation of R2R3 MYB subfamily and the downstream gene expression in cassava.
RNA from three independent biological samples were reverse transcribed and used for real-time qRT-PCR with SYBR Green I (Carlsbad, CA) detection on a STEP-ONE system. Gene-specific primers were designed using IDT Primerquest tools (http://www.idtdna.com/Primerquest/Home/Index) 25 . The SYBR Green PCR kit (Applied Biosystems) was used to perform QPCR 44 . The real-time PCR primer sequences are listed in Additional table 6 and 7. To avoid non-specific amplification from other family genes, the primers were designed to span intron-exon boundaries or target the untranslated regions and the primer pairs were confirmed using the cassava genome database to ensure their specificity. The expression ratios of R2R3 MYB genes and the downstream genes at each time point are presented using T1 as a control first in SC5, and then to cut the T1 control expression ratios with those of SC5 at each time point. The QPCR procedures were as follows: 10 min of denaturation at 95 °C, followed by 40 cycles of amplification with 15 sec of denaturation at 94 °C, 30 sec of annealing according to the melting temperatures provided in Supplemental Table S2, 35 sec of extension at 72 °C, and the fluorescence data collection at 72 °C. After a final extension at 72 °C for 10 min, the specificity of the amplified product was evaluated by melting curve 44 .