The complete mitochondrial genome of the early flowering plant Nymphaea colorata is highly repetitive with low recombination

General features of Nymphaea mitogenome

The Nymphaea mitogenome is assembled into a single circular molecule of 617,195 bp (Fig. 1), a size larger than ca. 80% of the currently sequenced vascular plant mitogenomes (as of July 2018). The relatively large size of Nymphaea mtDNA is primarily due to its abundant repetitive sequences, which add up to 301,676 bp and account for nearly half (49%) of the mitogenome, in contrast to most of other vascular plant mitogenomes with repeat ratio generally below 30% (Additional file 1: Table S1). The Nymphaea mitogenome encodes 41 protein genes, three rRNA genes (rrn5, rrn18 and rrn26), and 20 tRNA genes (13 mitochondrial native and seven plastid derived) (Table 1). Intergenic spacers constitute the largest part (519,361 bp, 84%) of the Nymphaea mtDNA, and protein coding sequences comprise only 6% (35,961 bp) of the total length. In general, the gene content of Nymphaea is very similar to the other published angiosperm mitogenomes, especially to Amborella [48] and Liriodendron [6]. Nymphaea mt gene set differs from Amborella only by its presence of the functional protein-coding gene rps10 that is pseudogenized in Amborella, whereas differs from Liriodendron by its presence of plastid derived tRNA gene trnL(CAA)-pt and absence of trnV(TAC). Repeat-induced duplicated genes are widespread in vascular plants [49], such as Nelumbo nucifera possesses six duplicated protein genes [50] and maize (CMS-C) contains 10 duplicated protein genes [51]. In Nymphaea mitogenome, rps19 and atp6 each presents as two copies. The duplicated rps19 are identical, while the two copies of atp6 are different in length, with one copy 114 bp longer at the 3′ terminal. The shorter version of Nymphaea atp6 is still longer than that of Amborella and Liriodendron by 36 bp and 75 bp at the 5′ terminal. Blastn and Blastp searches of the 114-bp nucleotide sequence and the translated amino acid sequence against the NCBI database do not return any hits, suggesting a probably chimeric origin of atp6_D2 (the longer copy), via gene fusion of atp6 with Nymphaea specific intergenic spacer sequence at some evolutionary stage. Considering the majority of the two atp6 copies located in a pair of identical inverted repeats of 3293 bp at a distance of 196 Kb, the identical 882 bp of the two atp6 copies may be indicative of the result of repeat recombination in homogenization of the gene copies carried [7]. We further checked all intergenic spacers for possible pseudogene pieces using 68 annotated Nymphaea coding regions as queries. Altogether, we identified 52 pseudogene fragments ranging from 28 bp to 182 bp, which matched nine protein coding genes (nad5.× 2.× 5, rpl2.× 1, co× 1, atp6, ccmC, nad6, atp8, rrn18, rrn26) with identities ranging from 85 to 100%. Two largest pseudogene pieces of atp6 (182 bp) and rpl2.× 1 (142 bp) formed Nymphaea specific chimeric ORFs with parts of the adjacent intergenic spacer sequences, which, in some cases, may cause cytoplasmic male sterility (CMS) [52]. Blastn search of all these pseudogene fragments against the NCBI nucleotide database yielded much lower similarities with any other species than Nymphaea, indicating the origin of these gene vestiges from intragenomic recombination events [25] rather than horizontal gene transfers from other plants.

Genome feature	Amborella trichopoda	Nymphaea colorata	Liriodendron tulipifera
Accession	KF754799–KF754803	KY889142	KC821969
Size (bp)	3,866,039	617,195	553,721
GC%	45.9%	45.1%	47.7%
Genes	63	64	64
tRNAs	20	20	20
rRNAs	3	3	3
Protein coding	40	41	41
Cis-spliced introns	19	19	20
Trans-spliced introns	6	6	5
Gross length of repeats (> 19 bp)	914,403 (23.65%)	301,777 (48.89%)	85,993 (15.53%)
Plastid-derived (bp)	138,313 (3.58%)	13,170 (2.13%)	24,807 (4.48%)
Total gene length (bp)	85,152 (2.20%)	97,834 (15.85%)	82,284 (14.86%)
Protein exons (bp)	33,866 (0.88%)	35,961 (5.83%)	35,324 (6.38%)
Cis-spliced intron length (bp)	38,896 (1.01%)	54,840 (8.89%)	33,631 (6.07%)

Nymphaea shares 27% (168,686 bp) of its mtDNA with other sequenced plant mitogenomes with nearly half occurred in the genic region, and the other half (95,941 bp) in the intergenic region, accounting for 15% of the mitogenome. Nymphaea shares its intergenic spacer sequences the most with Amborella (58,049 bp), and Liriodendron (28,715 bp), then Phoenix (26,790 bp). As multiple lines of evidence suggested a divergence time of Nymphaea from the rest of angiosperms at 180 Mya (www.timetree.org), the seemingly low level of sequence sharing between Nymphaea and other angiosperms fits well to the regression line generated by analyzing 14 phylogenetically independent seed plant taxa [31], suggesting the generally high divergence nature of angiosperm mitogenomes. For example, Citrullus lanatus [14] shares with Vitis vinifera [49] 72,313 bp of its intergenic spacers despite a divergence time of 105–115 Mya; Carica papaya [53] shares with Nicotiana tabacum [54] 66,327 bp with a divergence time of 110–124 Mya.

The Nymphaea mitogenome contains 25 group II introns, including 19 cis-spliced and six trans-spliced introns (nad1i394g2, nad1i669g2, nad1i728g2, nad2i542g2, nad5i1455g2, nad5i1477g2), which is similar to the intron set of Amborella [48] and Phoenix [55], but differs from Liriodendron by its presence of the trans-splicing nad1i728g2, which is a cis-spliced intron in Liriodendron. It is noteworthy that cox2i373g2 of Nymphaea reaches a length of 11.4 Kb, making it the longest organellar intron reported in plants to date. We checked the coverage of the genome assembly on this intron region. A continual and even coverage of cox2i373 and its cox2 exon regions indicated that the presence of this intron is unlikely yielded from an artifactual assembly result (Additional file 2: Figure S3). We mapped the transcriptomic reads to the mitogenome, but due to the low coverage of the transcriptome data we cannot figure out whether this intron is continually transcribed (Additional file 2: Figure S3). We aligned the Nymphaea cox2i373g2 with that of Triticum timopheevii [56], five out of the six conserved domains of this group II intron were well aligned, except for the domain IV, indicating this domain may expanded in Nymphaea (Additional file 3: Figure S4). Although we recognized Nymphaea cox2i373 as a cis-spliced intron here, we still cannot rule out the possibility that this intron is trans-spliced, but the two parts of the trans-spliced intron happens to locate proximately in the genome and in an orientation consistent with cis-splicing. Besides, intron rpl2i846g2 and nad4i976g2 exceed 6 Kb; intron nad2i1282g2, nad2i156g2, and nad7i917g2 exceed 3 Kb in length. Overall, the total length of the 19 cis-spliced introns add up to 55 Kb, comprising 9% of the whole mitochondrial genome, which is substantially higher than any other angiosperm mitogenomes sequenced to date in both absolute and percentage terms [13]. The highly repetitive nature of the Nymphaea mtDNA accounts for a large portion of its intron size expansion (Fig. 2). About 40% to 80% of the six large introns of the Nymphaea mitogenome (> 3 Kb) are made of repetitive elements, a phenomenon similar to what observed in ferns [9] (Additional file 4: Table S2).

Repeats and homologous recombinations

Blastn searches identified 1,188,860 repeated sequences that are longer than 30 bp and with unique begin-end coordinates in an overlapping fashion, accounting for nearly half (49%, 301,676 bp) of the Nymphaea mitogenome (Additional file 1: Table S1). These numerous imperfect and partially overlapping repeated sequences in Nymphaea constitute 886,983 repeat pairs, with the length distribution mainly in the range of 100–200 bp and identity distribution mainly between 80 and 95% (Additional file 5: Figure S1). Cd-hit-est as implemented in the cdhit suite [57] recovered 290 families with an identity threshold of 0.8 and a word size of five out of the total repeated sequences using a greedy incremental clustering algorithm method. The representatives of these repeat families were subsequently checked for occurrences using blastn searches against NCBI nucleotide database. Most (252, 87%) of these repeat families are restricted to Nymphaea and are unique in plants, only 38 are shared with other plant mitogenomes, such as 22 with Amborella, 18 with Liriodendron, 11 with Arabidopsis, 11 with Gymnosperms, four with ferns, eight with bryophytes, and nine with charophycean green algae. The observed low repeated sequence sharing of Nymphaea with other plant mitogenomes reflected a commonplace phenomenon of wild divergence of intergenic spacers as has been exemplified by the remarkable intraspecific variation in four mitogenomes of Silene vulgaris [25].

Benefited from the deep sequencing of PacBio long reads (601×, average 7294 bp, Additional file 6: Figure S2), we were able to detect minor recombinations at a frequency as low as 1/1200. A total of 886,983 repeat pairs with length ranging from 30 to 3293 bp and blast identity above 80% were examined for recombination activity (Additional file 5: Figure S1). Unexpectedly, only ten repeat pairs show evidence of recombinations with one to 48 recombined reads detected for each repeat pair (Table 2). Three direct repeats and seven inverted repeats recombined at frequencies ranging from 0.07 to 8.18%, which could possibly give rise to a set of alternative mtDNA configurations and subgenomes via inversions and subdivisions of the master conformations (Fig. 3). According to our observations, a majority of the repeats (R3–R10) recombined at a frequency below 1%; only two repeats yielded more than 10 recombined reads, including the longest inverted repeats of 3293 bp with 48 recombined reads and a self-inverted repeats of 128 bp with 13 recombined reads detected, suggesting alternative conformations (ACs) with a full set of genes rearranged are more abundant than subgenomes with reduced gene set in Nymphaea mitochondria, which resembles that observed in fern Ophioglossum with predominant ACs harboring inversions induced by the longest 4-Kb inverted repeats recombining at a frequency of 24.5% and a small number of subgenomes generated by recombinationally less active medium-sized repeats recombining at a frequency less than 2.5% [9].

Repeat	Length (bp)	Identity (%)	Direction	Position	Reads support master circle conformation	Reads support alternative conformation
R1	3293	100.00	+	153,171–156,463	539 (91.82%)	48 (8.18%)
			–	355,403–352,111
R2	128	82.81	+	414,017–414,144	1000 (98.72%)	13 (1.28%)
			–	414,144–414,017
R3	462	100.00	+	187,030–187,491	874 (99.66%)	3 (0.34%)
			+	307,094–307,555
R4	1538	97.33	+	386,370–387,906	822 (99.76%)	2 (0.24%)
			–	389,671–388,149
R5	508	81.50	+	15,484–15,960	960 (99.79%)	2 (0.21%)
			–	388,795–388,149
R6	399	99.75	+	39,755–40,153	880 (99.77%)	2 (0.23%)
			+	466,916–467,314
R7	945	98.41	+	301,710–302,649	915 (99.89%)	1 (0.11%)
			+	357,890–358,834
R8	729	99.86	+	313,110–313,838	932 (99.89%)	1 (0.11%)
			–	315,533–314,805
R9	692	99.86	+	80,337–81,028	783 (99.87%)	1 (0.13%)
			–	309,412–308,722
R10	72	97.22	+	302,169–302,240	1362 (99.93%)	1 (0.07%)
			+	313,560–313,631

Recombination involving large repeats generally result in equimolar or nearly equimolar recombined molecules in the genome [58], as exemplified in Silene latifolia, Silene vulgaris [7], Mimulus guttatus [33], and Ginkgo [31]. In our study, two large repeats (out of 2224 repeats) with a length of 3293 bp and 1538 bp show evidence of recombination but with low recombination frequency at only 0.24% and 8.18%. Such low recombination frequency has also been seen in other plant mitogenomes, for example, in Silene conica and Silene noctiflora [7], tens of large repeats induced recombinations at a frequency around 5%; in Ginkgo biloba, a large repeat of 1.5 Kb recombined at a frequency of 9%; in fern Ophioglossum californicum, two large repeats with 4-Kb and 1-Kb induced recombinations at frequency of 24.5% and 0.1%, respectively [9]. In addition to these large repeats, all seven medium sized repeats of the Nymphaea mitogenome recombined even more rarely with recombination frequencies ranging from 0.11 to 0.34%, which is similar to the observation in the gymnosperms Ginkgo biloba and Welwitschia mirabilis [31], ferns [9], the flowering plants Cucumis sativus [13] and Vigna angularis [26], whereas significantly lower than those observed from Viscum scurruloideum [10], Silene latifolia and Silene vulgaris. For example, in Viscum, three medium-sized repeats result in recombination equilibrium, and another two recombined actively at frequencies of 11.2% and 38.7%.

The low recombination level of the Nymphaea mitogenome is further evidenced by a small number of recombined reads detected, i.e., only 74 reads were found to support alternative configurations resulted from repeat-mediated recombinations (Additional file 7: FASTA.fa), accounting for only 0.13% of the total reads (74 out of 56,849), which is apparently lower than those found in other plants, such as 10% in Mimulus guttatus [33], 6.6% in Ophioglossum, and 2.2% in Psilotum. The low level of recombination rate found in Nymphaea suggests the predominant existence of master configuration in vivo in this plant, with a low level of substoichiometric recombinant forms. The latter has been proved to exert profound effect on plant growth, such as cytoplasmic male sterility (CMS) [59–61] and abnormal growth phenotypes [62, 63].

Understanding the paradoxical coexistence of the low recombination and abundant repeats in mitogenomes, such as Nymphaea and Ophioglossum, must take into account the nucleus’ control over the accuracy of the repair of mitochondrial chromosomes by a series of nuclear-encoded and mitochondrial targeted factors [64, 65]. Disruption of these genes could initiate and promote mitochondrial intragenomic recombination [58], as have been documented in Physcomitrella [66] and Arabidopsis [24]. Such nuclear genes may be under different levels of selection pressure, resulting in distinctive stability of mitogenomes in specific plant groups. For example, in each of the major bryophyte lineage, mitochondrial genomes kept a high degree of structural conservation over long period of evolution [2], which is in contrast to the observations in Silene vulgaris [25] and Beta vulgaris [67] with remarkable intraspecific mitogenome rearrangements.

Plastid DNA insertions

Conserved gene clusters

Plant mitogenomes are highly fluid in genome structure due to the repeat mediated homologous recombinations, sequence duplications, genome expansion and shrinkage, and incorporation of foreign DNAs [71], whereas some gene clusters are conserved across large phylogenetic scale [6, 50, 72]. The relatively low recombination level observed in Nymphaea does not necessarily predict strictly conserved genome arrangement compared with the ‘fossilized’ angiosperm mitogenome of Liriodendron or the other early angiosperm Amborella, as we found 38 and 44 rearrangements between mitogenomes of Nymphaea-Amborella and Nymphaea-Liriodendron, respectively. Nevertheless, in comparison of gene order of Nymphaea with that of the 214 plant mitogenomes, we identified 11 conserved gene clusters in Nymphaea, of which, three (rpl2–rps19–rps3–rpl16, rps13–rps11, and rrn18–rrn5) could be dated back deeply to the origin of mitochondrion from its endosymbiont bacterial ancestor [72]. The cluster trnfM(CAU)–rrn26 is widely distributed in streptophytes. Four clusters (cox3–sdh4, nad3–rps12, rpl5–rps14–cob, and rps10–cox1) emerged since gymnosperms. The cluster trnP(UGG)–sdh3 shows a sporadic distribution pattern in bryophytes, Ginkgo, Cycas and many angiosperms, indicative of the secondary loss of the gene cluster in lycophytes and ferns. The angiosperm conserved cluster trnP(UGG)-pt–trnW(CAA)-pt does not show up in Amborella, suggesting its emergence in Nymphaea or even earlier in the ancestor of angiosperms then secondary loss of the cluster in Amborella. The gene cluster <nad5.× 4.× 5 > <trnE (TTC)–nad7 > is only shared by three angiosperm species, namely, Liriodendron, Nelumbo and Nymphaea, suggesting its emergence in the ancestor of angiosperms followed by fast degeneration as a consequence of extensive genome rearrangements. However, the sporadic distribution of the cluster could more likely indicate a coincidence of independent structural evolutions in the three lineages (Additional file 8: Table S3).

Mitochondrial genome assembly and annotation

The mitochondrial genome of Nymphaea colorata was obtained from the genome project of Nymphaea colorata led by Liangsheng Zhang (unpublished data). The genome sequencing was performed on a PacBio RSII platform (Pacific Biosciences, Menlo Park, CA). The Raw PacBio reads were corrected to accuracy above 99% using the RS_PreAssembler, and then assembled into contigs using the program Canu (github.com/marbl/canu). Two mitochondrial contigs of 527,532 bp and 157,672 bp were identified using the NCBI Blast program with the Liriodendron mitochondrial genome as a reference. The two contigs overlapped with each other at both ends by 34,745 bp and 16,132 bp, and finally formed a circular molecular of 617,195 bp, with an average depth of 601×. RNA-seq data of Nymphaea colorata were also obtained from the genome project of Nymphaea colorata (unpublished data).

The annotation for the Nymphaea mitogenome was performed as previously described [3, 70]. Protein coding genes and rRNA genes were annotated by blastn searches of the non-redundant database at National Center for Biotechnology Information (NCBI). The exact gene and exon/intron boundaries were further confirmed in Geneious software (v.10.0.2, Biomatters, www.geneious.com) by aligning each gene to its orthologs from available annotated plant mitochondrial genomes at the NCBI website (www.ncbi.nlm.nih.gov/genome/organelle). The tRNA genes were detected using tRNAscan-SE 2.0 [73].

Repeats and repeat-mediated homologous recombinations

Repeats identification of Nymphaea and other vascular plant mitogenomes were carried out using NCBI blastn searches by searching the Nymphaea mitogenome sequence against itself with an e-value cut-off of 1e^− 6, and a word size of 7 following Guo et al. [9]. All the repeat sequences were subsequently extracted and clustered into difference families using the program cd-hit-est as implemented in cdhit suite v4.6.7 [57], with a word size of 5 and sequence similarity threshold of 0.8. We estimated the number of repeats from the number of unique begin-end coordinates of hits from blastn search according to Alverson et al. [13]. To detect the active repeat-mediated intragenomic recombinations within the PacBio reads, we built up an mt read database using corrected genome sequencing PacBio reads. We used the Nymphaea mitogenome sequence as the reference to blast the total Nymphaea PacBio reads database with an e-value cut-off of 1e^− 100 for extraction of mt reads, the resultant mt reads was further searched against Nymphaea plastid sequence with the same parameters to remove putative plastid reads with overall alignment coverage > 85% of the read length. Finally, we got a mitochondrial read database of 75,863 reads with an average length of 7294 bp, and total length 553,358,387 bp.

Repeat-mediated homologous recombinations were evaluated for those repeat pairs ranging from 30 to 3293 bp with blast identity > 80% following Alverson et al. [13]. Specifically, for each repeat pair, we built four or six reference sequences, each with 200 bp up- and down-stream of the two template sequences (original sequences), and the two (for repeat pair with identity =100) or four (for repeat pair with identity < 100) recombined sequences (alternative configurations) constructed from the putative recombination products. Then, we searched the reference sequences against the Nymphaea mt reads database, and count the number of matching reads with a blast identity above 99.5%, and a hit coverage over 200 bp in both flanking regions of each repeat sequence. After that, the templates with evidence of recombination were extracted and elongated in both sides to 2000 bp and searched again to the Nymphaea mt reads database to remove the recombinants with undersized flanking regions. Finally, the best matched reads for all the recombinants were extracted and aligned with the Nymphaea mitogenome in Geneious v10.0.2 (https://www.geneious.com/) to authenticate the accuracy of the recombined reads.

Identification of plastid derived sequences

To identify plastid derived mitochondrial sequences, the Nymphaea mitogenome was searched against the plastid genome of Nymphaea colorata (data unpublished), and all plant mitogenome database with an e-value cut-off of 1e^− 6 and a word size of 7, simultaneously. The blastn output was then visualized in Geneious v10.0.2 (https://www.geneious.com/) and each of the identified plastid sequence insert was compared with its co-occurring mt homologs from all other plant mitogenomes to infer the putative origin of the intracellular transfer.

Identification of conserved gene clusters

The gene orders of Nymphaea, Amborella, and Liriodendron were compared with each other using UniMoG [74] to identify rearrangements among three mitogenomes. The conserved gene clusters were identified if they appeared in any two of the three early angiosperms and simultaneously presented in at least one major plant group, e.g., lycophytes, ferns, gymnosperms, or angiosperms.