Assembly and comparative analyses of the mitochondrial genome of Castanospermum australe (Papilionoideae, Leguminosae)

Additional keywords: genome skimming, legume, mitochondrion, Moreton Bay chestnut, phylogeny, read-pair mapping.

In contrast to the explosion of plastid and animal mitochondrial genome sequencing over the past 15 years, the pace of plant mitochondrial genome sequencing has been relatively slow owing to its extreme structural variation (Smith and Keeling 2015). Likewise, fewer studies in plant phylogenetics utilise mitochondrial genes because of their low substitution rate (Hiesel et al. 1994). More than 100 angiosperm mitochondrial genomes (http://www.ncbi.nlm.nih.gov/genome/organelle/, accessed 10 March 2019) have been sequenced, and show remarkable variation in genome size, content and structure (Knoop et al. 2011; Mower et al. 2012a; Guo et al. 2016). Sequenced angiosperm mitochondrial genomes range from 66 kb in Viscum scurruloideum Barlow (Skippington et al. 2015) to more than 11.3 Mb in Silene L. (Sloan et al. 2012). The almost 200-fold difference of mitochondrial genome size has been primarily attributed to the expansion or contraction of intergenic sequences. In angiosperms, pseudogene formation in mitochondrial genes as a result of functional transfer to the nucleus is frequent (Adams et al. 2002). The structure of most angiosperm mitochondrial genomes is extremely variable. This results from frequent homologous recombination mediated by repeats (Maréchal and Brisson 2010; Gualberto and Newton 2017), which can make proper assembly of mitochondrial genomes very challenging with short-read sequences. Recent studies have leveraged the increasing availability of paired-end libraries to complete assemblies of mitochondrial genomes (e.g. Naito et al. 2013; Guo et al. 2016, 2017; Yu et al. 2018). Here, we employ these techniques to sequence and assemble the mitochondrial genome of Castanospermum australe, an economically important legume species.

Only a few Leguminosae mitochondrial genomes have been reported, most included within subfamily Papilionoideae, which is the largest subfamily of Leguminosae, with ~503 genera and 14 000 described species (Legume Phylogeny Working Group 2017). Mitochondrial genome sizes in Papilionoideae range from 217 kb in Medicago truncatula Gaertn. (Bi et al. 2016) to 588 kb in Vicia faba L. (Negruk 2013). These mitochondrial genomes share 30 protein-coding genes and three rRNA genes (Shi et al. 2018), as well as the loss of three ribosomal protein genes (rps2, rps11 and rps13). Several other protein-coding genes (including cox2, rpl2, rps1, rps19, sdh3 and sdh4) are functional in some Papilionoideae mitochondrial genomes, but pseudogenised or lost in others, with a functional copy having been transferred to the nuclear genome (Palmer et al. 2000; Shi et al. 2018). Papilionoideae mitochondrial genomes also demonstrate great variation in the number and length of repeat sequences. Some mitochondrial genomes have high numbers of long repeats, such as in Glycine max (L.) Merr., which has many large repeats that have mediated recombination to produce an enriched molecular pool of 760 circles (Chang et al. 2013). Likewise, Vicia faba has 15 repeats >1 kb in length, with the longest repeat of 66 kb. In contrast, a few legume mitochondrial genomes have only short repeats of <1 kb in length, such as Vigna radiata (L.) R.Wilczek, which has only six repeats of >100 bp (Alverson et al. 2011a).

Castanospermum australe, the Moreton Bay chestnut or black bean, is the only species of the genus Castanospermum A.Cunn ex Hook. This species is native to north-eastern Australia, but has been introduced as an ornamental tree into India, South Africa and warm temperate regions of North America. The seeds of C. australe yield castanospermine, which is a major toxic alkaloidal constituent that inhibits replication of the human immunodeficiency virus (HIV; Walker et al. 1987) and other retroviruses (Sunkara et al. 1987), and reduces tumour growth in mice (Ostrander et al. 1988). Castanospermum australe belongs to the Aldinoid clade, one of the early branching lineages of Papilionoideae (Legume Phylogeny Working Group 2013). Mitochondrial genomes of a few major clades of Papilionoideae have been sequenced, including the Cladrastis clade, the Genistoid clade, the Hologalegina clade, and the Millettioid clade. However, no mitochondrial genome of the Aldinoid clade has been reported.

The purpose of the present study was to sequence and assemble the mitochondrial genome of C. australe and to compare it to other Papilionoideae mitochondrial genomes to help clarify the diversification of mitochondrial genomes in this subfamily. Given the challenges of assembling mitochondrial genomes, it is important to verify alternative conformations mediated by repeats. To accomplish this, we employed a computational approach that was developed to use read-pair counts to estimate alternative conformations (e.g. Alverson et al. 2011b; Mower et al. 2012a). Across Papilionoideae, we undertook comparative analyses of mitochondrial genome size, gene content, structural variation and RNA editing patterns. We also tested the utility of mitochondrial genomes in resolving the deep relationships of the subfamily.

Fresh leaves of C. australe were collected from the Brisbane Botanic Garden, Queensland, Australia. The voucher specimen (Yi14471) was deposited in the herbarium of the Kunming Institute of Botany, Chinese Academy of Sciences (KUN). No specific permits were required for the collection of the material. Total DNA was isolated from silica gel-dried leaves of C. australe by using a modified CTAB protocol described in Yang et al. (2014). The DNA sample was sequenced from a paired-end library with an insert size of 350 bp at Beijing Genomics Institute (Shenzhen, China) on an Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA), generating ~8.6 million 150-bp paired-end reads. Reads were deposited in the NCBI Sequence Read Archive (Accession SRR8648141).

Genome-sequence data were assembled de novo using Linux-OS SPAdes genome assembler (ver. 3.10.1, see http://cab.spbu.ru/software/spades/; Bankevich et al. 2012), with the following parameters: careful; k-mer values: 75, 89, 105, 113, 121; all other parameters were set to defaults. The use of iterative k-mer lengths in SPAdes leverages the full potential of paired-end reads. Assembled contigs were then filtered and binned into mitochondrial and plastid contigs by using the Python script ‘slim_fastg_by_blast’ (see https://github.com/Kinggerm/GetOrganelle, accessed August 2016). To avoid assembly errors at regions of shared homology among the different genomes, we carefully rechecked and deleted the plastid and nuclear DNA contigs, as described preciously (Mower et al. 2012a). The final hypothesised conformations (Fig. 1, more details in Results and discussion) were manually exported using Bandage Ubuntu dynamic (ver. 8.0, see http://rrwick.github.io/Bandage/; Wick et al. 2015).

Fig. 1.  Read-mapping strategy for confirming alternative mitochondrial genomic conformations (C1−C9). OC represents the map of the de novo assembly graph. Boxes of the same colour denote homologous repeats R1, R2, and R3. A, B, C, D, E, and F represent the fragments separated by repeats. Grey arrows in the circles denote the clockwise and anticlockwise orientation of fragment D. C1−C3 and C4−C8 show hypothesised subgenomic circles, and C4 and C9 show hypothesised completely assembled circular genomes. Dashed arrows indicate the rejection of hypothesised assemblies on the basis of paired-end read mapping. Solid arrows denote the two supported mitochondrial genome assemblies (MC1 and MC2).

To evaluate sequence-assembly quality and accuracy, paired-end reads were mapped to the assembled mitochondrial genomes using Bowtie2 (ver. 2.3.2, see http://bowtie-bio.sourceforge.net/bowtie2/index.shtml; Langmead and Salzberg 2012) with default parameters. Read maps were visually inspected in Geneious (ver. 9.1.14, see https://www.geneious.com/; Kearse et al. 2012). Seven subgenomic and two complete genomic conformations were hypothesised on the basis of three medium-size repeats (Fig. 1). Sequencing reads were set as pairs and mapped to these nine hypothetical structures. The following parameters were used for Bowtie2: end-to-end alignment type, medium sensitivity assembly, and a maximum insert size of 600 bp.

Genes, introns and tRNAs were annotated manually by using the previously published Vigna radiata (NC_015121), Glycine max (NC_020455) and Arabidopsis thaliana (L.) Heynh. (NC_001284) mitochondrial genomes as query sequences in Geneious, employing a basic local alignment search tool (BLAST)-like algorithm to search for annotations with ≥80% sequence similarity. Start and stop codons and intron splice sites were visually inspected and adjusted as necessary. A mitochondrial genome map was generated using OGDraw (see https://chlorobox.mpimp-golm.mpg.de/OGDraw.html; Lohse et al. 2013). The final annotated genome was deposited in GenBank (Accession MK426679).

The mitochondrial genome structure of C. australe was compared with all available Papilionoideae mitochondrial genomes, including those of Ammopiptanthus mongolicus (Kom.) S.H.Cheng (MF683210), Glycine max (NC_020455), Lotus japonicus (Regel) Larsen (NC_016743), Medicago truncatula (NC_029641), Pongamia pinnata (L.) Pierre (NC_016742), Styphnolobium japonicum (L.) Schott. (MG757109), Vicia faba (KC189947), Vigna radiata (NC_015121) and Vigna angularis (Willd.) Ohwi & H.Ohashi (NC_021092). Syntenic blocks shared by C. australe and other species were identified using both Mauve (ver. 2.3.1, see http://darlinglab.org/mauve/mauve.html; Darling et al. 2004), with default parameters as implemented in Geneious, as well as using BLASTn (ver. 2.8.1, see ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/) with a word size of seven and an e-value cutoff of 1 × 10⁻⁶, following Guo et al. (2016).

C-to-U RNA editing sites in the 30 mitochondrial protein coding genes of C. australe and other selected Papilionoideae mitochondrial genomes were predicted using the PREP-Mt online tool (see http://prep.unl.edu/; Mower 2009), with a cutoff value of 0.5, with the exception of six genes (cox2, rpl2, rps1, rps19, sdh3 and sdh4) that were not universally present in Papilionoideae mitochondrial genomes.

To test the utility of mitochondrial genome sequences in resolving the deep relationships of Papilionoideae, phylogenetic analyses were performed using 33 protein-coding genes (atp1, atp4, atp6, atp8, atp9, ccmB, ccmC, ccmFc, ccmFn, cob, cox1, cox2, cox3, matR, mttB, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7, nad9, rpl5, rpl16, rps1, rps3, rps4, rps10, rps12, rps14 and sdh4) of the 10 complete mitochondrial genomes of Papilionoideae, with the mitochondrial genome of Senna tora (L.) Roxb. (Caesalpinioideae; NC_038053) as the outgroup. Genes were aligned using MAFFT (ver. 7.1.2, see https://mafft.cbrc.jp/alignment/software/; Katoh and Standley 2013) using the L-INS-I algorithm. Gene alignments were subsequently adjusted manually and concatenated in Geneious, which produced a 35 360-bp nucleotide dataset. Phylogenetic analyses were performed with standard Bayesian inference (BI) and Maximum likelihood (ML) methods. The best-fit substitution model was estimated in PartitionFinder2 (ver. 2.1.1, see http://www.robertlanfear.com/partitionfinder/; Lanfear et al. 2017), by using the following settings: the all model, RAxML (ver. 8.2.12, see https://github.com/stamatak/standard-RAxML; Stamatakis 2014), the rcluster algorithm (Lanfear et al. 2014) with the rcluster-percent set to 10, and the Akaike information criterion (AICc) to compare the fit of the different models. The BI analysis was performed with MrBayes (ver. 3.2, see http://nbisweden.github.io/MrBayes/download.html; Ronquist et al. 2012), with four Markov Chain Monte Carlo (MCMC) runs using a random starting tree, an invgamma rate model with six discrete categories and 10 million generations, with a sampling frequency of 1 every 1000 generations, and the first 25% being discarded as burn-in. Stationarity was considered to be reached when the average standard deviation of split frequencies was <0.01. Maximum likelihood analysis was performed with RAxML. Branch support was estimated using 1000 rapid bootstrap replicates (‘-f a’ option) with the GTRGAMMA substitution model (Stamatakis 2006).

Repeats are likely to mediate intramolecular recombination within mitochondrial genomes, leading to their often-complex genome structure (Palmer and Shields 1984; Alverson et al. 2011b; Cole et al. 2018; Kovar et al. 2018). Such complexity often precludes confident mitochondrial genome assemblage using short inserts and reads, leading to the availability of few angiosperm mitochondrial genomes, including in Papilionoideae. For example, the large repeats of Glycine max are responsible for its complex mitochondrial genome structure (Chang et al. 2013). In sequencing the C. australe mitochondrial genome, we employed paired-end sequencing with a 350-bp insert size to try to maximise our ability to fully assemble the genome by using short (150 bp) reads. Only three repeats >100 bp in length were identified in C. australe, namely, R1 (119 bp), R2 (133 bp) and R3 (114 bp) (Fig. 1). Homologous recombination among these repeats yielded nine hypothetical alternative mitochondrial-genome arrangements (C1−C9; Fig. 1) on the basis of the original assembly graph (OC). Through mapping of paired-end reads, subgenomic conformations C1−C3 and C5−C8 were rejected (Fig. 1) because paired-end reads did not completely span these repeats. In contrast, two circular conformations (C4 and C9) that contained two copies of all three repeats and contained exactly the same genomic information were strongly supported. Conformations C4 and C9 differed only in the orientation of the 64-kb fragment D (Fig. 1). Read mapping against these two confirmations yielded 367 consistent paired-end reads that spanned Repeat R3 (Fig. 1) at the inversion breakpoint in the major conformation (C4, or MC1; Fig. 1), whereas only six paired-end reads consistently covered R3 at the inversion breakpoint in the minor conformation (C9, or MC2; Fig. 1). MC1 and MC2 are depicted as circular chromosomes with lengths of 542 079 bp. A circular map of MC1 is shown in Fig. 2; the sequence of MC1 was deposited in GenBank as accession MK426679. A linear map of both MC1 and MC2 is shown in Fig. 3. Repeats R1, R2 and R3 separate the MC1 and MC2 into six fragments, labelled A (303 789 bp), B (72 549 bp), C (2518 bp), D (64 174 bp), E (64 831 bp) and F (33 486 bp) (Fig. 1). The endpoints of fragment D, which is inverted between MC1 and MC2, are located between trnfM and Exon 1 of nad1 (Fig. 3).

Fig. 2.  Mitochondrial genome map (Master circular 1, MC1) of Castanospermum australe. Genes shown on the outside of the circle are transcribed clockwise, whereas those on the inside are transcribed counter-clockwise. Genes are colour-coded by functional groups; exons are shown as coloured boxes and introns are shown as white boxes. Asterisks indicate the presence of intron(s) in a gene. The map in the inner circle indicates the GC content of the mitochondrial genome.

Fig. 3.  Linear map of MC1 and MC2 in Castanospermum australe. Purple boxes indicate the 64-kb inversion between the two conformations.

Of the 8 610 483 paired-end reads, 704 104 mapped to the MC1 mitochondrial genome. The mean mitochondrial genome coverage was ~389.7× (Fig. S1, available as Supplementary material to this paper; coverage ranged from ~12× to ~911×). The GC content was 45.3%, which is typical for other selected mitochondrial genomes (Table S1, available as Supplementary material to this paper). The mitochondrial genome of C. australe is large among published Papilionoideae mitochondrial genomes, which a range from 217 618 bp in Medicago truncatula to 588 000 bp in Vicia faba (Table S1). The large size of the mitochondrial genome of C. australe is mainly due to the accumulation of species-specific non-coding sequences, which account for 87.7% (475 542 bp) of the mitochondrial genome. The length of the intergenic regions also contributes greatly to the variation in the mitochondrial genome size among other Papilionoideae species.

The mitochondrial genome of C. australe comprises 37 known protein-coding genes, three rRNAs and 18 tRNAs (Table 1). Of these, trnC–GCA and trnN are present in two copies. Of the protein-coding genes, 33 are also found in all other available Papilionoideae mitochondrial genomes (Table 1). The cox2 gene is absent from Vigna radiata and V. angularis, and rps1 is absent from Lotus japonicus. The full rRNA and tRNA gene set in C. australe is also shared with all other Papilionoideae species. The total length of the mitochondrial protein-coding gene sequences in C. australe is 30 135 bp, which is comparable to that in most other papilionoid species, but is 11 kb shorter than those of Ammopiptanthus mongolicus and Vicia faba. This is because the mitochondrial genome of A. mongolicus has seven chloroplast-derived protein-coding genes (atpA, atpB, ndhJ, ndhK, psbA, rbcL and rpoB), and the mitochondrial genome of V. faba possesses two or three copies of 10 mitochondrial genes (atp1, atp6, atp9, ccmC, ccmFc, mttB, nad7, nad9, rpl16 and rps3).

Table 1.  Comparisons of protein-coding gene content across Papilionoideae mitochondrial genomes
The 33 genes include atp [1, 4, 6, 8, 9], ccm [B, C, Fc, Fn], cob, matR, mttB, nad [1, 2, 3, 4, 4L, 5, 6, 7, 9], rpl [5, 16], rps [3, 4, 10, 12, 14], rrn [5, 18, 26] and cox [1, 3]. Y, gene present; N, gene absent; PS, putative pseudogene present. Species abbreviations: Cas, Castanospermum australe; Sty, Styphnolobium japonicum; Amm, Ammopiptanthus mongolicus; Lot, Lotus japonicus; Vic, Vicia faba; Med, Medicago truncatula; Pon, Pongamia pinnata; Gly, Glycine max; Vig_rad, Vigna radiata and Vig_ang, Vigna angularis

Intronic regions comprise 5.4% of the C. australe mitochondrial genome, including nine cis-spliced introns and eight trans-spliced introns (Table S1) across eight protein-coding genes (ccmFc, nad1, nad2, nad4, nad5, nad7, rps3 and rps10). These introns are homologous to those in other Papilionoideae mitochondrial genomes. Variation in total intron length among Papilionoideae (ranging from 27 066 bp in Lotus japonicus to 32 553 bp in Glycine max) was not high compared with variation in mitochondrial genome length. Five introns (ccmFc intron, nad2 Intron 2, nad4 Intron 3, rps10 intron and rps3 intron) showed high intron-length variations across Papilionoideae (Table S1, Fig. S2, available as Supplementary material to this paper); introns within the ccmFc gene of C. australe, Styphnolobium japonicum, Ammopiptanthus mongolicus, L. japonicus, Medicago truncatula and Vicia faba are only ~950 bp, but they are ~4 kb in other sequenced mitochondrial genomes of Papilionoideae. The intron of the rps10 gene also shows high length variation in Papilionoideae; the intron is 1.9 kb shorter in the Hologalegina clade (L. japonicus, Vicia faba and M. truncatula) than in remaining Papilionoideae. The length variations of ccmFc intron and rps10 intron are similar to those found in studies of Chang et al. (2013) and Shi et al. (2018).

Plant mitochondrial genomes have a high incidence of C-to-U RNA editing (Hiesel et al. 1989; Takenaka et al. 2008; Alverson et al. 2010; Suzuki et al. 2013), with ~200−800 editing sites in the protein-coding mitochondrial genes of most angiosperms (Mower 2008; Sloan et al. 2010; Richardson et al. 2013; Guo et al. 2016). However, the incidence of RNA editing across lineages of Papilionoideae is largely unknown (Kovar et al. 2018; Shi et al. 2018). Within C. australe, 487 C to U RNA editing sites were predicted among the 30 protein-coding genes (Table S2, available as Supplementary material to this paper). Among these RNA editing sites, 34.5% (168 sites) occurred in the first base position of the codon, 65.5% (319 sites) occurred in the second base position, and none occurred in the third position. Among these protein-coding genes, nad4 possessed the highest number of RNA editing sites (42 sites) and atp1 possessed the fewest (1 sites). In total, 52.8% (257 sites) of the RNA editing sites involved amino acids that were converted from hydrophilic to hydrophobic, whereas 7.6% (37 sites) involved conversion from hydrophobic to hydrophilic amino acids. In contrast, 30.6% (149 sites) involved hydrophobic to hydrophobic amino acid conversions and 8.2% (40 sites) involved hydrophilic to hydrophilic amino acid conversions. Two glycines and two arginines were converted into stop codons. The ccmB, nad4L and ccmC genes had the highest editing frequency, with 5.31, 4.29 and 4.05 edits per 100 nt respectively, whereas the editing frequencies of other genes were lower than four editing sites per 100 nt (Table S3, available as Supplementary material to this paper).

Papilionoideae mitochondrial genes showed RNA editing patterns that were largely similar to those of other angiosperms (Alverson et al. 2010). The positions and amino acid conversions of predicted RNA editing of the 10 mitochondrial genomes in Papilionoideae were provided in Tables S4–S13 (available as Supplementary material to this paper). Across 30 genes that were present in all 10 Papilionoideae mitochondrial genomes, the lowest levels of RNA editing were predicted in atp1 and atp9, with only one or two editing sites. The ndh4 gene had the highest level of editing, with an average of 40.8 editing sites. Variation in the RNA editing-site number across Papilionoideae was highest in atp6, where it ranged from two to 20, whereas four genes (atp4, atp8, rps14, rps4) had the same number of editing sites across all Papilionoideae mitochondrial genomes (Fig. 4, Tables S2, S–S13). Papilionoid ribosomal protein-coding genes had fewer edits, whereas ccmB, nad4L and ccmC had the highest editing frequency.

Fig. 4.  Boxplot diagram of RNA editing-site number for 30 genes showing the ranges of RNA editing-site number in genes across 10 Papilionoideae mitogenomes. Genes with missing data are not included.

Previous studies have shown that closely related taxa generally share more RNA-editing sites (Chen et al. 2011), and we detected a similar phenomenon in Papilionoideae; of the respective 487 and 489 RNA editing sites that were found in the C. australe and Styphnolobium japonicum mitochondrial genes, 473 sites were shared. In Vigna radiata and V. angularis, all 475 detected RNA editing sites were shared. In the four species of the Millettioid clade (Fig. 4), 455 sites were shared between Pongamia pinnata and Glycine max, 456 sites were shared between G. max and Vigna radiata or V. angularis, and 445 sites were shared between P. pinnata and V. radiata or V. angularis.

Analyses based on 33 mitochondrial genes resolved the deep phylogenetic relationships of Papilionoideae, with maximum bootstrap support and posterior probability at all nodes (Fig. 5). The Aldinoid clade and Cladrastis clade were successively sister to remaining Papilionoideae. These early diverging Papilionoideae relationships were not resolved in Legume Phylogeny Working Group (2017) on the basis of plastid matK sequence, although recent phylogenomic analyses using 81 plastid-coding genes recovered the same relationships (R. Zhang, Y. H. Wang, J. J. Jin, A. Bruneau, D. Cardoso, L. P. de Queiroz, M. Moore, S. D. Zhang, S. Y. Chen, J. Wang, D. Z. Li, T. S. Yi, unpubl. data). Hence, concatenation of the mitochondrial genes does in fact provide useful information for resolving deep phylogenetic relationships of Papilionoideae.

Fig. 5.  Phylogenetic relationships inferred using maximum likelihood (ML) analyses from 33 mitochondrial protein-coding genes. Bootstrap proportions and Bayesian posterior probabilities are 100% and 1.0 on all branches. The newly sequenced species is marked with a star.

Angiosperm mitochondrial genomes usually have low structural conservation (Knoop 2012; Mower et al. 2012b; Guo et al. 2016). Synteny analyses of Papilionoideae mitochondrial genomes showed extreme structural variations, including a large amount of re-arrangement and inversions (Fig. S3, available as Supplementary material to this paper). Blocks of ~136−261 kb in the mitochondrial genome of C. australe were homologous to those in the remaining nine genomes compared here. However, closely related taxa shared more syntenic blocks. For example, 48.3% (261 763 bp) of the 542 079-bp mitochondrial genome of C. australe was homologous to the more closely related Styphnolobium japonicum, but only 33.5% (181 693 bp) and 33.1% (179 660 bp) were syntenic to the more distantly related Ammopiptanthus mongolicus and Glycine max respectively (Fig. 6). This variation is not surprising, given that recombination events in plant mitochondria are known to play a major role in the structural variation of mitochondrial genomes (Arrieta-Montiel and Mackenzie 2011; Knoop et al. 2011; Kühn and Gualberto 2012), including within Papilionoideae mitochondrial genomes (Chang et al. 2013; Shi et al. 2018).

Fig. 6.  Syntenic blocks shared between mitochondrial genomes of Castanospermum australe and nine additional Papilionoideae species as shown by MAUVE alignments. In each comparison, arrows on the left indicate the length and percentage of the top genome that is homologous to the bottom genome, whereas those on the right indicate the length and percentage of the bottom genome that is homologous to the top genome. Species abbreviations are same as those of Table 1.

The comparison of the C. australe mitochondrial genome with nine other mitochondrial genomes of subfamily Papilionoideae has provided a much-improved understanding of mitochondrial genomic evolution across papilionoid legumes. Although these genomes share strong conservation of gene and intron content as well as some syntenic blocks and similarities in RNA editing, their variable genome size and structure indicate that recombination within Papilionoideae mitochondrial genomes has, nevertheless, been extensive. Equally important, our study demonstrated the value of mitochondrial genes for resolving deep-level phylogenetic relationships in legumes. Future research should continue to target legume mitochondrial genomes for sequencing and assembly, to further develop our understanding of the evolution of this overlooked genome, and to mine it for useful phylogenetic characters.

The authors declare that they have no conflicts of interest.

This work was supported by grants from the Strategic Priority Research Program of Chinese Academy of Sciences (XDB31010000), the National Natural Science Foundation of China [key international (regional) cooperative research project number 31720103903], and the Large-scale Scientific Facilities of the Chinese Academy of Sciences (number 2017-LSF-GBOWS-02).