The R2R3-MYB gene family in banana (Musa acuminata): Genome-wide identification, classification and expression patterns

The R2R3-MYB genes comprise one of the largest transcription factor gene families in plants, playing regulatory roles in plant-specific developmental processes, defense responses and metabolite accumulation. To date MYB family genes have not yet been comprehensively identified in the major staple fruit crop banana. In this study, we present a comprehensive, genome-wide analysis of the MYB genes from Musa acuminata DH-Pahang (A genome). A total of 285 R2R3-MYB genes as well as genes encoding three other classes of MYB proteins containing multiple MYB repeats were identified and characterised with respect to structure and chromosomal organisation. Organ- and development-specific expression patterns were determined from RNA-Seq data. For 280 M. acuminata MYB genes for which expression was found in at least one of the analysed samples, a variety of expression patterns were detected. The M. acuminata R2R3-MYB genes were functionally categorised, leading to the identification of seven clades containing only M. acuminata R2R3-MYBs. The encoded proteins may have specialised functions that were acquired or expanded in Musa during genome evolution. This functional classification and expression analysis of the MYB gene family in banana establishes a solid foundation for future comprehensive functional analysis of MaMYBs and can be utilized in banana improvement programmes.


Introduction
Banana (Musa spp.), including dessert and cooking types, is a staple fruit crop for a major world population, especially in developing countries. The crop is grown in more than 100 countries throughout the tropics and sub-tropics, mainly in the African, Asia-Pacific, and Latin American and Caribbean regions [1]. Bananas provide an excellent source of energy and are rich in certain minerals and in vitamins A, C and B6. Furthermore, this perennial, monocotyledonous plant provides an important source of fibre, sugar, starch and cellulose (used for paper, textiles). Bananas have also been considered as a useful tool to deliver edible vaccines [2]. Certain agronomic traits, such as stress and pest resistance as well as fruit quality, are thus

Search for MYB protein coding genes in the M. acuminata genome
A consensus R2R3-MYB DNA binding domain sequence [12] (S1 Table) was used as protein query in tBLASTn [13] searches on the M. acuminata DH-Pahang genome sequence (version 2) (https://banana-genome-hub.southgreen.fr/sites/banana-genome-hub.southgreen.fr/files/ data/fasta/version2/musa_acuminata_v2_pseudochromosome.fna) in an initial search for MYB protein coding genes. To confirm the obtained amino acid sequences, the putative MYB sequences were manually analysed for the presence of an intact MYB domain. All M. acuminata MYB candidates from the initial BLAST were inspected to ensure that the putative gene models encode two or more (multiple) MYB repeats. The identified gene models were analysed to map them individually to unique loci in the genome and redundant sequences were discarded from the data set to obtain unique MaMYB genes. The identified MaMYB genes were matched to the automatically annotated genes from the Banana Genome Hub database [14]. The open reading frames of the identified MaMYB were manually inspected and, if available, verified by mapping of RNA-Seq data to the genomic sequence. Resulting MaMYB annotation is provided in the supporting information: Multi-FASTA files with MaMYB CDS sequences (S1 File) and MaMYB peptide sequences (S2 File), general feature format (GFF) file (S3 File) describing the MaMYB genes for use in genome viewers/browers.

Genomic distribution of MaMYB genes
The MaMYB genes were located on the corresponding chromosomes by the MapGene2chromosome web v2 (MG2C) software tool (http://mg2c.iask.in/mg2c_v2.0/) according to their position information of a physical map, available from the DH-Pahang genome annotation (V2).

Phylogenetic analyses
Protein sequences of 133 A. thaliana MYBs were obtained from TAIR (http://www. arabidopsis.org/). We also considered other multiple MYB-repeat proteins from A. thaliana in the phylogenetic analysis to determine orthologs in the M. acuminata genome: five AtMYB3R, AtMYB4R1 and AtCDC5. Additionally, 43 well-known, functionally characterised landmark plant R2R3-MYB protein sequences were collected from GenBank at the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/). Phylogenetic trees were constructed from ClustalO [15] aligned MYB domain sequences (293 MaMYBs, 132 AtMYBs and 43 plant landmark MYBs) using MEGA7 [16] with default settings. A mojority rule Maximum Likelihood (ML) consensus tree inferred from 1000 bootstrap replications was calculated. M. acuminata MYB proteins were classified according to their relationships with corresponding A. thaliana and landmark MYB proteins.

Results and discussion
The first annotated reference genome sequence of M. acuminata (A genome) became available in 2012 [3]. It was obtained from a double haploid (DH) plant of the Pahang cultivar, derived through haploid pollen and spontaneous chromosome doubling from the wild subspecies Musa acuminata ssp. malaccensis [3]. This wild subspecies was involved in the domestication of the vast majority of cultivated bananas and its genetic signature is commonly found in dessert and cooking bananas (ProMusa, http://www.promusa.org). An improved version (DH-Pahang version 2) of the genome assembly and annotation was presented in 2016 comprising 450.7 Mb (86% of the estimated size) from which 89.5% are assigned to one of the 11 chromosomes, predicted to contain 35,276 protein encoding genes [18]. This M. acuminata DH-Pahang version 2 genome sequence provides the platform of this study.

Identification and genomic distribution of M. acuminata MYB genes
A consensus R2R3-MYB DNA binding domain sequence (deduced from Arabidopsis, grape and sugarbeet R2R3-MYBs, S1 Table) was used as protein query in tBLASTn searches on the DH-Pahang version 2 genome sequence to comprehensively identify MYB protein coding genes in M. acuminata. The resulting putative MYB sequences were proven to map to unique loci in the genome and were confirmed to contain an intact MYB domain. This ensured that the gene models contained two or more (multiple) MYB repeats. We identified a set of 285 R2R3-MYB proteins and nine multiple repeat MYB proteins distantly related to the typical R2R3-MYB proteins: six R1R2R3-MYB (MYB3R) proteins, two MYB4R proteins and one CDC5-like protein from the M. acuminata genome sequence ( Table 1).
The number of R2R3-MYB genes is one of the highest among the species that have been studied to date, ranging from 45 in Ginkgo biloba [10] over 157 in Zea mays [25] and 249 in Brassica napus [26] to 360 in Gossypium hirsutum [11]. This is probably due to three wholegenome duplications (γ 100 Myr ago and α, β 65 Myr ago) that occurred during Musa genome evolution [3,27]. The number of atypical multiple repeat MYB genes identified in M. acuminata is in the same range as those reported for most other plant species, up to six MYB3R and up to two MYB4R and CDC5-like genes.
The 285 MaR2R3-MYB genes identified constitute approximately 0.81% of the 35,276 predicted protein-coding M. acuminata genes and 9.0% of the 3,155 putative M. acuminata transcription factor genes [18]. These were subjected to further analyses. The identified MaR2R3-MYB genes were named following the nomenclature of the locus tags provided in the DH-Pahang version 2 genome annotation (Table 1). A keyword search in the NCBI database (http://www.ncbi.nlm.nih.gov/) revealed no evidence for M. acuminata MYB genes not present in Table 1.
Six publications dealing with M. acuminata MYB genes were identified: one study describes the elevated expression of nine MaMYB genes in transgenic banana plants overexpressing the NAC domain transcription factor MusaVND1 (vascular related NAC domain) indicating a role of these MaMYBs in the regulation of secondary wall deposition [22]. MYBS3 (Ma08_25960) was found to be differentially expressed between cold-sensitive and cold-tolerant bananas [24]. Another publication described the M. acuminata R2R3-MYB gene Ma05_03690 being upregulated in the early response to the endoparasitic root-knot nematode Meloidogyne incognita in roots [28]. MusaMYB31 (Ma01_02850) was identified as a negative regulator of lignin biosynthesis and the general phenylpropanoid biosynthesis pathway [21]. MaMYB3 (Ma06_11140) was found to repress starch degradation in fruit ripening [23] and MaMYB4 (Ma01_19610) was recently described to control fatty acid desaturation [29]. On the basis of the DH-Pahang version 2 annotation, 291 of the 294 MaMYB genes could be assigned to the eleven chromosomes. The chromosomal distribution of MaMYB genes on the pseudochromosomes is shown in Fig 1 and revealed that M. acuminata MYB genes are distributed across all chromosomes.
Gene structure analysis revealed that most MaR2R3-MYB genes (155 of 285; 54.4%) follow the previously reported rule of having two introns and three exons, and display the highly conserved splicing arrangement that has also been reported for other plant species [25,30,31]. A total of 23 (8.1%) MaR2R3-MYB genes have two exons and seven (2.5%) were single exon genes. 67 (23.5%) MaR2R3-MYB genes have four exons, 21 (7.4%) five exons, eight (2.8%) six exons and two of each with seven and twelve exons, respectively. The complex exon-intron structure of Ma06_12110 and Ma10_18840 is conserved in their A. thaliana orthologs AtMYB88 and AtMYB124/FOUR LIPS (FLP) containing ten and eleven exons, respectively. This supports their close evolutionary relationship, but also suggests the conservation of this intron pattern in evolution since the monocot-dicot split 140-150 Myr ago [32].
Chromosomes are drawn to scale. The positions of centromeres (grey ovals) are roughly estimated from repeat distribution data. The chromosomal positions of the MaMYB genes (given in DH Pahang version 2 annotation ID) are indicated. R2R3-MYB genes are given in black letters, MYB genes with more than two MYB repeats are given in blue letters. The number of MYB genes on each chromosome is given below the respective chromosome. The genes are ordered by DH Pahang version 2 [18] pseudochromosomes, from north to south. The annotation-version specific gene code describing the chromosomal assignment and position on pseudochromosomes is given in the first column. "Ma00" indicates genes located in sequences without chromosomal assignment.

Phylogenetic analysis of the M. acuminata MYB family
With the aim to explore the putative function of the predicted M. acuminata MYBs, we assigned them to plant MYB proteins with known function. For this, we chose primarily data from A. thaliana, which is the source of most functional MYB characterisations. From comparable studies, MYB function appears conserved across MYB clades, suggesting that closely related MYBs recognise similar/same target genes and possess cooperative, overlapping or redundant functions.
To unravel the relationships, we constructed a phylogenetic tree with 468 MYB domain amino acid sequences of MYB proteins. We used 293 MaMYBs (omitting the CDC5-like MaMYB), the complete A. thaliana MYB family (132 members, including 126 R2R3-MYB, five MYB3R and one MYB4R) and 43 functionally well characterised landmark R2R3-MYBs from other plant species. The phylogenetic tree topology allowed us to classify the analysed MYBs into one MYB3R clade, one MYB4R clade and 42 R2R3-MYB protein clades (Fig 2).

PLOS ONE
Most R2R3-MYB clades (31 of 42) include variable numbers of MYB proteins from Arabidopsis and banana, indicating that the appearance of most MYB genes in these two species predates the monocot-dicot split as observed in other studies [33]. Several of these clades also contain landmark MYBs from other plant species. The two clades 4 and 28 only contain Arabidopsis R2R3-MYB members, while seven clades (9,13,19,21,26,29 and 32, displayed with white background in Fig 2) only contain banana R2R3-MYB members. Additionally, two clades, 24 with

PLOS ONE
monocot anthocyanin biosynthesis regulators and 27 with the strawberry landmark flavonoid biosynthesis repressor FaMYB1 [34], do not contain any A. thaliana MYB (Fig 2).
ML consensus tree inferred from 1000 bootstraps with 468 MYB domain amino acid sequences of MYB proteins from Musa accuminata (Ma), Arabidopsis thaliana (At) and landmark MYBs from other plant species built with MEGA7. Clades are labeled with different colors and functional annotations are given. The numbers at the branches give bootstrap support values from 1000 replications. SCW, secondary cell wall.
The lineage specificity of some MYB clades could indicate that these clades may have been lost or gained in a single order or species during plant evolution, as indicated by other studies [35][36][37]. For example, clade 28 lacked M. acuminata orthologs, but includes the A. thaliana R2R3-MYBs AtMYB0/GLABRA1 and AtMYB66/WEREWOLF, which have been identified as being involved in the formation of trichomes and root hairs from epidermal cells [38,39]. Similar observations have been made in maize (monocot) and sugarbeet (eudicot, caryophyllales) [30,31], both not containing clade 28 orthologs, while grape (eudicot, rosid) and poplar (eudicot, rosid) do [12,40]. It has been hypothesized that GLABRA1-like MYB genes have been acquired in rosids after the rosid-asterids split [41]. The absence of MaMYBs in clade 28 is consistent with this hypothesis, since monocots branched off before the separation of asterids and rosids in eudicots. Clade 4 also lacks banana R2R3-MYBs. This clade contains the glucosinolate biosynthesis regulators AtMYB28, AtMYB34 and AtMYB51 [42,43]. The absence of MaMYBs in this clade is concordant with the fact that glucosinolates are only present in the Brassicaceae family. This clade is thought to have originated from a duplication event before the divergence of the genus Arabidopsis from Brassica [44].
The seven clades containing only MaMYBs were manually inspected by applying BLAST searches at the NCBI protein database in order to identify high similarity to functionally characterized landmark plant R2R3-MYBs. In no case could a landmark MYB be identified. Consequently, these clades could be described as a lineage-specific expansion in M. acuminata, reflecting a species-, genus-or order-specific evolutionary change. These MaMYB proteins may have specialised functions that were acquired or expanded in M. acuminata during genome evolution. Further research will be needed to decipher the biological roles of these MaMYB genes.

Expression profiles for M. acuminata MYB genes in different organs and developmental stages
Since it is not unusual for large transcription factor families in higher organisms to have redundant functions, a particular transcription factor needs to be studied and characterised in the context of the whole family. In this regard, the gene expression pattern can provide a hint to gene function. In this study we used publicly-available RNA-Seq reads from the Sequence Read Archive (S2 Table) to analyse the expression of the 293 MaMYB genes in different organs and developmental stages: embryonic cell suspension, seedling, root and young, adult and old leaf, pulp stages S1-S4 and peel stages S1-S4. Filtered RNA-Seq reads were aligned to the reference genome sequence and the number of mapped reads per annotated transcript were quantified and compared across the analysed samples to calculate normalised RNA-Seq read values which are given in Table 2.
Our expression analyses revealed that M. acuminata MYBs have diverse expression patterns in different organs. Many of the MaMYBs exhibited low transcript abundance levels with expression in only one or a few organs. This is consistent with other transcription factor genes, typically found to be expressed in this manner due to functional specificity and diversity. The highest number of expressed MaMYB genes (221; 75.4%) is observed in roots, followed by pulp (209; 71.3%), leaf (203; 69.3%), peel (196; 66.9%) and embryonic cell suspension (129; 44%). The fewest MaMYB genes are expressed in seedlings (98; 33.4%) in the considered dataset. 40 MaMYB genes (13.7%) were expressed in all samples analysed (albeit with varying expression levels), which suggested that these MaMYBs play regulatory roles at multiple developmental stages in multiple tissues. 13 MaMYB genes (4.4%) lacked expression information in any of the analysed samples, possibly indicating that these genes are expressed in other organs (e.g. pseudostem, flower, bract), specific cells, at specific developmental stages, under special conditions or are pseudogenes. 280 MaMYBs (95.6%) are expressed in at least one analysed sample, although the transcript abundance of some genes was very low. Some MaMYB genes were expressed in all analysed RNA-Seq samples at similar levels (e.g. R2R3-MYBs Ma06_g33440 and Ma11_g06880) while others show variance in transcript abundance with low (no) levels in one or several organs and high levels in others (or vice versa). For example, Ma02_g16570, Ma03_g09310, Ma07_g11110, Ma10_g01730, Ma10_g05260 and Ma10_g09100 show organ-specific expression, as their transcripts were exclusively detected in leaves, which hints to leaf-specific functionality. Ma03_g07840 and Ma07_g19470 were found to be predominantly expressed in pulp, showing expression also in other analysed organs, but not in the seedling. Overall, these results suggests that the corresponding MaMYB regulators are limited to distinct organs, tissues, cells or conditions. Some paralogous MaMYB genes clustered in the genome (Tables 1 and 2) showed different expression profiles, while other clustered paralogous MaMYB genes did not. For example, the clustered R2R3-MYB genes Ma07_g19880 and Ma07_g19890 (both in the proanthocyanidinrelated cluster 25): while Ma07_g19880 is expressed in pulps, peels, roots, and leaves, Ma07_g19890 is nearly not expressed in the analysed organs. The expression pattern of Ma03_g07840 and Ma03_g07850 (also both in the proanthocyanidin-related cluster 25) is, in contrast, very similar, with low expression in root, leaf, pulp and peel, but no expression in embyonic cells and roots. These results could point to functional redundancy of the genes Ma03_g07840 and Ma03_g07850, while Ma07_g19880 and Ma07_g19890 could be (partly) involved in distinct, tissue-specific aspects of specialised mebabolite biosynthesis. The functional categorisation and its expression domains make Ma07_g19880 a good candidate to encode a proanthocyanidin biosynthesis regulator in developing banana fruits (pulp and peel), and thus an excellent target for genetic manipulation of proanthocyanidin content, known to influence senescence of banana fruit [46], to maintain the freshness of harvested banana fruit.

Conclusions
The present genome-wide identification, chromosomal organisation, functional classification and expression analyses of M. acuminata MYB genes provide a first step towards cloning and functional dissection to decode the role of MYB genes in this economically interesting species. Further, knowledge of selected MYB genes can be biotechnologically utilized for improvement of fruit quality, yield, disease resistance, tolerance to biotic and abiotic stresses, and the biosynthesis of pharmaceutical compounds.
Supporting information S1 Table. MYB domain consensus sequence of R2R3-MYBs, used for tBLASTn searches.