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

Plant materials

Sterile plantlets of the wild diploid M. balbisiana variety ‘Pisang Klutuk Wulang’ (‘PKW’ , BB genome), were obtained from the International Transit Centre, KU Leuven (reference number ITC 1587), and healthy leaf material collected for gDNA isolation from 8-week old, greenhouse grown plants. RNA was extracted from lyophilised fruit pulp samples from the following cultivars; ‘Yangambi-km5’ (AAA), ‘Gros Michel’ (AAA), ‘Batard’ (AAB), ‘Mbouroukou-3’ (AAB), ‘Iholena Lele sub var 'Long Peduncular'’ (AAB), and ‘Karat’ (Fe’i, Australimusa).

Banana root tissues were prepared from clonal tissue-cultured 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.

PKW gDNA Read Mapping to the A-genome

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–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 A-genome 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. For example, direct transfer of annotations from the A-genome to the new PKW B-genome results in the transfer of 36,483 gene sequences (99.8%), indicating that regions homologous to essentially all genic regions of the A-genome are present in the PKW B-genome. Intergenic/non-transcribed 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. For example, a comparison of 1.4 Mb of orthologous BAC clones from M. acuminata (cv. ‘Madang’), and M. balbisiana (cv. ‘Pisang Klutuk Wulung’), showed a 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

Microsatellite detection

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.

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].

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

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).

Sequence Name	Reference length/bp	Normalised mapped transcripts count
		CDS	Unigenes	Ygi	GM	Mbk	Btd	Iho	Karat
A_chr1	27,573,629	76.5	75.7	79.1	71.4	54.6	52.9	54.9	47.0
A_chr2	22,054,697	63.7	59.4	57.4	53.5	40.9	40.5	39.7	30.0
A_chr3	30,470,407	89.6	85.6	90.7	81.6	62.8	60.7	61.9	50.1
A_chr4	30,051,516	92.5	90.0	88.6	81.9	63.1	61.7	57.8	45.3
A_chr5	29,377,369	81.6	75.6	76.9	69.7	53.9	56.3	53.3	42.4
A_chr6	34,899,179	100.7	103.2	102.2	92.0	73.1	70.4	71.0	54.8
A_chr7	28,617,404	75.6	75.1	74.9	68.7	24.5	25.3	49.8	37.8
A_chr8	35,439,739	94.1	90.3	98.8	86.5	69.4	68.0	68.7	51.8
A_chr9	34,148,863	84.6	79.3	80.6	76.2	60.0	57.3	64.2	43.6
A_chr10	33,665,772	86.6	88.2	83.9	75.9	68.7	67.9	60.5	42.7
A_chr11	25,514,024	73.1	67.1	70.0	64.6	49.0	47.3	50.9	34.7
A_chrUn_random	141,147,818	79.7	51.6	58.1	51.3	37.7	37.8	50.4	31.7
B_chr1	22,038,404	0.3	4.4	3.1	10.4	29.6	32.1	32.7	42.1
B_chr2	17,349,238	0.1	3.8	2.3	7.6	22.4	24.1	23.2	29.6
B_chr3	24,161,952	0.3	4.8	3.7	11.8	33.5	35.4	34.4	47.6
B_chr4	24,656,528	0.1	4.9	3.5	11.1	28.5	29.2	35.1	43.4
B_chr5	23,648,591	0.1	5.2	3.0	10.6	27.2	28.0	29.7	41.7
B_chr6	27,831,592	0.2	5.9	4.4	13.9	36.4	36.2	36.1	47.6
B_chr7	22,212,853	0.1	5.9	3.9	10.2	61.5	59.6	28.8	39.2
B_chr8	27,665,716	0.1	6.6	2.8	12.2	25.6	26.3	24.6	46.3
B_chr9	25,900,723	0.1	5.9	3.8	10.8	21.8	23.0	4.7	41.4
B_chr10	25,230,959	0.1	4.1	3.8	11.6	10.0	10.8	22.2	43.5
B_chr11	20,721,546	0.1	4.8	2.5	9.0	24.3	26.9	26.5	35.6
B_chrUn_random	141,147,818	0.2	2.7	2.1	7.6	21.6	22.1	18.9	30.2
% A	54.0	99.8	94.1	96.1	87.3	65.8	64.6	68.3	51.2
% B	46.0	0.2	5.9	3.9	12.7	34.2	35.4	31.7	48.8

Read counts normalised per million total reads to allow direct comparison between cultivars. CDSm ‘Pahang’, AA variety CDS reference set; Unigenes’, EST set derived from AAA cultivar ‘Grande Naine’; Ygi, AAA cultivar ‘Yangambi-km5’; GM, AAA cultivar ‘Gros Michel’; Mbk-3, AAB cultivar ‘Mbouroukou-3’; Btd, AAB cultivar ‘Batard’; Iho, AAB cultivar ‘Iholena lele’; Karat, diploid Fe’i cultivar ‘Karat’.

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). Somewhat surprisingly though, only 87.3% of the ‘Gros Michel’ (AAA) transcripts map to the A-genome suggesting either that the ‘Ibota’ subgroup is more closely related to ‘Malacenssis’ subgroup to which ‘Pahang’ belongs than to the ‘Gros Michel’ subgroup, or that there is a larger proportion of M. balbisiana-like gene sequences in the ‘Gros Michel’ genome. Importantly the proportion of AAB transcripts mapping to the B-genome approaches the theoretically expected 33%, at 34.2 – 35.4% for the African plantains, but is slightly different than the 31.5% observed for Pacific plantain cultivar ‘Iholena lele’. Despite the fact that ‘Karat’ belongs to a completely different species (Australimusa), diploid ‘Karat’ transcripts map to both reference genomes equally well (51.2:48.8, A:B respectively), indicating that these resources could be useful to study gene expression even in Australimusa cultivars.

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 ‘Iholena lele’ (Pacific plantain) transcripts compared to the African plantains, ‘Mbouroukou-3’ and ‘Batard’ supports their different classifications. i.e. even though all are hybrids between the M. acuminata subsp. banksii and M. balbisiana genomes [79], ‘Iholena lele’ is classified as a member of the Iholena and Maoli-Popoulu subgroup of plantains while the African plantains are considered to have arisen as a result of somaclonal mutations following a limited number of introductions to the African continent 2,300 – 2,750 years ago.

Summary

The majority of modern cultivated banana cultivars are interspecific sterile triploids formed by hybridization 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 B-genome sequence, which based on size estimates, represents ~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 B-genome 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.

Data accession

The raw Illumina PKW gDNA sequence data generated was submitted to SRA with accession no. SRR956987. The assembled PKW consensus chromosome sequences, with annotations, as well as the CDS set and list of gDNA contigs are hosted at the Banana Genome Hub and publically available at http://banana-genome.cirad.fr/content/musa-balbisiana-pisang-klutuk--wulung.