Complete sequence of kenaf (Hibiscus cannabinus) mitochondrial genome and comparative analysis with the mitochondrial genomes of other plants

Plant mitochondrial (mt) genomes are species specific due to the vast of foreign DNA migration and frequent recombination of repeated sequences. Sequencing of the mt genome of kenaf (Hibiscus cannabinus) is essential for elucidating its evolutionary characteristics. In the present study, single-molecule real-time sequencing technology (SMRT) was used to sequence the complete mt genome of kenaf. Results showed that the complete kenaf mt genome was 569,915 bp long and consisted of 62 genes, including 36 protein-coding, 3 rRNA and 23 tRNA genes. Twenty-five introns were found among nine of the 36 protein-coding genes, and five introns were trans-spliced. A comparative analysis with other plant mt genomes showed that four syntenic gene clusters were conserved in all plant mtDNAs. Fifteen chloroplast-derived fragments were strongly associated with mt genes, including the intact sequences of the chloroplast genes psaA, ndhB and rps7. According to the plant mt genome evolution analysis, some ribosomal protein genes and succinate dehydrogenase genes were frequently lost during the evolution of angiosperms. Our data suggest that the kenaf mt genome retained evolutionarily conserved characteristics. Overall, the complete sequencing of the kenaf mt genome provides additional information and enhances our better understanding of mt genomic evolution across angiosperms.


Results
Kenaf mitochondrial genome sequencing and assembly. Isolated kenaf mitochondrial DNA (mtDNA) was used to construct a library for sequencing using PacBio RS II single-molecule real-time sequencing technology (SMRT), which generated 1.12 G of raw data, with an average read length of 4.6 kb, and the longest read was 32 kb. In total 67,152 reads (363,717,023 bp) were obtained after removing the adapter and low-quality regions, and the average coverage depth, read length and read quality were 605×, 5.4 kb, and 0.81, respectively (Supplementary Table S1). In total 1,819 reads (12,114, 267 bp, average length of 6.7 kb) were obtained after correcting by mapping the short reads to the long seed reads. After filtering the chloroplast reads, 1,762 reads (11,733,852 bp, average length of 6.7 kb) were used for the assembly process. Finally, the kenaf mt genome was assembled into a single circular molecule with a total length of 569,915 bp and an overall GC content of 44.9% (Fig. 1, Supplementary Table S2).
Gene content in kenaf. Sixty-two genes including 36 core protein-coding genes, conserved among all plant mt genomes were annotated by comparing the assembled kenaf mt sequence with known plant mt sequences in the NCBI public DNA database using BLASTn. The kenaf mt genome contains 20, 7, 4, 4, 3, 29, 1 and 1 genes responsible for electron transport, oxidative phosphorylation, small ribosomal proteins, large ribosomal proteins, cytochrome C maturation protein, rRNAs, tRNAs, and matR and mttB, respectively (Supplementary Table S3). Most protein-coding genes, except for sdh3, rps13 and rps19, were identical to those in the mt genome of the Gossypium species (Fig. 2). Twenty-three tRNA genes specifying 18 amino acids were identified in the kenaf mt genome. Of these genes, 15 tRNAs had a mt origin, and eight tRNA had a chloroplast origin (Table S3). The presence and locations of these genes in the kenaf mt genome and comparisons with other plant mt genomes are shown in Fig. 1 and Supplementary Table S4.
In most spermatophytes, the genes responsible for the electron transport chain and oxidative phosphorylation are conserved, except for mitochondrial complex II, which contains sdh3 and sdh4. Notably, the high diversity in the gene content among the higher plant mt genomes was a primary contributor to the variety of ribosomal protein genes (Fig. 2). The mt genome of plants is known to contain genes encoding products involved in electron transport, oxidative phosphorylation, ATP synthesis, cytochrome c biogenesis, ribosomes, and the translation of proteins (Fig. 2).
The protein-coding genes in the kenaf mt genome account for 6.9% of the genome and a total length of 39,534 bp. In addition, 126 open reading frames (ORFs) larger than 100 amino-acid residues in size were annotated in the kenaf mt genome (Supplementary Table S2). However, none of these ORFs could be assigned a function based on sequence similarity at either the nucleotide or protein level. Most of ORFs were considered hypothetical proteins. A putative protein of 295 amino-acid residues encoded by ORF295 has unknown functions, although its first 183 nucleotides were similar to those of rps4 from the 5′ to 3′ end ( Supplementary Fig. S2). The nucleotide sequence in the coding region had no similarity to any other plant mt genomes, except for a part of the Gossypium mt genome. This chimeric characteristic of ORF295 may have resulted from a horizontal gene transfer (HGT) event between angiosperm mt genomes. Interestingly, despite the large size differences among the mt genomes of various higher plant species, these genomes share a similar set of functional genes (protein, rRNA and tRNA genes), which is consistent with the results reported by Mower et al. 27 . However, the additional ORF e.g., ORFs identified in the kenaf mt genome, were not shared even among closely related plants, suggesting that many ORFs likely do not encode functional proteins and may have unidentified species-specific functions.
Repeat sequences of kenaf mitochondrial DNA. Repeat sequences are extensively found in the plant mt genome, are characterized primarily as forward repeat and palindromic repeats, and exhibit high levels of polymorphism. In the present study, we identified 584 repeat sequences that ranged from 20 to 7,782 bp and accounted for 11.71% of the total kenaf mt genome (Supplementary Table S2, Fig. S3). Most repeats (approximately 95%) were between 20 and 100 bp in length, accounting for 6.63% of the total genome; approximately 5% (28) of the repeats were larger than 100 bp, and three repeats were larger than 1 kb (R1, 7,782 bp; R2, 1, 877 bp; and R3 1,528 bp) ( Table 1). Most repeat sequences (≥60 bp) contained 2 copies of the repeat, and eight repeat sequences contained three copies (Table 1). Introns. In the kenaf mt genome, nine mt genes, composed of 25 introns ranging from 41 to 2,878 bp in size were identified, and occupied 7.8% of the total kenaf mt genome (Supplementary Table S4). Five of the nine mt genes were nad1, nad2, nad4, nad5 and nad7, which are components of mitochondrial complex I, and the remaining four genes were cox2, ccmFc, rps3 and rps10. In addition, five trans-spliced introns observed in nad1, nad2 and nad5 were fragmented into separate coding regions, which is consistent with angiosperm plants. Twenty cis-spliced intron sequences were observed in the remaining mt genes. The intron locations and splicing were highly similar to those observed in other higher plant mt genomes (Fig. 2).
Chloroplast-like sequences. BLASTn was used to identify chloroplast-like sequences in the kenaf mt genome. Twelve such sequence fragments were identified and showed >97% nucleotide sequence identity with the corresponding chloroplast sequences, and the segments ranged from 73 bp to 2,653 bp with a total length of 11,281 bp (accounting for 1.98% of the genome size). These chloroplast-derived fragments included eight tRNA-related sequences. Moreover, three intact chloroplast-related genes, i.e., ndhB, psaA and rps7, were identified in the kenaf mt genome (Supplementary Table S5).
Gene organization and gene clusters in plant mt genomes. The gene organization greatly differs among plant mt genomes. In this study, we compared the gene orders in the 28 mt genomes and counted the number of syntenic gene clusters (genes that remain in the same order). Four gene clusters (i.e., rrn5-rrn18, nad1-matR, rps12-nad3, and rps3-rpl16) were found to be highly conserved in the plant mt genomes (Fig. 3, Supplementary Table S6). The gene cluster cox3-sdh4 was widely distributed in most dicotyledonous species, except for Brassicaceae, while the conserved rpl5-rps14 gene cluster was scattered in the other dicotyledonous species, but present in all Brassica spp. Understandably, species that have close evolutionary relationships share more clusters. Each gene cluster is transcribed from the same strand, implying that the genes may undergo co-transcription as a polycistronic mRNA. Distribution of tRNAs and DNA transfer from the plastid to mitochondrial DNA. A complete set of tRNAs is essential for protein translation in the plant mt genome. However, many tRNAs undergo loss, migration and inactivation during mt genome evolution in higher plants 27 . To evaluate the origin and distribution of the tRNA genes, tRNA scan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) was used to predict the number and types of tRNA genes in the kenaf mt genome. In total, 23 tRNA genes were identified, and these genes recognized 18 amino acids (i.e., Asp, Gly, Met, Ser, His, Phe, Pro, Glu, Cys, Asn, Tyr, Trp, Asp, Lys, Ser, Leu, Ile, and Val). Thus, tRNA genes for two amino acids (i.e., Ala and Thr) were not identified and appeared to be missing from the kenaf mt genome (Fig. 4). Of these 23 tRNAs, eight had a plastid origin, and twenty-one had a mt origin.
The mt genomes of twenty-eight land plants and the fungal species G. lucidum were analysed to explore the patterns of tRNA loss during the evolution of plant mt genomes. Only G. lucidum has a complete set of tRNAs (Fig. 4). The trnA gene was lost from gymnosperms to angiosperms, indicating that trnA was lost early in the evolution of land plants. The trnG gene was absent from monocots, but existed in dicotyledons, suggesting that this gene was specifically present in dicotyledons. Although trnL, trnR, trnT and trnV were lost during the evolution of angiosperms, trnR and trnV existed in certain dicotyledons, suggesting that these genes may have been subsequently regained. Interestingly, most of the tRNAs in C. melo exhibited a pattern of plastid-like origin, suggesting that frequent exchanges occurred between the mt genome and the chloroplast genome.
BLASTn was used to assess the mt sequence fragments that originated in the chloroplast. Four chloroplast-derived fragments (trnH, trnM, trnN and trnW) were found to be conserved in all analysed mt genomes, and one (trnD) and two (trnC and trnF) chloroplast-derived fragments were found to be conserved in dicots and monocots, respectively. In contrast, other chloroplast-like tRNA genes exhibited scattered distributions, and certain native tRNA genes were irregularly lost among the higher plant mt genomes, suggesting that the gain and loss events of the tRNA genes occurred multiple times during evolution. Overall, trnC, trnE, trnK, trnM,   Table 1. Repeats (≥60 bp) in the kenaf mt genome. a Compare with copy-1 as control. b The letters F and P represent forward repeats and palindromic repeats, respectively. The numbers listed in the starting and ending points refer to positions in the kenaf mt genome sequence (GenBank accession MF163174). and gymnosperms, and bryophytes was chosen as the outgroup. The chloroplast-derived sequences and non-protein-coding sequences were removed before blasting against the other mt genomes. First, these mt functional genes were concatenated in a head-to-tail format. Maximum likelihood method was used to complete the phylogenetic tree analysis. As shown in Fig. 5, the Hibiscus cannabinus and Gossypium species belonging to the Malvaceae family were classified into one clade with a high bootstrap support value of 100. In addition, the species share a high sequence similarity, as supported by the higher bootstrap support values. Species belonging to different groups share less sequence similarity and have reduced bootstrap support value. The phylogenetic tree strongly supported the separation of monocot plants and dicot plants, and the separation of angiosperms from gymnosperms. Additionally, the evolutionary relationship of these 28 plant species was analysed using the plant taxonomy method and used to construct an NCBI taxonomy common tree (Fig. 6). The phylogenetic relationships based on mt genome homologous sequences are consistent with the species taxonomy despite the exceptional variability among these mt genomes.
To further explore the utility of these mt genes in phylogenetic reconstruction, twenty-two mt genes were divided into five groups according to the function of their proteins (Supplementary Table S3), and the genes in each group were assembled in a head-to-tail arrangement. Among the five groups of phylogenetic trees, the set of mitochondrial complex I genes was congruent with a previous reconstruction based on 22 functionally related genes (Supplementary Fig. S4). The trees of mitochondrial complex III and complex IV reconstructing the divergence of monocots and dicots resulted in topologies that differed from those obtained by the previous reconstruction base on 22 functionally related genes, but the species fell into clades that belonged to the same family (Supplementary Figs S5, S6). The phylogenetic trees of the mitochondrial complex V and cytochrome c biogenesis genes revealed evolution relationships that slightly differed from those obtained with the previous reconstruction base on 22 functional related genes (Supplementary Figs S7, S8). Therefore, the phylogenetic analysis base on the function of mitochondrial genes revealed slightly different topologies, but the species fell into clades that were consistent with their family designations. In contrast, the phylogenetic tree based on the whole set of genes was congruent with the species taxonomic tree (Fig. 6).

Discussion
Third-generation SMRT sequencing technology based on the PacBio RS II platform can produce substantially longer reads (>5 k/read) than second generation sequencing 28 . Which also can be used to closed genome gaps, whole-genome sequencing projects for any species [29][30][31][32][33][34][35] and non-SNP DNA variations identification 36 . In our study, the SMRT sequencing technology was used to sequence the whole mt genome sequence of kenaf (Hibiscus  Table S2 and include 17 respiratory complex genes, four cytochrome c biogenesis genes and a cob gene, and the tree was rooted with Marchantia paleacea.  . We obtained the kenaf mt sequence with a high accuracy, and the genome size was 569,915 bp. The longest read was 32 kb, which is much longer than the usual reads obtained using other sequencing technologies.

Characteristics of plant mitochondrial genes. Plant mtDNAs primarily comprise of protein-coding
genes, tRNAs and rRNAs. The kenaf mt genes included only 36 of the 41 protein-coding genes present in ancestral land plant mt genomes 27 , indicating that several protein-coding genes were lost or transferred to other organelles during the evolution of kenaf mitochondria. The frequent loss and functional transfer of ribosomal protein genes and succinate dehydrogenase (sdh) genes to the nuclear genome were the main causes of the variable gene contents among the plant mt genomes. This finding is consistent with previous results that have been confirmed by Southern blot hybridization 21 . Our results revealed the presence of only sdh4 in the kenaf mt genome, while both sdh3 and sdh4 were identified in the closely related Gossypium species 37 . Thus, the presence of the succinate dehydrogenase genes is highly variable, even among evolutionary close angiosperm species. However, two atp9 copies were identified in the kenaf mt genome, which may have resulted from HGT events or mtDNA recombination during the evolution of the kenaf mt genome.
Repeat sequences in the genome. The mt genomes of land plants, particularly angiosperms, are frequently characterized by repeat sequences 38 , which could explain most of the variation in the mt genome size. Moreover, these sequences are sites of intragenomic recombination, underlining the evolutionary changes in the mt genome organization in vivo 39,40 . Tandem simple and scattered repeat sequences are extensively found across plant mt genomes and exhibit high levels of polymorphisms 41,42 . The repetitive sequences in the Cucumis melo mt genome had a size of 2,738 kb and comprised primarily small repeats, accounting for 42.7% of the mt DNA 13 . In contrast, other genomes contain fewer larger segmental duplications 41,43 . The Vitis mt genome (773 kb) has only 6.8% repetitive DNA sequences 11 , while the moderately sized Legume vigna genome (401 kb) has fewer and smaller repeats that account for 2.7% of the mt genome. These data suggest that the genome size is not a good indicator of repeat content in angiosperm mt genomes. In the present study, the repetitive structure of the kenaf mt genome accounted for 11.71% of the genome. These sequences are poorly conserved across species and have a high proportion of smaller repeats, indicating that the increased size of the kenaf mt genome was primarily due to the duplication of short sequences. 44 to 34 in the hornwort Phaeoceros laevis 45 . According to the present study, the kenaf mt genome retained 25 introns, disrupting 9 protein genes, which is consistent with the common ancestor theory of gymnosperms and angiosperms 27 and suggests that the introns were lost or gained during plant evolution. However, cis-splicing is ubiquitous in most introns of seed plants, while the mt genes nad1, nad2 and nad5 evolved a split structure that requires trans-splicing and were highly consistent with the sequenced angiosperm mt genomes 27 . Thus, the transcription process of introns was conserved among the angiosperm mt genomes.

Introns in the mitochondrial genome. Land plant mt genomes contain a large and variable number of introns, ranging from 19 in Silene latifolia
Conservation of gene clusters. Genome recombination can disrupt clusters, while multiple recombination events can generate similar syntenic gene clusters, leading to vast differences in the gene order among plant mt genomes 46 . In general, evolutionarily close species have more similar gene orders and clusters. Two gene clusters (rrn18-rrn5 and nad5-nad1-matR) are conserved in all land plant mt genomes, and may date back to the original plant mt genomes of liverworts, mosses, and most charophytes 27 . In contrast, the rps12-nad3 and (rps19)-rps3-rpl16 gene clusters were evolutionary conserved in most land plant mt genomes but absent from M. paleacea and S. suchowensis, respectively 47,48 . Two other gene clusters (rps10-cox1 and sdh4-cox3-(atp8)) were specifically conserved in dicots, except for Brassica 49 . The cluster of atp4-nad4L exists in all surveyed dicots, except for the species of Gossypium, H.cannabinus and C. capsularis. These exceptions were likely due to frequent recombination events during plant mt genome evolution. In our study, the characteristics of the gene clusters in the kenaf mt genome were consistent with the general conservation of most dicotyledon, indicating that the gene clusters in the kenaf mt genome were more conserved during plant mt genome evolution.

DNA transfer in the mitochondrial genome.
HGT is thought to be the main process during the acquisition of exogenous sequences [50][51][52] . The transfer of DNA sequences among plastid, nuclear and mt genomes is a common phenomenon that has been observed in the fully sequenced mt genomes of land plants 53,54 . Although the amount of plastid DNA in the mt genome is 3-6% in most examined species, plastid DNA varies from 2 kb (0.5%) in the Legume Vigna to 113 kb (11.5%) in Cucurbita pepo 12 . In many cases, these plastid-to-mitochondrion transfers have resulted in the insertion of plastid genes into the mt genome, but most genes are clearly nonfunctional 53 . Occasionally, the only plastid genes that are transferred into the mt genome and remain functional encode tRNAs. The absence of chloroplast-derived tRNAs from liverworts, moss, hornworts, and bryophytes indicates that DNA transfer from the chloroplast genome to the mt genome might have occurred after the divergence of gymnosperms and angiosperms. Eight chloroplast-derived tRNAs identified in the kenaf mt genome can be traced to the retention of an earlier HGT event.
The plastid-derived trnW (GTT) and trnH (GTG) genes are frequently observed in angiosperms but are absent from C. taitungensis, indicating that these tRNA genes may have been transferred after the separation of angiosperms. Identifying the numbers and types of tRNA genes in the kenaf mt genome may be helpful for evaluating the origin and evolution of tRNA genes in higher plants. These results suggest that the intracellular transfer of tRNA and ribosomal genes from the chloroplast to the mitochondria was a frequent process.

Conclusion
Plant mt genomes are intriguing due to their highly conserved genic content and slow rate of genic evolution. In contrast, features, such as the genomic structure, the genome size and repeat sequences, are highly variable. In this study, we determined the complete sequence of the kenaf mt genome. The comparison of the kenaf mt genomic features with those of other plant mt genomes should provided a more comprehensive understanding of mt genome evolution in higher plants. The complete mt genome of kenaf shares many common genomic characteristics with other plant mt genomes, such as the conservation of genic content, gene clusters, certain intergenic sequences and tRNA gene origin and distribution. These observations suggest that the evolution of mt genomes is consistent with the species relationships in plant taxonomy. However, the highly dynamic genome structures (genome size and gene order) suggest that the recombination of higher plant mt genomes is independent and random among species. The sequencing of the kenaf mt genome contributes to our understanding of the characteristics of the mt genome across angiosperm evolution.

Materials and Methods
Mitochondrial DNA isolation and sequencing. Mitochondria were isolated from the kenaf maintainer line UG93B and purified from 7-day-old etiolated seedlings using differential centrifugation and discontinuous (18%, 23% and 40%) Percoll density gradient centrifugation according to the methods described by Wilson and Chourey 55 . The mitochondrial DNA isolation was performed as described by Sue 56 with modifications. The purified mitochondria were lysed with cetyltrimethylammonium bromide (CTAB) supplemented with 2% polyvinylpyrrolidone and 0.7% β-mercaptoethanol (Solarbio, Beijing) at 65 °C for 30 min. The lysis solution was extracted two to three times with chloroform/isoamyl alcohol (24:1), and absolute ethyl alcohol was used to precipitate the mtDNA. DNase-free water (50 μL) was added to resuspend the DNA pellets. The integrity, quality and concentration of the UG93B mtDNA were analysed using agarose gel electrophoresis, a NanoDrop 2000 (Thermo Scientific, USA) and a Qubit fluorometer (Thermo Scientific, USA).
DNA sequencing and genome assembly. In total, 20 μg of UG93B mtDNA were randomly sheared to fragments using a Covaris S220 (Thermo Scientific, USA). Large fragments with an average size of 20 kb were purified by magnetic bead enrichment. SMRTbell templates were obtained by ligating the hairpin adaptors to the end of a double-stranded DNA molecule and removing the failed ligation products with exonuclease. An Agilent 2100 Bioanalyzer High Sensitivity Kit (Agilent Technologies, USA) was used to assess the quality of the library. Subsequently, eight SMRT cells were sequenced using P4-C2 reagents on a PacBio RS II sequencing platform 57 . The sequencing and de novo assembly were performed at Nextomics Biosciences Co., Ltd, Wuhan, China. The clean reads were obtained by filtering out the sequencing adapters and low-quality sequences using SMRT Analysis 2.3.0 with the default settings. The kenaf mt genome sequence was extracted from filtered reads containing both chloroplast and mt genomes. Blat 58 was used with the default parameters against the NCBI chloroplast genome data to filter reads containing chloroplast genomes, and reads with a match greater more than 90% were moved. The kenaf mt genome sequence was assembled using a Hierarchical Genome Assembly Process (HGAP) workflow, including preassembly, error correction, Celera assembly and polishing using Quiver 59 . The long reads were selected as "seed" reads, to recruit all other subreads and construct highly accurate preassembled reads using a directed acyclic graph-based consensus procedure. This procedure was followed by assembly using off-the-shelf long-read assemblers. A basic local alignment with successive refinement (BLASR) was used to align the short reads to the seed reads and improve the sequence accuracy 60 . The Celera assembler software and overlap-layout-consensus (OLC) algorithm were used to assemble all corrected contigs 61 . Finally, Quiver was used to improve the site-specific consensus and generate the gap-free kenaf mt genome. Additionally, a specific prime pair was designed to verify the circular mitochondrial genome of kenaf ( Supplementary Fig. S1).
Genome annotations and analyses. The mt genomes were annotated using BLASTn and MITOFY 12 using previous angiosperm mt genes to query sequences in the NCBI database (https://www.blast.ncbi. nlm.nih.gov). The tRNA genes were identified using the tRNA scan-SE software (http://lowelab.ucsc.edu/ tRNAscan-SE/). ORFs that contained more than 100 amino-acid residues and started with methionine were predicted and annotated using ORF-Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). The repeat sequences were analysed using REPuter software (http://bibiserv.techfak.uni-bielefeld.de/reputer) with the following parameters: the repeat sequence was at least 20 bp in length and the repeat identity was greater than 90% 62 . The circular map and syntenic gene cluster maps of the plant mt genomes were created using OGDRAW v1.2 (http://ogdraw.mpimp-golm.mpg.de/) 63  nih.gov/genome/organelle/). These mt genome sequences were selected because they are available for analysis in NCBI and are clearly taxonomically classified. Phylogenetic analyses were performed using concatenated exon sequences from 22 conserved protein-coding genes (atp1, atp4, atp6, atp8, atp9, ccmB, ccmC, ccmFc, ccmFn, cob,  cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7 and nad9) extracted from these 28 plant mt genomes. These nucleotides were aligned using ClustalW and manually modified to eliminate gaps and missing data. Finally, the maximum likelihood (ML) method was used to construct original phylogenetic trees by MEGA 6.0 64 . The bootstrap replications were performed with 1000 according to Felsenstein 65 . The evolutionary distances were computed using the Kimura 2-parameter method 66 and the tree was rooted with Marchantia paleacea. The NCBI taxonomy common tree was described by Federhen 67 and constructed using the online NCBI taxonomy database (https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/wwwcmt.cgi).