Three New Genome Assemblies Support a Rapid Radiation in Musa acuminata (Wild Banana).

Edible bananas result from interspecific hybridization between Musa acuminata and Musa balbisiana, as well as among subspecies in M. acuminata. Four particular M. acuminata subspecies have been proposed as the main contributors of edible bananas, all of which radiated in a short period of time in southeastern Asia. Clarifying the evolution of these lineages at a whole-genome scale is therefore an important step toward understanding the domestication and diversification of this crop. This study reports the de novo genome assembly and gene annotation of a representative genotype from three different subspecies of M. acuminata. These data are combined with the previously published genome of the fourth subspecies to investigate phylogenetic relationships. Analyses of shared and unique gene families reveal that the four subspecies are quite homogenous, with a core genome representing at least 50% of all genes and very few M. acuminata species-specific gene families. Multiple alignments indicate high sequence identity between homologous single copy-genes, supporting the close relationships of these lineages. Interestingly, phylogenomic analyses demonstrate high levels of gene tree discordance, due to both incomplete lineage sorting and introgression. This pattern suggests rapid radiation within Musa acuminata subspecies that occurred after the divergence with M. balbisiana. Introgression between M. a. ssp. malaccensis and M. a. ssp. burmannica was detected across the genome, though multiple approaches to resolve the subspecies tree converged on the same topology. To support evolutionary and functional analyses, we introduce the PanMusa database, which enables researchers to exploration of individual gene families and trees.

Introduction

Bananas are among the most important staple crops cultivated worldwide in both the tropics and subtropics. The wild ancestors of bananas are native to the Malesian Region (including Malaysia and Indonesia) (Simmonds 1962) or to northern Indo-Burma (southwest China). Dating back to the early Eocene (Janssens et al. 2016), the genus Musa currently comprises 60–70 species divided into two sections, Musa and Callimusa (Häkkinen 2013). Most of modern cultivated bananas originated from natural hybridization between two species from the section Musa, Musa acuminata, which occurs throughout the whole southeast Asia region, and Musa balbisiana, which is constrained to an area going from east India to south China (Simmonds and Shepherd 1955). While no subspecies have been defined so far in M. balbisiana, M. acuminata is further divided into multiple subspecies, among which at least four have been identified as contributors to the cultivated banana varieties, namely banksii, zebrina, burmannica, and malaccensis (reviewed in Perrier et al. 2011). These subspecies can be found in geographical areas that are mostly nonoverlapping. Musa acuminata ssp. banksii is endemic to New Guinea. Musaa. ssp. zebrina is found in Indonesia (Java island), M. a. ssp. malaccensis originally came from the Malay Peninsula (De Langhe et al. 2009; Perrier et al. 2011), while M. a. ssp. burmannica is from Burma (today’s Myanmar) (Cheesman 1948).

While there are many morphological characters that differentiate M. acuminata from M. balbisiana, the subspecies of M. acuminata have only a few morphological differences between them. For instance, M. a. ssp. burmannica is distinguished by its yellowish and waxless foliage, light brown markings on the pseudostem, and by its compact pendulous bunch and strongly imbricated purple bracts. Musaa. ssp. banksii exhibits slightly waxy leaf, predominantly brown-blackish pseudostems, large bunches with splayed fruits, and nonimbricated yellow bracts. Musaa. ssp. malaccensis is strongly waxy with a horizontal bunch, and bright red nonimbricated bracts, while M. a. ssp. zebrina is characterized by dark red patches on its dark green leaves (Simmonds 1956).

Previous studies based on a limited number of markers have been able to shed some light on the relationships among M. acuminata subspecies (Sardos et al. 2016; Christelová et al. 2017). Phylogenetic studies have been assisted by the availability of the reference genome sequence for a representative of M. acuminata ssp. malaccensis (D’Hont et al. 2012; Martin et al. 2016) and a draft M. balbisiana genome sequence (Davey et al. 2013). However, the availability of large genomic data sets from multiple (sub)species are expected to improve the resolution of phylogenetic analyses, and thus to provide additional insights on species evolution and their specific traits (Bravo et al. 2018). This is especially true in groups where different segments of the genome have different evolutionary histories, as has been found in Musaceae (Christelová et al. 2011). Whole-genome analyses also make it much easier to distinguish among the possible causes of gene tree heterogeneity, especially incomplete lineage sorting (ILS) and hybridization (Folk et al. 2018).

Moreover, the availability of multiple reference genome sequences opens the way to so-called pangenome analyses, a concept coined by Tettelin et al. (2005). The pangenome is defined as the set of all gene families found among a set of phylogenetic lineages. It includes 1) the core genome, which is the pool of genes common to all lineages, 2) the accessory genome, composed of genes absent in some lineages, and 3) the species-specific or individual-specific genome, formed by genes that are present in only a single lineage. Identifying specific compartments of the pangenome (such as the accessory genome) offers a way to detect important genetic differences that underlie molecular diversity and phenotypic variation (Morgante et al. 2007).

Here, we generated three de novo genomes for the subspecies banksii, zebrina and burmannica, and combined these with existing genomes for M. acuminata ssp. malaccensis (D’Hont et al. 2012) and M. balbisiana (Davey et al. 2013). We thus analyzed the whole genome sequences of five extant genotypes comprising the four cultivated bananas’ contributors from M. acuminata, that is, the reference genome “DH Pahang” belonging to M. acuminata ssp. malaccensis, “Banksii” from M. acuminata ssp. banksii, “Maia Oa” belonging to M. acuminata ssp. zebrina, and “Calcutta 4” from M. acuminata ssp. burmannica, as well as M. balbisiana (i.e., “Pisang Klutuk Wulung” or PKW). We carried out phylogenomic analyses that provided evolutionary insights into both the relationships and genomic changes among lineages in this clade. Finally, we developed a banana species-specific database to support the larger community interested in crop improvement.

Materials and Methods

Plant Material

Banana leaf samples from accessions “Banksii” (Musa acuminata ssp. banksii, PT-BA-00024), “Maia Oa” (Musa acuminata ssp. zebrina, PT-BA-00182), and “Calcutta 4” (Musa acuminata ssp. burmannica, PT-BA-00051) were supplied by the CRB-Plantes Tropicales Antilles CIRAD-INRA field collection based in Guadeloupe. Leaves were used for DNA extraction. Plant identity was verified at the subspecies level using SSR markers at the Musa Genotyping Centre (MGC, Czech Republic) as described in (Christelová et al. 2011) and passport data of the plant is accessible in the Musa Germplasm Information System (Ruas et al. 2017). In addition, the representativeness of the genotypes of the four subspecies was verified on a set of 22 samples belonging to the same four M. acuminata subspecies of the study (supplementary fig. 3, Supplementary Material online).

Sequencing and Assembly

Genomic DNA was extracted using a modified MATAB method (Risterucci et al. 2000). DNA libraries were constructed and sequenced using the HiSeq2000 (Illumina) technology at BGI (supplementary table 1, Supplementary Material online). “Banksii” was assembled using SoapDenovo (Luo et al. 2012), and PBJelly2 (English et al. 2012) was used for gap closing using PacBio data generated at the Norwegian Sequencing Center (NSC) with Pacific Biosciences RS II. “Maia Oa” and “Calcutta 4” were assembled using the MaSuRCA assembler (Zimin et al. 2013) (supplementary table 2, Supplementary Material online). Estimation of genome assembly completeness was assessed with BUSCO plant (Simão et al. 2015) (supplementary table 3, Supplementary Material online).

Gene Annotation

Gene annotation was performed on the obtained de novo assembly for “Banksii,” “Maia Oa,” and “Calcutta 4,” as well as on the draft Musa balbisiana “PKW” assembly (Davey et al. 2013) for consistency and because the published annotation was assessed as low quality. For structural annotation we used EuGene v4.2 (http://eugene.toulouse.inra.fr/) (Foissac et al. 2008) calibrated on M. acuminata malaccensis “DH Pahang” reference genome v2, which produced similar results (e.g., number of genes, no missed loci, good specificity, and sensitivity) as the official annotation (Martin et al. 2016). EuGene combined genotype-specific (or closely related) transcriptome assemblies, performed with Trinity v2.4 with RNAseq data sets (Sarah et al. 2017), to maximize the likelihood to have genotype-specific gene annotation (supplementary table 4, Supplementary Material online). The estimation of gene space completeness was assessed with Busco (supplementary table 3, Supplementary Material online). Because of its high quality and to avoid confusing the community, we did not perform a new annotation for the M. a. malaccensis “DH Pahang” reference genome but used the released version 2. Finally, the functional annotation of plant genomes was performed by assigning their associated generic GO terms through the Blast2GO program (Conesa et al. 2005) combining BLAST results from UniProt (E-value 1e-5) (Magrane and UniProt Consortium 2011).

Gene Families

Gene families were identified using OrthoFinder v1.1.4 (Emms and Kelly 2015) with default parameters based on BLASTp (e-value 1e-5). Venn diagrams were made using JVenn online (http://jvenn.toulouse.inra.fr) (Bardou et al. 2014) and alternate visualization was produced with UpsetR (https://gehlenborglab.shinyapps.io/upsetr) (Lex et al. 2014).

Tree Topology from Literature

A species tree was initially identified based on previous studies (Janssens et al. 2016; Sardos et al. 2016). Those two studies included all M. acuminata subspecies, and had the same tree topology (supplementary fig. 1, Supplementary Material online). In the first study, Sardos et al. (2016) computed a Neighbor-Joining tree from a dissimilarity matrix using biallelic GBS-derived SNP markers along the 11 chromosomes of the Musa reference genome. Several representatives of each subspecies that comprised genebank accessions related to the genotypes used here were included (Sardos et al. 2016). We annotated the tree to highlight the branches relevant to M. acuminata subspecies (supplementary fig. 2, Supplementary Material online). In the second study, a maximum clade credibility tree of Musaceae was proposed based on four gene markers (rps16, atpB-rbcL, trnL-F, and internal transcribed spacer, ITS) analyzed with Bayesian methods (Janssens et al. 2016).

Genome-Scale Phylogenetic Analyses and Species Tree

Single-copy OGs (i.e., orthogroups with one copy of a gene in each of the five genotypes) from protein, coding DNA sequence (CDS), and genes (including introns and UTRs) were aligned with MAFFT v7.271 (Katoh and Standley 2013), and gene trees were constructed using PhyML v3.1 (Guindon et al. 2009) with ALrT branch support. All trees were rooted using Musa balbisiana as outgroup using Newick utilities v1.6 (Junier and Zdobnov 2010). Individual gene tree topologies were visualized as a cloudogram with DensiTree v2.2.5 (Bouckaert 2010).

Single-copy OGs were further investigated with the quartet method implemented in ASTRAL v5.5.6 (Mirarab and Warnow 2015; Zhang et al. 2018). In parallel, we carried out a Supertree approach following the SSIMUL procedure (http://www.atgc-montpellier.fr/ssimul/) (Scornavacca et al. 2011) combined with PhySIC_IST (http://www.atgc-montpellier.fr/physic_ist/) (Scornavacca et al. 2008) applied to a set of rooted trees corresponding to core OGs (including single and multiple copies), and accessory genes for which only one representative species was missing (except outgroup species). Finally, single-copy OGs (CDS only) were used to generate a concatenated genome-scale alignment with FASconCAT-G (Kück and Longo 2014) and a tree was constructed using PhyML (NNI, HKY85, 100 bootstrap).

Search for Introgression

Ancient gene flow was assessed with the ABBA-BABA test or D-statistic (Green et al. 2010; Durand et al. 2011) and computed on the concatenated multiple alignment converted to the MVF format and processed with MVFtools (Pease and Rosenzweig 2018), similar to what is described in Wu et al. (2017) (https://github.com/wum5/JaltPhylo). The direction of introgression was further assessed with the D2 test (Hibbins and Hahn 2018). The D2 statistic captures differences in the heights of genealogies produced by introgression occurring in alternate directions by measuring the average divergence between species A and C in gene trees with an ((A, B), C) topology (denoted [dAC|A, B]), and subtracting the average A–C divergence in gene trees with a ((B, C), A) topology (denoted [dAC|B, C]), so that D2 = (dAC|A, B)−(dAC|B, C). If the statistic is significantly positive, it means that introgression has either occurred in the B→C direction or in both directions. D2 significance was assessed by permuting labels on gene trees 1,000 times and calculating p values from the resulting null distribution of D2 values. The test was implemented with a Perl script using distmat from EMBOSS (Rice et al. 2000) with Tajima–Nei distance applied to multiple alignments associated with gene trees fitting the defined topologies above (https://github.com/mrouard/perl-script-utils).

Results

Assemblies and Gene Annotation

We generated three de novo assemblies belonging to M. acuminata ssp. banksii, M. a. ssp. zebrina, and M. a. ssp. burmannica. The M. a. ssp. zebrina and M. a. ssp. burmannica assemblies contained 56,481 and 47,753 scaffolds (N50 scaffold of 37,689 and 22,183 bp) totaling 623 Mb and 526 Mb, respectively. The M. a. ssp. banksii assembly, which benefited from long-read sequencing, contained 9,467 scaffolds (N50 scaffold of 435,833 bp) for a total of 464 Mb (78.2% of the genome) (supplementary tables 1 and 2, Supplementary Material online).

The number of predicted protein coding genes per genome within different genomes of Musa ranges from 32,692 to 45,069 (supplementary table 3, Supplementary Material online). Gene number was similar for M. a. ssp. malaccensis “DH Pahang,” M. balbisiana “PKW,” and M. a. ssp. banksii “Banksii” but higher in M. a. ssp. zebrina “Maia Oa” and M. a. ssp. burmannica “Calcutta 4.” According to BUSCO (supplementary table 4, Supplementary Material online), the most complete gene annotations are “DH Pahang” (96.5%), “Calcutta 4” (74.2%), and “Banksii” (72.5%), followed by “PKW” (66.5%) and “Maia Oa” (61.2%).

The percentage of genes in orthogroups (OGs), which is a set of orthologs and recent paralogs (i.e., gene family), ranges from 74 in M. a. zebrina “Maia Oa” to 89.3 in M. a. malaccensis “DH Pahang” with an average of 79.8 (table 1). Orthogroups have a median size of 4 genes and do not exceed 50 (supplementary table 5, Supplementary Material online). A pangenome here was defined on the basis of the analysis of OGs in order to define the 1) core, 2) accessory, and 3) unique gene set(s). On the basis of the five genomes studied here, the pangenome embeds a total of 32,372 OGs composed of 155,222 genes. The core genome is composed of 12,916 OGs (fig. 1). Among these, 8,030 are composed of only one sequence in each lineage (i.e., are likely single-copy orthologs). A set of 1,489 OGs are specific to all subspecies in M. acuminata, while the number of genes specific to each subspecies ranged from 14 in the M. acuminata “DH Pahang” to 110 in M. acuminata “Banksii” for a total of 272 genes across all genotypes. No significant enrichment for any Gene Ontology (GO) category was detected for subspecies-specific OGs.

Table 1

Summary of the Gene Clustering Statistics Per (Sub)Species

	Musa acuminata malaccensis “DH Pahang”	M. acuminata burmannica “Calcutta 4”	M. acuminata banksii “Banksii”	M. acuminata zebrina “Maia Oa”	M. balbisiana “PKW”
# genes	35,276	45,069	32,692	44,702	36,836
# genes in orthogroups	31,501	34,947	26,490	33,059	29,225
# unassigned genes	3,775	10,122	6,202	11,643	7,611
% genes in orthogroups	89.3	77.5	81	74	79.3
% unassigned genes	10.7	22.5	19	26	20.7
# orthogroups containing species	24,074	26,542	21,446	25,730	23,935
% orthogroups containing species	74.4	82	66.2	79.5	73.9
# species-specific orthogroups	6	46	47	11	9
# genes in species-specific orthogroups	14	104	110	23	21
% genes in species-specific orthogroups	0	0.2	0.3	0.1	0.1

. 1.

—Intersection diagram showing the distribution of shared gene families (at least two sequences per OG) among M. a. banksii “Banksii,” M. a. zebrina “Maia Oa,” M. a. burmannica “Calcutta 4,” M. a. malaccensis “DH Pahang,” and M. balbisiana “PKW” genomes. The figure was created with UpsetR (Lex et al. 2014).

Variation in Gene Tree Topologies

Phylogenetic reconstruction performed with single-copy genes (n = 8,030) showed high levels of discordance among the different individual gene trees obtained, both at the nucleic acid and protein levels (fig. 2 and supplementary data 1, Supplementary Material online). Considering M. balbisiana as outgroup, there are 15 possible bifurcating tree topologies relating the four M. acuminata subspecies. For all three partitions of the data—protein, CDS, and gene (including introns and UTRs)—we observed all 15 different topologies (table 2). We also examined topologies at loci that had bootstrap support >90 for all nodes, also finding all 15 different topologies (table 2). Among trees constructed from whole genes, topologies ranged in frequency from 13.12% for the most common tree to 1.92% for the least common tree (table 2) with an average length of the 1,342 aligned nucleotide sites for CDS and 483 aligned sites for proteins. Based on these results, gene tree frequencies were used to calculate concordance factors on the most frequent CDS gene trees (table 2), demonstrating that no split was supported by >30% of gene trees (fig. 2). Therefore, in order to further gain insight into the subspecies phylogeny, we used a combination of different approaches described in the next section.

. 2.

—Illustration of gene tree discordance. (A) Cloudogram of single copy OGs (CDS) visualized with Densitree. The blue line represents the consensus tree as provided by Densitree. (B) Species tree with bootstrap-like support based on corresponding gene tree frequency from table 2 (denoted topology number 2). PKW, M. balbisiana “PKW”; C4, M. acuminata burmannica “Calcutta 4”; M, M. acuminata zebrina “Maia Oa”; DH, M. acuminata malaccensis “DH Pahang”; B, M. acuminata banksii “Banksii”.

Table 2

Frequency of Gene Tree Topologies of the 8,030 Single Copy OGs

No.	Topology	# CDS (%)	# Protein (%)	# Gene (%)	# Gene Bootstrap >90 (%)
1	(PKW,(C4,(M,(DH, B))))	11.9	10.58	13.12	13.72
2	(PKW,(C4,(DH,(B, M))))	10.8	10.48	11.92	14.88
3	(PKW,((DH, C4),(B, M)))	9.59	7.28	12.73	17.52
4	(PKW,(M,(C4,(DH, B))))	9.53	12.51	7.78	5.91
5	(PKW,(C4,(B,(DH, M))))	8.02	7.37	8.89	8.44
6	(PKW,((DH, B),(C4, M)))	7.67	6.55	9.16	12.56
7	(PKW,(M,(B,(DH, C4))))	6.66	8.21	5	3.06
8	(PKW,(B,(M,(DH, C4))))	5.58	5.23	4.61	2.53
9	(PKW,(DH,(C4,(B, M))))	5.41	5.21	5.18	4.96
10	(PKW,(B,(C4,(DH, M))))	5.26	4.45	6.2	7.07
11	(PKW,(B,(DH,(C4, M))))	5.02	6.82	3.36	1.9
12	(PKW,(M,(DH,(B, C4))))	4.23	4.68	2.84	1.16
13	(PKW,((DH, M),(B, C4)))	4.037	3.61	4.79	5.06
14	(PKW,(DH,(B,(C4, M))))	3.85	4.18	2.44	0.63
15	(PKW,(DH,(M,(B, C4))))	2.38	2.77	1.92	0.52

Note.—In bold, the most frequent topology.

PKW, Musa balbisiana “PKW”; C4, Musa acuminata burmannica “Calcutta 4”; M, Musa acuminata zebrina “Maia Oa”; DH, Musa acuminata malaccensis “DH Pahang”; B, Musa acuminata banksii “Banksii”.

Inference of a Species Tree

We used three complementary methods to infer phylogenetic relationships among the sampled lineages. First, we concatenated nucleotide sequences from all single-copy genes (totaling 11,668,507 bp). We used PHYML to compute a maximum likelihood tree from this alignment, which, as expected, provided a topology with highly supported nodes (fig. 3). Note that this topology (denoted topology number 1 in table 2) is not the same as the one previously proposed in the literature (denoted topology number 7 in table 2) (supplementary figs. 1 and 2, Supplementary Material online).

. 3.

—Species topologies computed with three different approaches. (A) Maximum likelihood tree inferred from a concatenated alignment of single-copy genes (CDS). (B) Supertree-based method applied to single and multilabelled gene trees. (C) Quartet-based model applied to protein, CDS, and gene alignments.

Next, we used a method explicitly based on individual gene tree topologies. ASTRAL (Mirarab and Warnow 2015) infers the species tree by using quartet frequencies found in gene trees. It is suitable for large data sets and was highlighted as one of the best methods to address challenging topologies with short internal branches and high levels of discordance (Shi and Yang 2018). ASTRAL found the same topology using ML gene trees from single-copy genes obtained from protein sequences, CDSs, and genes (fig. 3).

Finally, we ran a supertree approach implemented in PhySIC_IST (Scornavacca et al. 2008) on the single-copy genes and obtained again the same topology (fig. 3). PhySIC_IST first collapses poorly supported branches of the gene trees into polytomies, as well as conflicting branches of the gene trees that are only present in a small minority of the trees; it then searches for the most resolved supertree that does not contradict the signal present in the gene trees nor contains topological signal absent from those trees. Deeper investigation of the results revealed that ∼ 66% of the trees were unresolved, 33% discarded (pruned or incorrectly rooted), and therefore that the inference relied on fewer than 1% of the trees. Aiming to increase the number of genes used by PhySIC_IST, we included multicopy OGs of the core genome, as well as some OGs in the accessory genomes using the pipeline SSIMUL (Scornavacca et al. 2011). SSIMUL translates multilabeled gene trees (MUL-trees) into trees having a single copy of each gene (X-trees), that is, the type of tree usually expected in supertree inference. To do so, all individual gene trees were constructed on CDSs from OGs with at least 4 M. acuminata and M. balbisiana genes (n = 18,069). SSIMUL first removed identical subtrees resulting from a duplication node in these trees, it then filtered out trees where duplicated parts induced contradictory rooted triples, keeping only coherent trees. These trees can then be turned into trees containing a single copy of each gene, either by pruning the smallest subtrees under each duplication node (leaving only orthologous nodes in the tree), or by extracting the topological signal induced by orthology nodes into a rooted triplet set, that is then turned back into an equivalent X-tree. Here, we chose to use the pruning method to generate a data set to be further analyzed with PhySIC_IST, which lead to a subset of 14,507 gene trees representing 44% of the total number of OGs and an increase of 80% compared with the 8,030 single-copy OGs. This analysis returned a consensus gene tree with the same topology as both of the previous methods used here (fig. 3).

Evidence for Introgression

Although much of the discordance we observe is likely due to incomplete lineage sorting, we also searched for introgression between subspecies. A common approach, performed in other plant genomes (Eaton and Ree 2013; Eaton et al. 2015; Novikova et al. 2016; Choi et al. 2017), relies on the use of the ABBA-BABA test (or D statistics) (Green et al. 2010). This test allows to differentiate admixture from incomplete lineage sorting across genomes by detecting an excess of either ABBA or BABA sites (where “A” corresponds to the ancestral allele and “B” corresponds to the derived allele state). An excess of each of this patterns is indicative of ancient admixture. Therefore, we applied it in a four-taxon phylogeny including three M. acuminata subspecies as ingroups and M. balbisiana as outgroup. Because there were five taxa to be tested, analyses were done with permutation of taxa denoted P1, P2, and P3 and Outgroup (table 3). Under the null hypothesis of ILS, an equal number of ABBA and BABA sites are expected. However, we always found an excess of sites grouping malaccensis (“DH”) and burmannica (“C4”) (table 3). This indicates a history of introgression between these two lineages.

Table 3

Four-Taxon ABBA-BABA Test (D-Statistic) Used for Introgression Inference from the Well-Supported Topology from Fig. 3

P1	P2	P3	BBAA	ABBA	BABA	Disca	Db	p valuec
Malaccensis (DH)	Banksii (B)	Burmannica (C4)	12185	4289	8532	0.51	−0.33	<2.2e-16
Malaccensis (DH)	Zebrina (M)	Burmannica (C4)	9622	5400	9241	0.6	−0.26	< 2.2e-16
Zebrina (M)	Banksii (B)	Burmannica (C4)	11204	6859	6782	0.54	0.005	0.5097
Malaccensis (DH)	Banksii (B)	Zebrina (M)	10450	7119	6965	0.57	0.02	0.1944

Discordance=(ABBA+BABA)/Total

D =(ABBA−BABA)/(ABBA+BABA)

Based on Pearson chi-squared.

To test the direction of introgression, we applied the D2 test (Hibbins and Hahn 2018). While introgression between a pair of species (e.g., malaccensis and burmannica) always results in smaller genetic distances between them, the D2 test is based on the idea that gene flow in the two alternative directions can also result in a change in genetic distance to other taxa not involved in the exchange (in this case, banksii). We computed the genetic distance between banksii and burmannica in gene trees where malaccensis and banksii are sister (denoted dAC|A, B) and the genetic distance between banksii and burmannica in gene trees where malaccensis and burmannica are sister (denoted dAC|B, C). The test takes into account the genetic distance between the species not involved in the introgression (banksii) and the species involved in introgression that it is not most closely related to (burmannica). We identified 1,454 and 281 gene trees with dAC|A, B = 1.15 and dAC|B, C = 0.91, respectively, giving a significant positive value of D2=0.23 (p < 0.001 by permutation). These results support introgression from malaccensis into burmannica, though they do not exclude the presence of a lesser level of gene flow in the other direction.

PanMusa, a Database to Explore Individual OGs

Since genes underlie traits and wild banana species showed a high level of incongruent gene tree topologies, access to a repertoire of individual gene trees is important. This was the rationale for constructing a database that provides access to gene families and individual gene family trees in M. acuminata and M. balbisiana. A set of web interfaces are available to navigate OGs that have been functionally annotated using GreenPhyl comparative genomics database (Rouard et al. 2011). PanMusa shares most of the features available on GreenPhyl to display or export sequences, InterPro assignments, sequence alignments, and gene trees (fig. 4). In addition, new visualization tools were implemented, such as MSAViewer (Yachdav et al. 2016) and PhyD3 (Kreft et al. 2017) to view gene trees.

. 4.

—Overview of available interfaces for the PanMusa database. (A) Homepage of the website. (B) List of functionally annotated OGs. (C) Graphical representation of the number of sequence by species. (D) Consensus InterPro domain schema by OG. (E) Individual gene trees visualized with PhyD3. (F) Multiple alignment of OG with MSAviewer.

Discussion

Musa acuminata Subspecies Contain Few Subspecies-Specific Families

In this study, we used a de novo approach to generate additional reference genomes for the three subspecies of Musa acuminata; all three are thought to have played significant roles as genetic contributors to the modern cultivars. Genome assemblies produced for this study differ in quality, but the estimation of genome assembly and gene annotation quality conducted with BUSCO suggests that they were sufficient to perform comparative analyses. Moreover, we observed that the number of genes grouped in OGs were relatively similar among subspecies, indicating that the potential overprediction of genes in “Maia Oa” and “Calcutta 4” was mitigated during the clustering procedure. Indeed, overprediction in draft genomes is expected due to fragmentation, leading to an artefactual increase in the number of genes (Denton et al. 2014).

Although our study is based on one representative per subspecies, Musa appears to have a widely shared pangenome, with only a small number of subspecies-specific families identified. The pangenome analysis also reveals a large number of families shared only among subsets of species or subspecies (fig. 1); this “dispensable” genome is thought to contribute to diversity and adaptation (Tettelin et al. 2005; Medini et al. 2005). The small number of species-specific OGs in Musa acuminata also supports the recent divergence between all genotypes including the split between M. acuminata and M. balbisiana.

Musa acuminata Subspecies Show a High Level of Discordance between Individual Gene Trees

Gene tree conflict has been recently reported in the Zingiberales (Carlsen et al. 2018) and Musa in not an exception. By computing gene trees with all single-copy genes OG, we found widespread discordance in gene tree topologies. Topological incongruence can be the result of incomplete lineage sorting, the misassignment of paralogs as orthologs, introgression, or horizontal gene transfer (Maddison 1997). With the continued generation of phylogenomic data sets over the past dozen years, massive amounts of discordance have been reported, first in Drosophila (Pollard et al. 2006) and more recently in birds (Jarvis et al. 2014), mammals (Li et al. 2016; Shi and Yang 2018), and plants (Novikova et al. 2016; Pease et al. 2016; Choi et al. 2017; Copetti et al. 2017; Wu et al. 2017). Due to the risk of hemiplasy in such data sets (Avise et al. 2008; Hahn and Nakhleh 2016), we determined that we could not accurately reconstruct either nucleotide substitutions or gene gains and losses among the genomes analyzed here.

In our case, the fact that all possible subspecies tree topologies occurred, and that ratios of minor trees at most nodes were equivalent to those expected under ILS, strongly suggests the presence of ILS (Hahn and Nakhleh 2016). Banana is a paleopolyploid plant that experienced three independent whole genome duplications (WGD), and some fractionation is likely still occurring (D’Hont et al. 2012) (supplementary table 6, Supplementary Material online). But divergence levels among the single-copy OGs were fairly consistent (fig. 2), supporting the correct assignment of orthology among sequences. However, we did find evidence for introgression between malaccensis and burmannica, which contributed a small excess of sites supporting one particular discordant topology (table 3). This event is also supported by the geographical overlap in the distribution of these two subspecies (Perrier et al. 2011).

Previous studies have attempted to resolve the topology in the Musaceae, but did not include all subspecies considered here, and had very limited numbers of loci. In Christelova et al. (2011), a robust combined approach using maximum likelihood, maximum parsimony, and Bayesian inference was applied to 19 loci, but only burmannica and zebrina out of the four subspecies were included. Jarret et al. (1992) reported sister relationships between malaccensis and banksii on the basis of RFLP markers, but did not include any samples from burmannica and zebrina. However, the resolved species tree supported by all methods used here is a new topology compared with species trees comprising at least one representative of our 4 subspecies (Janssens et al. 2016; Sardos et al. 2016; Christelová et al. 2017) (supplementary fig. 1, Supplementary Material online). Indeed “Calcutta 4” as representative of M. acuminata ssp. burmannica was placed sister to the other Musa acuminata genotypes in our study, whereas those studies indicates direct proximity between burmannica and malaccensis. The detected introgression from malaccensis to burmannica may be an explanation for the difference observed but increasing the sampling with several genome sequences by subspecies would enable a better resolution.

More strikingly considering previous phylogenetic hypotheses, malaccensis appeared most closely related to banksii, which is quite distinct from the other M. acuminata spp. (Simmonds and Weatherup 1990) and which used to be postulated as its own species based on its geographical area of distribution and floral diversity (Argent 1976) (fig. 5). However, on the bases of genomic similarity, all our analyses support M. acuminata ssp. banksii as a subspecies of M. acuminata.

. 5.

—Area of distribution of Musa species in Southeast Asia as described by Perrier et al. (2011); including species tree of Musa acuminata subspecies based on results described in figure 4. Areas of distribution are approximately represented by colors; hatched zone shows area of overlap between two subspecies where introgression may have occurred.

Gene Tree Discordance Supports Rapid Radiation of Musa acuminata Subspecies

In their evolutionary history, Musa species dispersed from “northwest to southeast” into Southeast Asia (Janssens et al. 2016). Due to sea level fluctuations, Malesia (including the nations of Indonesia, Malaysia, Brunei, Singapore, the Philippines, and Papua New Guinea) is a complex geographic region, formed as the result of multiple fusions and subsequent isolation of different islands (Thomas et al. 2012; Janssens et al. 2016). Ancestors of the Callimusa section (of the Musa genus) started to radiate from the northern Indo-Burma region toward the rest of Southeast Asia ∼30 Ma, while the ancestors of the Musa (formerly Eumusa/Rhodochlamys) section started to colonize the region ∼10 Ma (Janssens et al. 2016). The divergence between M. acuminata and M. balbisiana has been estimated to be ∼5 Ma (Lescot et al. 2008). However, no accurate dating has yet been proposed for the divergence of the Musa acuminata subspecies. We hypothesize that after the speciation of M. acuminata and M. balbisiana (ca. 5 Ma) rapid diversification occurred within populations of M. acuminata. This hypothesis is consistent with the observed gene tree discordance and high levels of ILS. Such a degree of discordance may reflect a near-instantaneous radiation between all subspecies of M. acuminata. Alternatively, it could support the proposed hypothesis of divergence back in the northern part of Malesia during the Pliocene (Janssens et al. 2016), followed by introgression taking place among multiple pairs of species as detected between malaccensis and burmannica. While massive amounts of introgression can certainly mask the history of lineage splitting (Fontaine et al. 2015), we did not find evidence for such mixing.

Interestingly, such a broad range of gene tree topologies due to ILS (and introgression) has also been observed in gibbons (Carbone et al. 2014; Veeramah et al. 2015; Shi and Yang 2018) for which the area of distribution in tropical forests of Southeast Asia is actually overlapping the center of origin of wild bananas. Moreover, according to Carbone et al. (2014), gibbons also experienced a near-instantaneous radiation ∼5 Ma. It is therefore tempting to hypothesize that ancestors of wild bananas and ancestors of gibbons faced similar geographical isolation and had to colonize and adapt to similar ecological niches, leading to the observed patterns of incomplete lineage sorting.

In this study, we highlighted the phylogenetic complexity in a genome-wide data set for Musa acuminata and Musa balbisiana, bringing additional insights to explain why the Musaceae phylogeny has remained controversial. Our work should enable researchers to make inferences about trait evolution, and ultimately should help support crop improvement strategies.

Click here for additional data file.