“A draft Musa balbisiana genome sequence for molecular genetics in polyploid, inter- and intra-specific Musa hybrids”

Background Modern banana cultivars are primarily interspecific triploid hybrids of two species, Musa acuminata and Musa balbisiana, which respectively contribute the A- and B-genomes. The M. balbisiana genome has been associated with improved vigour and tolerance to biotic and abiotic stresses and is thus a target for Musa breeding programs. However, while a reference M. acuminata genome has recently been released (Nature 488:213–217, 2012), little sequence data is available for the corresponding B-genome. To address these problems we carried out Next Generation gDNA sequencing of the wild diploid M. balbisiana variety ‘Pisang Klutuk Wulung’ (PKW). Our strategy was to align PKW gDNA reads against the published A-genome and to extract the mapped consensus sequences for subsequent rounds of evaluation and gene annotation. Results The resulting B-genome is 79% the size of the A-genome, and contains 36,638 predicted functional gene sequences which is nearly identical to the 36,542 of the A-genome. There is substantial sequence divergence from the A-genome at a frequency of 1 homozygous SNP per 23.1 bp, and a high degree of heterozygosity corresponding to one heterozygous SNP per 55.9 bp. Using expressed small RNA data, a similar number of microRNA sequences were predicted in both A- and B-genomes, but additional novel miRNAs were detected, including some that are unique to each genome. The usefulness of this B-genome sequence was evaluated by mapping RNA-seq data from a set of triploid AAA and AAB hybrids simultaneously to both genomes. Results for the plantains demonstrated the expected 2:1 distribution of reads across the A- and B-genomes, but for the AAA genomes, results show they contain regions of significant homology to the B-genome supporting proposals that there has been a history of interspecific recombination between homeologous A and B chromosomes in Musa hybrids. Conclusions We have generated and annotated a draft reference Musa B-genome and demonstrate that this can be used for molecular genetic mapping of gene transcripts and small RNA expression data from several allopolyploid banana cultivars. This draft therefore represents a valuable resource to support the study of metabolism in inter- and intraspecific triploid Musa hybrids and to help direct breeding programs.


Background
Bananas and plantains (Musaceae, Zingiberlaes, Musa ssp., 'bananas') are giant monocotoyledenous herbs, which originated in Southeast Asia and the western Pacific. They were one of the first crops to be domesticated, and are now widely distributed throughout the subtropics where they constitute a major staple food for millions of people [1]. Indeed, bananas are the World's most important fruit and rank fourth in the list of global food crops [2], but with the important caveat that 80% of the production is for local consumption [3].
There are four genomes present within Musa spp., corresponding to the genetic constitutions of the four wild Eumusa species i.e. Musa acuminata (A-genome, 2n = 2x = 22), Musa balbisiana (B-genome, 2n = 2x = 22), Musa schizocarpa (S genome, 2n = 2x = 22) and the Australimusa species (T genome, 2n = 2x = 20). Within M. acuminata, a number of subspecies are also recognized which have been further classified into four subgroups on the basis of DNA markers i.e. banksii, zebrina, malaccensis and burmannica/burmannicoides [4][5][6]. Despite evidence of interspecific genetic variation, there has been no subspecies classification within M. balbisiana [7]. Members of the Australimusa sections have a basic chromosome number of 2n = 20, and comprise seven species one of which contains edible parthenocarpic types known as 'Fe'i' bananas. Members of this Fe'i group are notable for their upright fruit and in some cases exceptionally high fruit provitamin A carotenoid contents [8,9].
Cultivated bananas are parthenocarpic, generally seedless and vegetatively propagated hybrids that have arisen primarily as a result of hybridizations between wild diploid M. acuminata and M. balbisiana species [10][11][12]. Intraspecific hybridizations within M. acuminata, and interspecific hybridizations between M. acuminata and M. balbisiana have naturally resulted in various combinations of these A-and B -genomes which are classified into six groups (AA, AAA, AB, AAB, ABB and ABBB) [11]. The majority of edible cultivars though are allopolyploid triploids with a genome constitution of AAA (dessert banana), AAB (plantains) and ABB (cooking bananas).
It has also become increasingly clear that there is no simple division of parental A-and B-chromosomes during hybridization, and several research groups have provided evidence for pairing and recombination between homeologous A-and B-chromosomes [13,14]. Similar events have been observed in allopolyploids of other species [15]. It seems likely therefore that current interspecific triploids have arisen from one or more steps of (re)combination and exchange of chromosomal segments between the A-and B-genomes [13,16,17]. As a result, most if not all Musa cultivars probably have genomes consisting of different proportions of the A-and B-genome. A similar process of hybridization between subspecies of M. acuminata probably also underlines the evolution of the edible AA and AAA types. A consequence of these recombination events is that the hybrid genomes contain an unbalanced number of A-and B-genome alleles [18]. This clearly complicates genetic studies of trait inheritance, as well as the development and application of molecular marker technologies [17].
Currently, over 90% of commercial export dessert bananas are produced from a single cultivar, namely 'Cavendish' (AAA). Unsurprisingly, this dependency on a single cultivar and the consequent lack of genetic variation in production systems, has resulted in a crop which is potentially highly susceptible to disease pandemics. Pandemics are also not without precedent in banana, and in the 1950's, 'Panama Disease' ('Fusarium wilt') caused by the soil-born fungus Fusarium oxysporum was responsible for entirely wiping out commercial banana production, which at that time was dependent on the cultivar 'Gros Michel' (AAA). By necessity, 'Gros Michel' was replaced by the Fusarium-resistant cultivar 'Cavendish', but unfortunately, in the intervening years new strains of F. oxysporum have evolved, (F. oxysporum f. sp. cubense race 4, 'Tropical race 4'), that have overcome this natural Cavendish resistance [19]. Currently, the other major disease threat is that posed by Mycosphaerella fijiensis ('Black Sigatoka'). Chemical control of Fusarium wilt is ineffective, and even though 'Black Sigatoka' can be controlled by frequent (up to 50) fungicide applications per year, these are socio-economically and environmentally inappropriate, and require integrated strategies to avoid the development of fungicide resistance in the pathogen [20,21].
Consequently, there is currently much interest in developing sustainable solutions to these disease threats through the introgression of novel resistance loci and the development of new, disease resistant varieties. However the largely sterile nature of the majority of commercial varieties means that banana breeding is both lengthy and time-consuming. Therefore the application of molecular marker technologies for the identification of trait-marker associations, and high throughput genotyping technologies can greatly accelerate the whole breeding cycle through marker assisted selection (MAS) [22]. In this regard the recent publication of a high quality 523 Mb draft reference genome for the doubled haploid genome of the wild M. acuminata variety 'Pahang' (AA), represents a major advance in this field [14]. The variety 'Pahang' is member of the subspecies malaccensis that contributed one of the three M. acuminata genomes for 'Cavendish' [4]. This invaluable resource therefore provides for the first time a complete catalogue of all predicted genes, transcripts and markers in Musa, and can greatly facilitate and accelerate the search for novel genes, transcripts, allelic variants etc., for important biological processes. It also opens up the prospects for the rapid development of high-throughput molecular tools for Musa improvement. Despite this, there is also an urgent need for the development of similar resources for the Musa B-genome, to be able to identify and exploit M. balbisiana accessions with resistance to various abiotic and biotic stresses including osmotic and cold stresses [23][24][25][26] as well as for vigour [27,28].
The aim of this work therefore was to sequence and assemble a draft genome of the wild M. balbisiana diploid variety 'Pisang Klutuk Wulung' ('PKW', B-genome) for use in comparative transcriptomics and genomics studies of interspecific triploid and tetraploid hybrids. PKW is one of the possible ancestral parents of the B-genome present in cultivated triploids, and has for example has been shown to have very strong partial resistance to black leaf streak virus (http://www.proMusa.org/tiki-index. php?page=Pisang+Klutuk+Wulung). Other M. balbisiana cultivars have demonstrated resistance to Xanthomonas [29], and are considered to be more drought tolerant [24,26]. The utility of the PKW B-genome sequence generated here was validated by examining the distribution of CDS, EST and RNA-read mappings from a range of genetically diverse interspecific triploid cultivars across the combined A-and B-genomes, by characterising the predicted miRNA sequences encoded in both A-and B-genomes and by prediction of miRNA targets in both genomes.
Improved insights into genome structure, allelic diversity and regulatory elements within M. balbisiana spp. will help in the design and application of breeding strategies for novel banana cultivars with improved stress resistance and quality traits.
Banana root tissues were prepared from clonal tissuecultured plantlets of M. acuminata cultivar 'Berangan' (AAA), when 6-8 cm tall, with three fully-expanded leaves and healthy roots. Root tissues were pooled from 3-4 plantlets for RNA isolation and small RNA sequencing. Banana embryogenic suspension cell samples were prepared from embryogenic callus induced from immature male flowers of 'Berangan' as described by [30].
Newly initiated suspension cells were used for RNA isolation and small RNA sequencing.

gDNA extraction
Genomic DNA (gDNA) was extracted from young leaves essentially as according to Michiels et al. [31].

RNA extraction
Total nucleic acids were extracted from lyophilized fruit pulp samples, from 3 individual fruit per cultivar using the Tris-LiCl procedure of Tattersall et al. [32], modified as previously described [33]. Equivalent amounts of RNA from each sample were combined per cultivar. Total nucleic acids were isolated from banana root tissues and embryogenic suspension cells using a modified CTAB nucleic acid isolation method [34].

RNA and DNA quality
Concentrations of purified nucleotides were determined at 260 nm using a NanoDrop 2000 Spectophotometer (Thermo Scientific) and purity assessed at an absorbance ratio of 260/280 nm and 260/230 nm. RNA integrity was confirmed by agarose gel electrophoresis and on an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.). Only samples with high RNA integrity number (RIN ≥ 8) were used for RNA sequencing. A total of 2 μg of purified gDNA and of each combined RNA/DNA sample was precipitated in ethanol and used for sequencing.

RNA and DNA sequencing
Sequencing of both RNA and DNA samples was carried out at the Genomics Core sequencing facilities of the Katholieke Universiteit Leuven (http://www.uzleuven.be/ genomicscore/), using an Illumina HiSeq 2000 II instrument. Small RNA libraries were sequenced using an Illumina HiSeq 2000 II at BGI, Shenzhen.

Illumina paired-end cDNA library construction and sequencing
The cDNA libraries were constructed using the TruSeq™ RNA Sample Preparation Kit (Illumina, Inc.) according to the manufacturer's instructions. Poly-A containing mRNA was purified from 2 μg of total RNA using oligo (dT) magnetic beads and fragmented into 200-500 bp pieces using divalent cations at 94°C for 5 min. The cleaved RNA fragments were copied into first strand cDNA using SuperScript II reverse transcriptase (Life Technologies, Inc.) and random primers. After second strand cDNA synthesis, fragments were end repaired, a-tailed and indexed adapters were ligated. The products were purified and enriched with PCR to create the final cDNA library. The 6 tagged cDNA libraries were pooled in equal ratios and used for 2 × 100 bp paired-end sequencing on a single lane of the Illumina HiSeq2000 II.

Sequence data processing
Sequencing data was provided as fastq files and unless otherwise mentioned, all data processing steps were carried out using the CLC Genomics Workbench software package v 6.01. Raw reads were uploaded to GenBank and are accessible via accession number SAMN02333823.
Mapping PKW gDNA reads to the reference A-genome The raw data was first trimmed to remove low quality bases and the trimmed PKW gDNA reads then aligned to the publically available M. acuminata 'Pahang' doubled haploid A-genome assemblies available from http:// banana-genome.cirad.fr/download.php. This consists of 11 chromosomes together with one sequence containing concatenated unassembled contigs, each separated by 100 'N's'. Reads were aligned using the settings; mismatch cost (2), insertion cost (3) deletion cost (3), length fraction (0.9) and similarity fraction (0.9). Reads mapping equally well to two positions were assigned randomly. Following mapping, the consensus sequences were extracted and served as the PKW consensus reference genome for further treatments. For chromosomes 1 -11, during extraction of the consensus PKW genome sequence, regions of 0 read coverage (e.g. gaps etc.) are removed to produce a single continuous sequence. For the large unmapped chromosome, the consensus PKW sequence was extracted using 'N' ambiguity symbols to fill in gap regions, as otherwise unrelated genic sequences could be concatenated together allowing 'bridging' of reads across unrelated genic sequences.
Mapping RNA de novo assembled transcripts, CDS and unigenes to gDNA contigs and genome sequences RNA reads (de novo assembled transcripts, CDS, unigenes), were aligned to the PKW genome or gDNA contig data using the 'large gap' mapping function within the CLC Genomics Workbench, using the following settings; Maximum number of hits for a segment = 10, Maximum distance from seed = 50,000, Key Match _Mode = random, Mismatch cost = 2, Insertion cost = 3, Deletion cost = 3, Similarity = 0.8, Length fraction = 0.9. The large gap mapper function aligns reads to a reference sequence, while allowing for large gaps in the mapping. It is therefore able to map reads that span introns without requiring prior transcript annotations or for the detection of large deletions in genomic data. Additional details can be found http:// www.clcbio.com/white-paper.

B-genome annotation
Ab initio gene prediction was carried out using the FGENESH software, available online from http://linux1. softberry.com/all.htm and using the default parameters and the 'monocot model plant' parameters. The list of predicted PKW gene models was then blasted against the NCBI nr protein database and gene ontology terms assigned using the Blast2Go software (Conesa et al., [35]; Götz et al., [36]). Repeats were annotated by BLAST against the repetitive part of the Musa genome containing 1902 sequences which were retrieved from a published report [37]. Evaluation of the PKW B-genome gene model set took place by large-gap mapping of available CDS, and EST resources within CLC Genomics Workbench. These resources consisted of the 'Pahang' consensus CDS set, an in-house Musa unigene set of 22,205 sequences derived from the Syngenta M. acuminata 3' EST database [33], transcript sets generated from the de novo assembly of Illumina 100 bp paired end RNA reads from 6 Musa cultivars.

De novo assembly
All the trimmed, PKW gDNA reads were de novo assembled using the default settings in CLC Genomics Workbench with the settings as follows; Word size: 25, Bubble size: 50, Minimum contig length = 200, Mismatch cost = 2, Insertion cost = 3, Deletion cost = 3, Length fraction = 0.5, Similarity fraction = 0.8. In addition the 'scaffolding' option was used to take advantage of the paired read data, and reads were mapped back to the contigs generated to validate the sequences. The same parameters were also used for the de novo assembly of PKW gDNA reads that could not be mapped to the 'Pahang' reference genome, and for the de novo assembly of all RNA transcript data from the 5 triploid hybrids and the diploid Fe'i variety.

Variant analysis
Sequence variant analysis was carried out on both the RNA-and DNA-mappings using the 'Probabilistic Variant Detection' plugin in CLC Genomics Workbench, with settings specifying a minimum read coverage of 10 and a variant probability of >90%. The maximum expected number of variants was set as 2 or 3 according to the ploidy level of the samples and variants were calculated using either all the mapped reads, or only using the uniquely mapped reads.

Repeat annotation
This was performed on both the 'Pahang' A-genome and the PKW consensus B-genomes (our data) as well as de novo assembled contigs with Repeat Masker V4.0.3 software tool (http://www.repeatmasker.org), [38] using RMBLAST 2.2.27 as the engine and using the customized library of M. acuminata repeats (1903 sequences) from Hribova et al. 2010 [37]. SSR detection was carried out using Tandem Repeat Finder (TRF) software [39] and TRAP [40] using the default parameters. microRNA (miRNA) and target prediction miRNA prediction was performed using miRDeep2 tool [41] using scripts modified according to the criteria set for plant genomes [42]. A non-redundant query set of small RNA reads was compiled from root and embryogenic cell suspension and included all 235 miRNA sequences reported for the Musa 'Pahang' doubled haploid A-genome retrieved from the banana genome database [14], and publically available small RNA data from M. acuminata 'Calcutta 4' leaf, flower and fruit tissues (sequenced within the framework of a NSF project http:// smallrna.udel.edu). Small RNA reads were trimmed to remove adapter sequences and low quality reads. Reads were then mapped to the Rfam database [43] using Bowtie [44]. Matches with tRNA, rRNA, small nucleolar RNA and sequences below 19 or above 24 nucleotides in length were not considered for further analysis. The query miRNA data set was mapped separately to the A-genome [14] and the draft PKW B-genome developed here. Regions of 300 nt surrounding the matched position of each read were excised from the genome sequences, and then RNAfold software [45] used to predict sequences able to form stem-loop structures, using default options, whilst 'Randfold' [46] was used to calculate p-values for potential miRNA precursors predicted by the 'miRDeep2' algorithm. The candidate miRNA precursors selected had the following features; a predicted stem loop structure of 75 nt and a bulge-loop size of less than 6 nt; the mature miRNA was within the stem region of the precursor; less than four mismatches were allowed between the mature miRNA:miRNA* duplex; miRNA and miRNA* were on opposite arms of the precursor forming a duplex with 3' overhangs; the predicted minimum folding energy (MFE) was between -15 kcal/mol to −47.2 kcal/mol. Known miRNA sequences/homologues were annotated by BLASTn comparison [47] to mature and stem loop miRNA sequences from miRBase v19 [48]. Predicted miRNA were considered novel if they had no match (allowing for a maximum of 2 mismatches i.e. n/n, (n-1)/n, (n-2)/n nucleotide matches, n = length of mature miRNA) to any entry in miRBase (release 19) and PMRD (accessed, February 2013). Novel Musa miRNA sequences not present in either miRBAse or PMRD databases, were arbitrarily named starting at '1' and using the miRBase species based name format. For miRNA families observed to be present in both A-and B-genomes, paralogous miRNA loci counts in each Musa genome were estimated based on the 300 nt precursor regions predicted by miRDeep2. miRNA targets were predicted with 'psRNAtarget' online server (http://plantgrn. noble.org/psRNATarget/) [49] with default options.

Results and discussion
The majority of edible cultivated banana varieties are inter-and intra-specific triploid hybrids between varieties of two wild diploid M. acuminata and M. balbisiana species. To be able to carry out molecular genetic studies in these cultivars therefore it is necessary to have a reference B-genome. For these reasons we carried out Illumina 100 bp paired end gDNA sequencing of the wild diploid M. balbisiana variety 'Pisang Klutuk Wulang' (PKW), which is considered by many in the Musa research community to be one of the possible ancestors of the modern AAB hybrids.

PKW gDNA Read Mapping to the A-genome
Over 281 M trimmed Illumina gDNA reads were mapped against the 12 reference A-genome chromosomal sequences. An overview of the quality parameters of these reads is supplied in Additional file 1. The results following alignment show that 86.9% of all the PKW gDNA reads successfully mapped to the A-genome, with a mean read coverage per reference A chromosome, of 41.4x ( Table 1). The resulting consensus PKW genome sequences derived from these mappings were on average 78.9% of the length of the A-genome, with the greatest difference being found for chromosome 10 which was only 74.9% of the size of the corresponding 'Pahang' chromosome 10 ( Table 1). The results from chromosome 12 (B-chrUn_random) were not included in this calculation because this 'chromosome' consists of non-assembled contigs, concatenated together with spacer 'N's'. Therefore, while the assembled PKW consensus sequence for 'B-chrUn_random' derived from read mapping is 80.01 Mb, for molecular genetic studies, all gaps were filled with 'N's to prevent artificial concatenation and 'bridging' of mapped reads across unrelated genic regions, giving a working sequence length of 141.13 Mbp, and a working PKW consensus genome size of 402.5 Mbp.
This difference in size between the A-and B-genomes was largely expected as flow cytometric analyses of five M. acuminata genotypes (including DH 'Pahang') and representing 4 different subspecies, and of four M. balbisiana genotypes (including PKW) indicate that the haploid B-genome size is on average only 90% of the size of the A-genome [50][51][52]. However according to recent work from Cizkova et al. [50] the PKW haploid genome is actually 93.3% of the size of DH 'Pahang' with a sizes of 554 and 594 Mb respectively. This value for DH Pahang is slightly higher than the 523 Mb reported by D'Hont et al. 2012 [14]. The sequenced 'Pahang' genome size of 473 Mb, represents an assembly of~90% of the total M. acuminata genome, of which only 323 Mb is anchored to the 11 chromosomes [14]. If the expected PKW M. balbisiana genome size is 93.3% of 'Pahang' then we would expect to generate a consensus length of in the region of 0.93 x 473 or 440 Mb. By analogy therefore, our consensus PKW B-genome size of 341.4 Mb represents~78% of the expected PKW B-genome size.
De novo assembly of unmapped gDNA reads A total of 36.8 M gDNA reads (36.3 Gbp), remained unmapped after alignment to the A-genome. These could represent reads from regions that are structurally highly divergent from the A-genome so to test for the presence of unique, genic B-genome regions, the unmapped reads were de novo assembled into 63,245 contigs (Additional file 2: Table S2), and the presence of genic sequences tested for by large gap mapping of Musa unigene and reference CDS sequences, followed by a round of transcript detection. In total, 58,746 reads were used, but only 28 sequences actually mapped to these contigs. We can therefore conclude that the unmapped gDNA reads do not contain any significant gene rich regions, and that essentially all genic regions are retained in the consensus PKW B-genome sequence. An overview of the repeats annotation of these contig sequences is provided in Additional file 3: Table S3.

De novo assembly of gDNA reads
We also carried out de novo assembly of all gDNA reads, independent of a reference sequence. Here, over 96% of the 281 M trimmed reads, representing 27.4 Gbp of nucleotide sequence were assembled into 180,175 contigs with a total length of 339.3 Mb, an N50 of 7,884 bp, and an average contig length of 1,883 bp (Additional file 4: Table S4). The accumulated assembled contig length of 339.3 Mb is very similar to the consensus read mapping length of 341 Mb, but due to its much more fragmented nature (78.2% of these contigs were less than 1 kb in length) this resource is much more difficult to utilize. To evaluate the set of PKW gDNA contig sequences, the Musa reference CDS set was mapped to the PKW contig set as well to as the consensus PKW B-genome. In the case of the consensus PKW B-genome 32,192 (88% of total) Musa CDS were successively mapped, corresponding to 25,565 individual transcripts. In the case of the gDNA contig set, 71% of the CDS could be mapped (25,694 CDS), and a total of 21,272 individual transcripts were identified (data not shown). These data indicate therefore that simply mapping the gDNA reads to the Agenome and extracting the consensus sequence is the most effective way to generate a draft working M. balbisiana genome.

Evaluation/characterisation of the PKW B-genome assembly
A visual inspection of the gDNA mappings to the reference A-genome clearly demonstrates that there are many regions of structural variance between the two genomes ( Figure 1). However in general, the gene-rich regions seem to be well-conserved, as evidenced by the higher percentage of unbroken paired reads (in blue) in these regions.  regions by comparison typically contain a much higher proportion of unpaired, broken reads and more sequence variants.
While these results demonstrate that there are large regions with a high degree of homology across the M. acuminata and M. balbisiana genomes, it is important to realize that with this read mapping approach we cannot determine whether large scale genomic rearrangements such as insertions, inversions, transversions etc. have taken place. Previous work in banana however, suggests that gene order is likely to have been preserved -at least over short regions.  high degree of microsynteny with preservation of gene order, and 96 -96.3% sequence identity within genic regions [53,54]. The same authors also reported that predicted gene structure was good for well-conserved homologous genes, but that discrepancies were detected in the gene predictions of those orthologous BACs whose protein products had no match in public databases (i.e. hypothetical protein genes).

Variant detection in the PKW B-genome
In total, 20,657,389 sequence variants were detected in the PKW B-genome relative to the reference doubled haploid 'Pahang' A-genome based on only the uniquely mapped reads (Additional file 5: Table S5). Of these 18,868,899 were single nucleotide variants (SNVs or SNPs), 815,805 insertions and 972,588 deletions. From the list of SNVs 8,738,760 were homozygous variants which therefore represent sequence differences from 'Pahang'. The remaining 10,130,236 were heterozygous variants and therefore represent allelic variation and the degree of heterozygosity present in the PKW B-genome. On the basis of the total consensus PKW genome size of 341,431,243 bp, this heterozygosity corresponds to a SNP frequency of 1 variant every 33.7 bp, or 2.97%. The number and densities of sequence variants present in heterozygotic eukaryotic genomes varies enormously according to the species, whether they are obligate out-crossers or not, the number and genetic diversity of the cultivars assessed, and whether coding or noncoding regions are being considered. For example, in maize Zea mays L.) the SNV density was 1 SNV per 124 bp of coding sequence, and 1 per 31 bp in non-coding regions [55]. Visual inspection of Figure 1 also confirms that SNV density is higher in the intronic regions of the PKW B-genome than in the exons, and that there is also a higher SNV frequency in the non-transcribed regions. Similar conclusions were reached by Boonruangrod et al. in their comparison of the rDNA sequences of M. acuminata type 'Calcutta4' and 'Yangambi KM5' with the wild type M. balbisiana accession 'Tani [5].
By comparison, our data indicate that the PKW gDNA sequence differs from the A-genome at a frequency of 1 (homozygous) SNV every 39.1 bp (2.6%). This is much higher than estimates for interspecies SNP variation in rice, where a comparison of the O. indica and O. japonica genome sequence data found a SNP frequency of around 0.4%, [56]. However differences between M. acuminata and M. balbisiana are comparable to frequencies measured in Eucalyptus species (between one SNP per 16 bp to one per 33 bp were reported [57]), Quercus crispula (one SNP in 25 bp [58]) and Populus tremula (one SNP per 60 bp [59]). Estimates for the degree of interspecific variation in Musa are likely to grow as more cultivars and species are sequenced.

Repeat detection and annotation
In total, repetitive regions were found to occupy 26.85% (108.1 Mbp) of the PKW consensus B-genome, which is similar to the 27.76% reported by D'Hont et al. for the A-genome [1]. Annotation of the repetitive sequences of the B-genome showed that overall, the numbers of repeat elements is slightly higher in the B-genome and that the Ty1/copia and Ty3/Gypsy repeats dominate, representing 18.8% and 6.3% of the genome respectively (Table 2). Whilst the numbers of Non-LTR transposons (LINE), DNA transposons (clDNA and DNAhat) and Satellite repeats (Type 1 and Type 2) are similar in both A-and B-genomes and represent less than 1% of the total consensus B-genome sequence, the LTR transposons (Ty1/copia and Ty3/Gypsy) are more abundant in the B-genome.

Microsatellite detection
Microsatellite (SSR) sequences have many advantages as molecular markers, due to their abundance, hypervariability, co-dominant nature, reliability, and ease of interpretation and several of groups have already identified SSR markers for Musa [37,[60][61][62][63]. Type 1 SSRs (i.e. repeat sequences with a length greater than 20 bp), are considered to be hypervariable and the most efficient loci for use as molecular markers [64]. Analysis of the type I SSRs present in the A-and B-genomes demonstrates that the density of SSRs is slightly higher in the PKW B-genome at 1 per 323 bp, versus 1 per 387 bp in the 'Pahang' A-genome ( Table 3), but that the proportions of microsatellite motifs were similar in both species. A comparison of the microsatelites identified by Hribova et al. [37] in~100 Mb of repetitive sequences from M. acuminata cv. 'Calcutta 4' indicates a much higher proportion of trimeric and tetrameric motifs.
As shown in Table 3, dimeric repeat motifs are the most abundant type of class I SSR present in both Musa genomes, and TA the most common motif, representing 61.3 and 55.3% of all dimer repeats in the 'Pahang' and PKW genomes respectively ( Table 4). The next most abundant dimer repeats were GA/TC at 35.5 and 41.3% and TG/CA at 3.1 and 3.35% respectively. A previous study of 'Pahang' BAC end sequences by Arango et al., also found AT/TA to be the most prevalent SSR representing~26% of all SSR motifs [63], and characterization of the repeat component of~30% of the M. acuminata cv. 'Calcutta 4′ genome by low depth 454 sequencing also found TA and GA to be the most common dinucleotide repeats [37]. By comparison an extensive in silico study of EST databanks by Victoria et al. found AG/CT and GA/TC also to be the most common dimer motifs amongst vascular plant species [64]. The most abundant trimeric motifs in both M. acuminata and M. balbisiana genomes were ATT/AAT, ATA/TAT, GAA/TTC, and AAG/CTT, while Hribova et al. found GAA to be the most abundant trinucleotide motif in their study of the repetitive regions of the 'Calcutta 4' genome [37].

miRNA prediction
A non-redundant set of plant miRNAs, which included the 37 miRNA families (representing 234 precursors) previously reported for the 'Pahang' A-genome [14] was used to predict miRNA precursors and families within both Musa genomes (Figure 2). The results show a slightly larger number of predicted known miRNA precursors for the B-genome (270 miRNA precursors compared to 266 for the A genome), but the diversity of known miRNA families was lower, with 42 families predicted for the B-genome compared to 47 families for the A-genome. All of the known miRNA families detected in the B-genome were also found to be present in the A-genome. Overall, 10 additional miRNA families were found compared to those reported by D'Hont et al. [14]. The additional miRNA families detected were miR415, miR529, miR1134, miR5021, miRf10125, miRf10576, miRf11033, miRf11036, miRf11143 and miRf11357. Of these, only miR415 and miRf11036 were not detected in the B-genome. These new families may be due to additional entries added to the PMRD database since the analysis by D'Hont and colleagues, but could also be due to the new M. acuminata small RNA sequence data used in our query dataset. These were large (averaging >11 million Illumina sRNAseq "clean reads" per library) libraries, and derived from several different banana tissues (leaf, root, flower, fruit and somatic embryogenic cultures) -Additional file 6: Table S6.
As expected, the 42 miRNA families common to both genomes are shared with common ancestors of Musa i.e. embrophytes, angiosperms and poales. For example, the miR528 family which had previously been reported only for poales genomes, and recently shown by D'Hont et al. to be present in the A-genome [14] and is demonstrated here to also be present in the B-genome. Among the newly-predicted known miRNA families present in both Musa genomes, miR1134 has been reported as being abiotic stress -related and is found in the monocots Triticum aestivum [65] and Festuca arundinacea [66]; the miR5021 family has been reported to be pollen specific in Arabidopsis [67] whilst the families miRf10125, miRf10576, miRf11033, miRf11143 and miRf11357 are all of unknown function but were also computationally predicted from the Arabidopsis, poplar and rice genome sequences [68]. In addition to the known miRNA families, there were also 32 Musa miRNA precursors predicted, that belong to 28 novel miRNA families with no significant match to any previously reported mature miRNA sequence (Table 5). These include sequences that were unique to either Musa A-or B-genomes in addition to 4 families common to both genomes.

Genome distribution of miRNA precursors
Although the overall B-genome size is smaller, it contains more repetitive elements (RE's) as well as more known and novel miRNA (loci), but fewer conserved miRNA families than the A-genome. The three known miRNA families that were most highly represented in both Musa A-and B-genomes showed similar patterns of distribution across the chromosomes of both genomes, (Figure 2) suggesting synteny of A-and B-genomes. However, as the B-genome was assembled using the A-genome as reference, the gene order in this draft are preliminary and validation by FISH or similar methods will be needed to confirm this. As would be expected, the more recently evolved, novel miRNA families that are unique to either the M. acuminata or M. balbisiana genome are distributed fairly evenly across the genomes, with the exception of A-Chr2, B-Chr2, A-Chr3, B-Chr3, and A-Chr9 and A-Chr10, which lack any these sequences ( Figure 3). The higher number of miRNA loci in the B-genome may be related to a higher number of transposable elements (transposons and retrotransposons) present as these are thought to have contributed to the generation of species-specific miRNA genes in plants [69]. The differences in retroelements and miRNA on homeologous chromosomes suggest therefore that some of these miRNAs have arisen after the whole genome duplication events, since chromosomes 9 and 10 (Musa block 2 of D'Hont et al. 2012) are among the regions thought not to have been involved in the paleopolyploidisation [14].

Targets of novel B-genome miRNA
Using a fairly stringent cut-off expect value of 2.0, 18 predicted novel miRNA families were found in the B-genome. Of these, seven (mba-miR3, mba-miR5, mba-miR8, mba-miR12, mba-miR13, mba-miR15 and mba-miR18) were found to have targets within coding regions of the B-genome sequence. None of these predicted miRNA families were present in the A-genome [14] and thus are presumed to be B-genome specific in function and to have evolved after the divergence of the M. balbisiana and M. acuminata species~4.6 Mya [54]. Predicted targets of these B-genome candidate miRNAs correlate with a range of functions across plant development and metabolism and notably, several are proposed to be involved in tolerance or responses to biotic and abiotic stress. These included casein kinase, a predicted target of mba-miR3 and previously reported to affect multiple developmental and stress response pathways in Arabidopsis [70]; dirigent, a predicted target of mba-miR15, which are a family of proteins associated with lignification, biotic and abiotic stress responses, also recently reported as responding to drought, salt and oxidative stresses in sugarcane [71] and mba-miR8 targets from the multidrug and toxic compound extrusion (MATE) family, which in plants are associated with tolerance to various xenobiotic compounds [72] and metals including zinc tolerance in Arabidopsis [73] and aluminium tolerance in sorghum [74]. An additional predicted target of mba-miR8, Sal-1 phosphatase, has been reported as a negative regulator of drought tolerance in Arabidopsis, acting via ABA-dependent and independent pathways [75].  Results from Hribova et al. [37].
Several of the B-genome specific miRNA predicted targets, are proteins with unknown function, and appear to be only present in the B-genome and thus may represent B-genome distinct functional networks. Given the wider stress and disease resistance reported for banana B genomes [23,76,77] further functional validation of these miRNA and target genes, and in particular those of unknown function, is of particular interest. A full list of predicted miRNA targets is provided in Additional file 7: Table S7.

Annotating the PKW consensus B-genome
The results above demonstrate that the consensus genome derived from PKW gDNA mapping to the A-genome encompasses essentially all the genic B-genome regions. Analysis of the sequence variants and type 1 SSRs present in this consensus genome also produced results that show SSR and variant counts are comparable to those reported for the A-genome, and are in broad agreement with values obtained for other eukaryotic plant species. However, to be able to use the PKW genome sequences for molecular genetic research, gene annotation of the B-genome is obviously desirable. Therefore protein coding sequences in the B-genome were predicted by ab initio gene identification using the FGENESH software (http://linux1.softberry. com/all.htm). This resulted in the identification of 39,914 unique gene models. This number is higher than the 36,483 predicted from direct transfer of annotated regions of the A-genome to the B-genome and suggests that the gene count has been overestimated to the extent of around 109%. The higher predicted gene count is largely due to a nearly 2-fold higher number of gene models located within the concatenated contig set 'B_chrUn_random', relative to the chrUn_random of the A-genome ( Table 6). The set of PKW B-genome gene models were descriptively annotated online using the Blast2Go software [35,36]. Blast results against the NCBI non-redundant protein database show that 38,886 (97.4%) of the sequences had a positive hit, of which 30,541 had an e-value of 0. Following annotation steps, GO terms could be assigned to 37,367 (93.6%) sequences and 34,044 (85.3%) were annotated by interproscan. On the basis of the annotations assigned here, we can see that the B-genome contains 3,276 transposable elements (TE's), of which 1,470 are located in the B-unChr_random sequence. If these TS's are removed, the final functional B-genome gene count is actually 36,638, which is almost identical to the A-genome count of 36,542 of [14]. An overview and comparison of the A-and B-genomes and gene counts and annotations is provided in Table 6.

RNA read and transcript mapping to the A-and B-genomes
To evaluate the usefulness and validity of the consensus PKW B-genome for expression studies, a total of 256 M paired-end, 100 bp RNA reads from two AAA, and three AAB cultivars, and one diploid Austalimusa Fe'i cultivar were mapped to the combined A-and B-genome sets. In addition we carried out de novo assembly of these reads to generate transcript sets for each cultivar, and these expressed gene transcripts were also mapped to the combined genomes set together with a 'Grande Naine' (AAA), Unigene set.
As expected, the proportion of mapped 100 bp reads generally decreased with the predicted genetic distance of the cultivar from the reference A-genome, with 'Karat' (Fe'i) having the lowest value (74.8%) and 'Yangambi- km5' (AAA) having the highest proportion at 90.7% (Additional file 8: Table S8). Although reads from the AAA cultivars should theoretically map only to the A-genome, we can see that 23 -27% of the reads from the two triploid AAA cultivars preferentially map (have a higher homology) to regions of the B-genome, despite the apparent absence of a B-genome in these cultivars (Figure 4, Additional file 9: Table S9). This could reflect the different sub-group origins of these M. acuminata genomes. i.e. 'Yangambi-km5', and 'Gros Michel' belong to the 'Ibota' and 'Gros Michel' M. acuminata subgroups respectively, whereas the reference 'Pahang' belongs to the Malaccencis subgroup. Also, differences between orthologous genes in these M. acuminata subgroup genomes could mean that only minor sequence divergence from 'Pahang' could lead to higher homology to the orthologous sequences present in the B-genome, particularly for highly-conserved 'core' genes. Alternatively the presence of foreign chromosomal fragments as a result of historical recombinations between the A-and B-genomes as demonstrated by Jeridi et al. [13], could result in the mapping of reads/transcripts to homologous regions of the homeologous B-chromosomes. By comparison, in the AAB plantains, we see that 36.4 -40.7% of the reads from the three plantains preferentially map to the B-genome. This is in accordance with the presence of a single B-genome in these triploids, and confirms the utility of our PKW consensus B-genome sequence for this type of study. Finally, reads from the diploid Fe'i banana cultivar (Musa, Australimusa, 2n = 20), and a species which is probably most closely related to the wild species M. maclayi, M. peekelii and M. lolodensis [78], mapped nearly equally to both genomes (48. 6 : 51.4, A B).
Interestingly, the normalized read coverage across the all chromosomes of both genomes was also found to differ between the cultivars. Of particular note are the differences between the two AAA hybrids (Figure 4, Additional file 9: Table S9), where 40.2% of all 'Yangambi-km5' mapped reads localized to the A-and B-chrUn_ random' sequences, compared to only 16.7% of the 'Gros Michel' reads. This suggests that there could be substantial differences between the AAA genome sequences of 'Yangambi km-5' and 'Pahang'. However we also see that 'B-chrUn_random' is also the sequence with the highest count of TE sequences, so that these differences could simply represent differences in the abundance of these elements. Indeed comparing the genes with the highest expression levels in this chromosome in these 2 cultivars, shows that 250,822 reads specifically map to a single 472 bp intronic region of sequence (KMMuB_chrUn_random_G39317)_in 'Yangambi km-5', while 'only' 60,891 map to this same sequence in 'Gros Michel'. No such phenomenon was observed elsewhere however, where we see that a much higher proportion of 'Gros Michel' reads map to A_chr3 (Figure 4), and where for all cultivars a proportionally higher number of reads mapped to A_chr6, which is also the chromosome with the highest gene count (3,794, Table 6). However, even though the homeologous B_chr6 gene also has the highest gene count of the B-genome, more reads from the African plantains ('Mbouroukou' and 'Batard'), mapped to B_chr7. These results indicate that there are different patterns of expression across the two chromosome sets (Figure 4), and suggest different and unique contributions of each genome to banana metabolism. The bias in gene coverage across the chromosomes of both genomes was examined further by looking at the mapping of the de novo assembled RNA contigs (transcripts) derived from the same RNA reads. Here the longer mean read length of the sequences improves the specificity of mapping and allows us to make a comparison with the results of the A-genome CDS and unigene mappings (Table 7, Figure 5).
Firstly, we see that 99.8% of the reference 'Pahang' CDS map to the A-genome, as expected. Similarly, 96.1% of the 'Yangambi-km5' (AAA) transcripts, and 94.1% of the 'unigene' sequences derived from "Grande Naine" (AAA, Cavendish sub-group) map to the B-genome (Table 6)   Looking at the mean normalized expressed transcript coverage per chromosome, it is clear that the high proportion of 100 bp reads mapping to the 'chrUn_random' sequences ( Figure 4) is no longer evident. This suggests that the majority of the differences between cultivars involves non-expressed regions. For all cultivars the highest proportion of mapped reads are to be found on chromosomes A_chr6, and A_chr8, and the lowest number of mapped transcripts occurs on the two chrUn_random sequences. The African plantains are interesting in having a proportionally higher number of reads mapping to B_chr7 despite the low gene count of this chromosome. Some care has to be taken in the interpretation of differences in read distributions across the chromosomes as variations in physiological status (e.g. ripening) could lead to shifts in the mean mapped read depths per chromosomes. Nonetheless these results could provide indications as to those chromosomes in which there has been a higher degree of homeologous recombination, or which are preferentially expressed under certain conditions. The difference in the distribution of mapped between M. acuminata and M. balbisiana diploid progenitor species. Therefore the study of these hybrids ideally requires the presence of reference sequences for both genomes. In particular for the introgression of beneficial resistance traits associated with the M. balbisiana genome. Here, we generated a consensus Bgenome sequence, which based on size estimates, repre-sents~78% of the expected B-genome length. This smaller size is presumably due to differences in the organisation of intergenic repetitive regions between the A-and B-genomes. This B-genome shows a high degree of sequence divergence from the A-genome as well as moderately high levels of heterozygosity. Ab initio gene prediction for the B-genome generates a total of 39,914 gene models, of which 3,540 are annotated as TE's, so that the number of functional gene sequences is actually 36,638, which is nearly identical to the gene count of 'Pahang'. We also identified several new B-genome specific miRNAs, some of which have predicted targets suggestive of novel stress related pathways that may have evolved separately in M. Balbisiana. To validate the usefulness of the consensus B-genome, mRNA-reads and de novo assembled transcripts from a series of genetically diverse interspecific triploid cultivars were mapped to the combined A-and B-genomes set. Results suggest significant genetic divergence between subgroups of AAA cultivars, and the presence of regions with high B-genome homology causing up to 12.7% of the transcripts to map preferentially to the B-genome. Comparison of the transcript mappings of the African and Plantain (AAB) varieties suggests structural diversity between these two groups, but importantly validate the quality of the B- genome, as~33% of all reads map exclusively to these Bgenome sequences. For all varieties there is evidence for differences in expression levels across homeologous chromosomes suggesting independent contributions of the two genomes to banana metabolism. While our results demonstrate the usefulness of the consensus B-genome sequence for molecular genetic studies, it is important to relalise that the structure and organization is heavily dependent on the quality of the A-genome used as a reference. In this respect it should be remembered that only~70% of the published A-genome sequence has actually been anchored to the 11 chromosomes. As a result, large scale structural differences in the B-genome are difficult to detect with our strategy, and intergenic regions of high diversity, and regions of the B-genome consisting of repeats will be difficult to assemble. In particular the position of the TEs, that represent around half of the Musa genome will be based on their position in the A-genome while in reality their position will in the vast majority of cases be different in the B genome.

Conclusions
There is much interest in the exploitation of the M. balbisiana genome sequence for the introgression of beneficial traits such as biotic and abiotic stress resistance within Musa breeding programs. The PKW B-genome sequence therefore represents a valuable resource for the molecular genetic studies, not only in AAB, and ABB cultivars, but also for AAA dessert bananas. In addition we have shown that it can be used obtain useful information on