Broussonetia belongs to the family Moraceae, the nitrogen-fixing clade in the APG IV and contains the non-hybrid species B. monoica, B. kaempferi, B. papyrifera, and one hybrid species (Broussonetia × kazinoki). These species can be found across East Asia (China, Korea, Japan), Southeast Asia (Vietnam, Laos, Thailand, Cambodia, Myanmar), India, island Southeast Asia, Melanesia, and the Polynesia islands (Chang et al., 2015; Chung et al., 2017). Broussonetia papyrifera was essential in the development of the papermaking technology in ancient China and barkcloth-making in the Pacific Islands and Central America (Barker, 2002; Chang et al., 2015; Peng et al., 2019). Due to its important economical, medicinal, ecological attributes, B. papyrifera is being developed into a novel model system for woody plant research (Peng and Shen, 2018; Peng et al., 2019). The species Broussonetia monoica is a shrub (2–4 m tall) which exhibits notable leaf morphological diversity (Chung et al., 2017; Kuo et al., 2022), and Broussonetia kaempferi is a scandent shrub with highly suitable bark fibers that are used for making paper (Zhou and Michael, 2003). Recently, several studies have characterized the genus Broussonetia by molecular phylogeny (Chung et al., 2017), provided complete chloroplast genome assembly (Kuo et al., 2022) and an in-depth investigation into the chromosomal-scale genome of B. papyrifera (Peng et al., 2019). However, the dynamic evolution of mitochondrial gene and intron content of Broussonetia has never been assessed until now.

It has been greatly accepted that plant organellar genomes originate from ancient endosymbiotic bacteria approximately a billion years ago and play vital roles in massive essential life activities including photosynthesis, cellular respiration, and ATP synthesis (Woloszynska, 2009; Wu et al., 2022). The plastid genome (plastome)is a circular molecule, highly conserved in structure and gene content across most lineages, making it ideal for phylogeny (Wu and Ge, 2011; Wang et al., 2022). The mitochondrial genome (mitogenome) is semiautonomous encoding some self-related proteins but is still under co-regulation of products from nuclear-encoded genes (Palmer and Herbon, 1988; Burki, 2016). Mitogenomes of most land plants typically contain 40 to 60 genes that are inherited from a single parent (Palmer and Herbon, 1988). Compared to plastome, the mitogenome is usually larger in molecular weight. It varies widely in size, ranging from 208 kb (Brassica juncea) (Knoop, 2004) to over 2400 kb (Cucumis melo) (Ward et al., 1981), then to the largest known (Larix sibirica) of nearly 12 Mb (Putintseva et al., 2020). This is mainly due to the redundancy (multiple copies) of mitochondrial protein coding genes (PCGs) and the frequent recombination and integration of foreign DNA in the mitogenome (Mackenzie and McIntosh, 1999). Gene order, genome structure, and mitogenome size are highly variable in plants (Sloan et al., 2012a). The “evolutionary paradox” is an important evolutionary feature of plant mitogenome, that is, the mutation rate of plant mitogenome sequence is very low, but the rate of mitogenome structure rearrangement is really high (Palmer and Stein, 1986; Wolfe et al., 1987). The rate of nucleotide synonymous substitution in a plant mitogenome is several to dozens of times lower than that of the plastome and nuclear genome, and even 50 to 100 times lower than that of the mammalian mitogenome (Wolfe et al., 1987). Structural complexity is another important feature of plant mitogenomes. Physical mapping and sequencing of mitochondria in some species suggest that the complex structure of mitochondria is shaped by self-recombinant gene transfer to the nucleus and other unclear factors (Sloan, 2013). Structural analyses reveal frequent intramolecular and intermolecular recombination, resulting in a structurally dynamic assemblage of genome configurations, which makes the mitogenome a powerful model for studying genome structure and dynamic evolutionary patterns (Drouin et al., 2008). In addition, gene transfer is another kind of common event in plant mitogenomes which provides an essential contribution to the complexity of mitogenome. Intracellular gene transfer (IGT) is generally known as the mutual transfer of DNA between the two organelle genomes and the nuclear genome within the cell (Timmis et al., 2004). The mitogenome size variation largely reflects differences in the amount of non-coding content, which comes from diverse sources including repeats and large duplications (Alverson et al., 2011), IGT of nuclear and plastid DNA and horizontal gene transfer (HGT) from other species (Park et al., 2015).

In this study, we decoded the mitogenomes of three Broussonetia spp. (B. monoica, B. kaempferi, and B. papyrifera) based on a strategy of combining next-generation and third-generation sequencing methods. We comprehensively compared their genomic characteristics in terms of gene content, intracellular sequence transfer, mitogenomic synteny, organellar phylogeny and RNA editing sites, to investigate the intrageneric heterogeneity in terms of this maternally inherited organelle. These established mitogenomes will provide a valuable resource for future evolutionary and functional research on these traditional valuable taxa.

The plastome matrix was 54,143 bp in alignment size with 1200 sites with alignment gaps or missing data. A total of 4584 (8.47%) sites were variable, of which 3161 were singleton variable (3073 two-variants, 87 three-variants, and 1 four-variants) and 1423 were parsimony informative (1285 two-variants and 138 three-variants) ( Table S5 ). The mitogenome matrix was 13,251 bp in total with 388 sites with alignment gaps or missing data, 276 (2.08%) sites were variable which consisted of 197 singleton variable (196 two-variants and 1 three-variants) and 79 parsimony informative sites (75 two-variants and 4 three-variants) ( Table S5 ). The phylogenetic relationships revealed by plastome and mitogenome were identical here. In the plastome tree, all nodes were fully supported with the sole exceptional node supported by BS = 99, while in the mitogenome tree, there were two nodes with BS < 90. However, this would not influence greatly the definite phylogenetic relationship. Broussonetia was monophyletic with Allaeanthus as its sister. Within Broussonetia, B. papyrifera was well-supported to be sister to the B. kaempferi and B. monoica clade. Much shorter branch lengths in the mitogenome tree compared with the plastome tree indicated slower evolution of the mitogenome ( Figure 5 ).

The tissue-specific RNA editing site abundance and efficiency were compared in the Broussonetia spp. mitogenome and plastome. The mitogenome contained more than times of RNA editing sites than plastome did. Among interspecies, B. kaempferi had the richest (385, 104 in RA, and 281 in RB) while B. monoica contained the least (89) abundance, in terms of either type in mitogenome or plastome. The C>T was predominate accounting for 85to 97% in mitogenome and 65.3to 90.0% in plastome among the three species ( Figure 6A ). When it came to tissues within individuals, in B. kaempferi, mitochondrial sites were much more abundant in root (225, 58.4%) and plastid sites were more abundant in old stem (21, 42.9%). In the B. monoica mitogenome, there was no massive heterogeneity among the four tissues. But for the plastome, almost no tissues other than young leaf contained RNA editing sites. In B. papyrifera, the stipule possessed most sites in the both mitogenome (99, 52.9%) and the plastome (17, 73.9%) while inflorescence contained 58 (31.0%) in the mitogenome but only one (4.1%) in the plastome ( Figure 6C ). Interestingly, young leaf and young stem contained similar abundance in B. kaempferi and B. monoica, much more than those in B. papyrifera in each organelle genome ( Figures 6A, C ). RNA editing events tended to occur in different regions among three species ( Figure 6B ). In B. kaempferi, slight differences in location were observed in the mitogenome among the four tissues. In root and young stem, more sites were detected in cds across the plastome and the mitogenome, respectively. In B. monoica, huge differences could be found among tissues and organelle genomes. In root, more sites were located in cds showing a great gap with non-cds of the mitogenome while in the other three tissues, the gaps were not so huge, but in young leaf, the situation was reverse in the plastome. In B. papyrifera, more sites were in cds of mitogenome to various degrees except for young leaf with more sites in non-cds. RNA editing efficiency varied among species and tissues as well ( Figure 6D ). In B. monoica, high efficiency was shared in the mitogenome among four tissues especially in root and bud. In B. kaempferi mitogenome, only young stem contained over 50% sites with 0.8–1 efficiency as the stipule did in B. papyrifera mitogenome.

In addition, several nuclear mitochondrial DNA segments (NUMTs) could be actively transcribed as components of nuclear genome regardless of different tissues. In B. monoica, there were several PCGs and six intergenic spacers (IGS) participating in the expression process, of which two were annotated as the residue of rrn23 and rrn26, and four were open reading frames ( Figure 7 ). Similarly, five transcribed NUMTs from IGS of the mitogenome were detected in B. papyrifera, and three and two were in B. kaempferi RA and RB, respectively ( Figures S6–S8 ).

Plant mitogenomes have evolved from with dramatic and rapid structural changes, thus many unique features can be found compared with conserved plant plastomes and compact animal mitogenomes. These features include huge genome size variations, active multipartite arrangements, low gene density, abundant post-transcribed RNA editing, gene sequence transfer or loss, and foreign sequence capture (Hong et al., 2021; Wu et al., 2022). Although plant mitogenomes are conventionally reported as circular structures, much like the circular chromosomes found in animal mitochondria and bacteria, many studies have shown that the simple circular model of genome structure that applies to most animal species is inadequate when trying to understand plant mitochondria that can have multiple circular replicons, branched, linear, or mixed forms of genomic structure, i.e., Silene noctiflora has a mitogenome that is arranged into numerous circular chromosomes, that is, one ring cannot rule them all (Sloan, 2013; Wu et al., 2015) based on mapping to a circular reference. Here in this study, we obtained the first complete record of Broussonetia mitogenomes by providing three gap-free assemblies based on a hybrid strategy combining Illumina and ONT long sequencing reads. During de novo assembly, heterogeneity could be readily founded in configuration mediated by repeats that have potential to mediate structural transformation as all the connections were verified by PCR experimental methods ( Figure 1 ). It also showed that the three Broussonetia mitogenomes showed wide differences in structure. In particular, such a configuration consisting of double molecules was observed in the B. kaempferi mitogenome rather than the single molecule in the other two Broussonetia mitogenomes. A similar pattern was detected in okra (Abelmoschus esculentus) with circular and a liner-branching structure concurrently as well as Silene vulgaris with two independent molecules (Sloan et al., 2012b). In some cases, three or more multichromosomal circular mitogenomes could be observed like Cucumis sativus (Alverson et al., 2011) and Silene conica (Sloan, 2013). These observations suggest that multi-molecules may be more frequent than expected. Few genes shared between the two molecules indicated their necessary cooperation to function as normal. The polymorphisms in the conformation of the plant mitochondrial genome always puzzle us. Our results also supported the representation that plant mitogenomes should be considered a dynamic genome, at least, this structure is a more complete description of complicated plant mitogenomes (Li et al., 2022). The most complicated structure was observed in the B. papyrifera mitogenome because it has six repeats that could mediate 19 kinds of connections, but the B. monoica mitogenome had two repeats controlling seven kinds of connections. All these connections indicated an alternative configuration. To perform further analyses, a conventional single molecular structure was resolved in B. papyrifera and B. monoica of 325,822 and 276,967 bp in size while a double-ring structure was resolved in B. kaempferi of 267,420 bp in total size with 151,895 bp in the larger ring (RA) and 115,525 bp in the smaller one (RB). The varying size and structure may in part be the result of a limited number of repetitive sequences, which are known to be involved with large structural rearrangements (Abdelnoor et al., 2003; Hofmann, 2011), as the relative larger mitogenome size of B. papyrifera was due to a pair of 33-kb repeats ( Figures 1 , 2 ). During the evolutionary history, massive original mitochondrial genes have transferred into nuclear region for more stable inner environment or functional efficiency across different lineages (Wicke et al., 2011) and about 30 to 50 encoding genes remained. The three Broussonetia mitogenomes shared 36 unique genes including 23 PCGs, two rRNA genes, and 11 tRNA genes, and several tRNA genes were lost in B. kaempferi and/or B. papyrifera ( Table 1 , Figure 2 ), which was common in other species such as Dalbergia odorifera ( Hong et al., 2021 ) and Abelmoschus esculentus (Li et al., 2022). The conservation in gene content across many plant lineages indicated the stable function of the mitogenome as it is responsible for the basic energy synthesis for life activities through oxidative phosphorylation (Mackenzie and McIntosh, 1999). And obviously it is the evidence of “evolutionary paradox” of plant mitogenomes. The subsequent phylogenetic analysis based on organellar genes supported a closer relationship between B. kaempferi and B. monoica with B. papyrifera diverging the earliest within Broussonetia ( Figure 5 ), and that was consistent with the synteny pattern wherein more bulks of synteny could be found between these two rather than with B. papyrifera indicating a more similar sequence component and arrangement ( Figure 3 ).

RNA editing events are pretty frequent in plant mitogenomes [20]. RNA editing in functional genes can result in massive diversity in gene sequence beyond what is encoded at DNA level. These diversities may be influenced by indels and substitutions of specific nucleotides at the RNA level (Sloan and Wu, 2016; Wu et al., 2022). The number of RNA editing sites in land plant mitogenomes can vary dramatically from zero in Marchantia polymorpha to 2152 sites in Selaginella moellendorffii (Zhang et al., 2020b). Here in the three Broussonetia spp., massive heterogeneity in RNA editing site abundance, type, efficiency, and location could be detected among species and tissues. The B. kaempferi mitogenome contained the most abundant sites, especially in the root and old stem, and in B. monoica mitogenome, the stipule was the most site-rich tissue. The C>T was dominant of the type of no matter in any organellar genome, species, or tissue. In several tissues, more sites in cds tended to high-effectively editing ( Figure 6 ). Such massive heterogeneity among tissues might reflect the differences in activities based on functional divergence. Additionally, RNA editing may change the biochemical nature of the protein products by influencing functional structures. RNA editing events in Arabidopsis mitochondria showed that 35% of the modifications of the codons altered the amino acids from hydrophilic to hydrophobic. Also, RNA editing may contribute to protein stability and quantity (Jiang et al., 2022). It should be noted that although the average coverage was not that deep due to the low expression levels or a small amount of RNA sequencing data, it still provided plenty of information to further understand the potential functional roles that RNA editing has in the mitogenome. Given that RNA editing sites were not thought to exist in algae until recently, it suggests more attention to investigate RNA editing in a greater species range and that RNA editing is probably more common in later‐diverged land plants (Wu et al., 2022).

IGT as an intercellular type of horizontal gene transfer (HGT), widely exists in the plant genomes, occurs between organellar and nuclear genomes, and the DNA fragments from organellar genomes can be integrated into nuclear genomes to become NUMTs and NUPTs (Wang et al., 2012). Between the two plant organellar genomes, the most frequent direction is from plastome to mitogenome. Massive mitochondrial plastid fragments (MTPTs) could be detected in plant mitogenomes. In general, plant mitogenomes have about 0.56% (Marchantia polymorpha) – 10.85% (Phoenix dactylifera) plastid-derived sequences (Zhao et al., 2019). Here, we detected MTPTs in the three Broussonetia spp. and found that B. monoica, B. kaempferi, and B. papyrifera integrated 24 (30,604 bp), 27 (21,545 bp), and 26 (26,037 bp) MTPTs, respectively. They did not show great heterogeneity in MTPTs length distribution, and although B. papyrifera possessed the largest mitogenome size, its MTPTs were not the longest in length. These results indicate that the plastid-derived integrations could limitedly contribute to the mitogenome expansion in size. However, more analyses based on wide-ranged Broussonetia spp. are necessary to draw a precise conclusion. In addition, the integration of transferred fragments could be linked to DNA double-strand breaks (DSB) repair mechanisms such as homologous recombination (HR) and nonhomologous end joining (NHEJ), wherein several elements including transposable elements (TE), microhomology, and tRNA actively function (Zhao et al., 2019). Among the length of top-ranked regions, many tRNA genes could be detected, indicating they might take part in the integration of MTPTs.

Although the plant mitogenome has unresolved characteristics, it obviously has important functions for evolution and life activities. In this paper, we assembled three Broussonetia high-quality mitogenomes and performed a comprehensive comparison in terms of structure, gene content, codon usage, and IGTs to provide a detailed genome landscape. Single or double circular forms showed intrageneric complexity of configuration. Despite gene content being quite conserved the numbers of copies changed, especially in some tRNAs. IGTs tended to transfer from certain regions of plastomes but integrated randomly flanking with tRNAs across the mitogenomes. Additionally, the abundance of RNA editing sites was uneven in different species and tissues, but C>T sites were more frequent in general, and most RNA editing happened in CDS, together with a relatively higher editing efficiency. Intergenic nuclear mitochondrial DNA segments were transcribed at various levels, most of which were ORFs. Our analysis here will promote the understanding of plant mitogenome structure and evolution, which will benefit plant mitogenome evolutionary research in the future.

The data presented in the study are deposited in the GenBank repository, accession number OP342744, OP342745, OP342746, and OP342747.

YZ, ZW, JW and CL conceived and designed the study. CL, JW, and SZ performed the experiments and data analysis. PL, CL and YZ contributed materials collection. JW, and SK wrote the paper. ZW, PL, WR, and YZ revised the paper. All authors contributed to the article and approved the submitted version.

This study was supported by the Research Funds for the Natural Science Foundation of Zhejiang Province (Grant No. LY21C030002) and the Scientific Research Project of Baishanzu National Park (Grant No. 2021KFLY06) to YZ. And it also co-funded by the National Natural Science Foundation of China (31970244 and 32170238), and the Science, Technology, and Innovation Commission of Shenzhen Municipality (RCYX20200714114538196) to ZW.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.