Comprehensive analysis of the Lycopodium japonicum mitogenome reveals abundant tRNA genes and cis-spliced introns in Lycopodiaceae species

Lycophytes and ferns represent one of the earliest-diverging lineages of vascular plants, with the Lycopodiaceae family constituting the basal clade among lycophytes. In this research, we successfully assembled and annotated the complete Lycopodium japonicum Thunb. (L. japonicum) mitochondrial genome (mitogenome) utilizing PacBio HiFi sequencing data, resulting in a single circular molecule with a size of 454,458 bp. 64 unique genes were annotated altogether, including 34 protein-coding genes, 27 tRNAs and 3 rRNAs. It also contains 32 group II introns, all of which undergo cis-splicing. We identified 195 simple sequence repeats, 1,948 dispersed repeats, and 92 tandem repeats in the L. japonicum mitogenome. Collinear analysis indicated that the mitogenomes of Lycopodiaceae are remarkably conserved compared to those of other vascular plants. We totally identified 326 RNA editing sites in 31 unique protein-coding genes with 299 sites converting cytosine to uracil and 27 sites the reverse. Notably, the L. japonicum mitogenome has small amounts foreign DNA from plastid or nuclear origin, accounting for only 2.81% of the mitogenome. The maximum likelihood phylogenetic analysis based on 23 diverse land plant mitogenomes and plastid genomes supports the basal position of lycophytes within vascular plants and they form a sister clade to all other vascular lineages, which is consistent with the PPG I classification system. As the first reported mitogenome of Lycopodioideae subfamily, this study enriches our understanding of Lycopodium mitogenomes, and sets the stage for future research on mitochondrial diversity and evolution within the lycophytes and ferns.

The four existing mitogenomes have yielded significant insights into the degree of mitochondrial genomic diversity across these lineages.For instance, the group I and group II introns were found in S. moellendorffii and I. engelmannii, several of which are trans-splicing introns, while the introns in the mitogenomes of P. squarrosus and H. crispata are cis-splicing and all belong to group II introns.Although our understanding of the phylogenetic relationships of lycophytes and ferns has developed, related studies are still lacking compared with that of seed plants.Consequently, it is urgent to sequence the complete mitogenome from additional Lycopodiaceae members to elucidate the evolutionary status of basal vascular plants and clarify the phylogenetic relationships with the Lycopodiaceae species.
In this research, we conducted the assembly and annotation of the complete Lycopodium japonicum Thunb.(L.japonicum) mitogenome, which belongs to Lycopodioideae subfamily.The genome composition, intron content, repeats, RNA editing sites, as well as mitochondrial plastid DNA (MIPT) (Mower et al., 2012) of L. japonicum were analyzed, filling a gap of the mitogenome data in Lycopodium.Furthermore, we analyzed the mitogenome using a comparative genomics approach and reconstructed the phylogenetic trees according to protein-coding genes (PCGs).These findings will offer valuable insights for the molecular classification, identification, and germplasm conservation of Lycopodioideae plants and enhance our understanding of the evolutionary significance of lycophytes.

Plant materials and sequencing
We gathered fresh leaves of L. japonicum from Kunming Institute of Botany, Chinese Academy of Sciences and stored them at -80°C for subsequent use.We utilized Hi-DNA secure Plant Kit method to extract High-quality total genomic DNA.The quality and purity were assessed using a Nanodrop 2000 spectrophotometer (ThermoFisher) and a 1.0% agarose gel electrophoresis.Subsequently, we constructed sequencing libraries using the high-integrity genomic DNA and the SMRTbell Express Template Prep Kit 2.0 (PacBio Biosciences, CA, USA).Ultimately, the HiFi sequencing data was generated using the PacBio Revio platform.

Assembly and annotation of the L. japonicum mitogenome
The HiFi sequencing data was fed into PMAT v1.31 (Bi et al., 2024a) to assemble the mitogenome of L. japonicum.The parameters and mode used for PMAT were '-st hifi -g 2.3G -cpu 50' and 'autoMito' mode, respectively.The genome size of L. japonicum was estimated using Lycopodium clavatum as the reference (Yu et al., 2023).Bandage was then used to visualize and disentangle the raw assembly graph of L. japonicum mitogenome (Wick et al., 2015).Online program IPMGA (http://www.1kmpg.cn/ipmga/)was used to annotate the L. japonicum mitogenome.Further, tRNAs and rRNAs were verified by BLASTn and tRNAscan-SE v2.0, respectively (Camacho et al., 2009;Chan et al., 2021).Afterwards, the introns were designated by their orthologous gene positions within Marchantia polymorpha mitogenome (Dombrovska and Qiu, 2004).All the annotations were manually checked and reviewed using MacVector v18.5.The mitogenome map and the intron contents of the L. japonicum was visualized using PMGmap with default parameters (Zhang et al., 2024).
2.4 Detection of RNA editing events in L. japonicum mitogenome RNA editing, referring to the insertion, deletion, or substitution of bases in genomic transcripts, plays significant roles in increasing the diversity of gene transcription and function (Fan et al., 2019).LncRNA sequencing data were generated from the same sample of L. japonicum to detect the RNA editing events in the L. japonicum mitogenome.We extracted coding sequences of each PCG with 100 bp upstream and downstream to construct reference sequences (Yang H. et al., 2023).Next, we used BWA to map strand-specific RNA sequencing data to reference sequences (Li and Durbin, 2009), with parameter settings 'mem -t 40 -k 55 -d 150 -r 1.0 -M'.Sambamba was used to sort BAM files and to filter secondary or supplementary alignments (Tarasov et al., 2015).REDItools v2 was utilized to detect the RNA editing sites according to mapping results (Picardi and Pesole, 2013).We employed minimap2 to align the HiFi sequencing data to the L. japonicum mitogenome (Li, 2018;2021), and subsequently extracted sites that were not recognized as SNPs, which indicated potential RNA editing locations.IGV program was finally used to manually verify and validate all possible RNA editing sites (Robinson et al., 2011;Kim et al., 2019).Subsequently, we also used HISAT2 to identify RNA editing events and compared the results with those of REDItools.

Analysis of MIPTs
Higher plant mitogenomes have numerous sequences migrating from their plastomes.To detect the MIPTs of L. japonicum, we assembled the complete L. japonicum plastome (accession number: NC_085262).Subsequently, we employed BLASTn to detect the homologous fragments between the mitogenome and plastome with parameters '-word_size 9 -evalue 1e-5 -reward 2 -gapopen 5gapextend 2 -penalty -3 -outfmt 6'.The results with lengths ≥50 bp and matching rate ≥70% were selected for further analysis and were manually annotated to check the gene location in the MIPTs.Ultimately, we visualized the MIPT results utilizing TBtools v1.132 with 'Advanced Circos' module.
We detected 32 introns within 15 PCGs (atp6, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad5, rpl2, rps3, rps10, and rps14), and all introns are categorized as cis-splicing introns (Supplementary Figure S2).Mitochondrial introns are classified into two types for the majority of eukaryotes, group I and group II, based on the mechanisms by which they splice and the secondary structures (Michel and Westhof, 1990;Saldanha et al., 1993).In the mitogenome of L. japonicum, all 32 introns are group II introns.To elucidate the relatively conserved pattern of intron evolution, we compared cis-and trans-spliced intron content in mitogenomes Assembly graphs and genome map of the mitogenome.S1C).

Analysis of repeats and MIPTs
In mitochondrial DNA, the dynamic restructuring of its architecture was largely mediated by the presence of repeats, including SSRs, dispersed repeats and tandem repeats.A total of 195 SSRs were detected in L. japonicum using web-based platform MISA, which were composed of 132 monomeric, 19 dimeric, 7 trimeric, 30 tetrameric, 5 pentameric, and 2 hexameric repeats (Figure 2A).It was observed that the repeat units of A/T (94 repeats), AT/AT (nine repeats), AAT/ATT (five repeats), and AAAG/CTTT (eight repeats) were more prevalent in the monomeric, dimeric, trimeric, and tetrameric categories, respectively.Additionally, 92 tandem repeats were identified (Figure 2B), with nearly half of these ranging between 50 and 99 bp in length, while only two exceeded 500 bp (651 bp and 650 bp, respectively).We identified 1,949 dispersed repeats in L. japonicum mitogenome with length ≥30 bp (total length: 124,620 bp, account for 27.42% of the whole mitogenome), including 1,001 forward repeats, 943 palindromic repeats, and only five reverse repeats.The majority of the dispersed repeats are shorter than 100 bp, with a count of 1,811 repeats, and only 10 repeats exceed 500 bp, with the maximum length being 16,967 bp (Figures 2C, D).
The mitogenomes of higher plants contain numerous sequences that migrated from nuclear and plastid genomes (Mower et al., 2012).In this study, we identified 39 MIPTs (12,767 bp in length) with a similarity threshold of 70%, which represent 8.40% of the plastome and 2.81% of the mitogenome, respectively (Figure 3).Five MIPTs are more than 1,000 bp in length, with four of them uniformly measuring 1,165 bp.However, among these homologous sequences, only partial plastid genes are annotated locating within the MIPTs, including atpA, ndhF, rps11, rrn16S, rrn23S, trnN-GUU, and trnF-GAA.S2 and S3, respectively. of the entire mitogenome of H. crispata.A comparison with P. squarrosa revealed LCBs (total length: 395,656 bp), which constitutes 95.68% of the mitogenome of P. squarrosa (Figure 4).It is observed that compared with other land plants mitogenomes, the number of LCBs, average length and collinear identity score, were consistently lower than those observed in P. squarrosa and H. crispata.Additionally, numerous sequence rearrangements were identified among the L. japonicum mitogenome and those of other species.

Analysis of RNA editing events
To estimate the extent of RNA editing sites in L. japonicum mitogenome, the BWA program was employed to identify RNA editing sites based on 14.57Gb RNA sequencing data.After manually checking using IGV and comparison with SNPs, we finally identified 326 RNA editing sites among 31 PCGs of L. japonicum mitogenome.Notably, most of these sites (299, accounting for 91.72%) involved C-U substitutions, while 27 sites (8.28%) exhibited U-C substitutions (Supplementary Table S5).Specifically,21 genes (atp1,atp9,cob,cox1,mttB,nad1,nad2,nad4L,nad9,rpl2,rpl5,rpl6,rpl10,rpl16,rps11,rps2,rps3,rps12,rps13,rps19, and sdh4) exclusively contain C-U RNA editing sites, while 10 genes (atp4, atp6, cox2, cox3, nad3, nad4, nad5, nad6, rps10, and sdh3) exhibit both C-U and U-C RNA editing sites, with cox1 containing the highest number of RNA editing sites, totaling 40 (Figure 5A).Notably, three genes (matR, rps4, and rps14) were found to lack any RNA editing sites.Out of the total RNA editing sites detected, 79 (24.23%) sites were located in the first codon position, 204 (62.58%) in the second codon position, and 43 (13.19%) in the third codon position.We observed that RNA editing events had reconstituted six start codons and five stop codons in 11 PCGs (Supplementary Table S1).Subsequently, we analyzed the synonymous and non-synonymous substitutions caused by RNA editing events (Figure 5B).We detected 33 substitution types (10 kinds of synonymous and 23 kinds of non-synonymous substitutions) in the L. japonicum mitogenome.The most frequent substitution involved the transformation of 82 amino acids from Ser to Leu.

Phylogenetic analysis
The ML tree was constructed based on 15 conserved mitochondrial PCGs from 23 plant species to investigate the phylogenetic position of L. japonicum.We selected two algal species as outgroups.The phylogenetic analysis showed that  S4.
L. japonicum clustered together with the other two Lycopodiaceae P. squarrosus and H. crispata (Figure 6).The three Lycopodiaceae species are a sister group of S. moellendorffii + I. engelmannii.The reconstructed ML tree demonstrates a wellsupported phylogeny of lycophytes (bootstrap values = 100), and the entire topology displays a high degree of concordance with the PPG I classification system (PPG I, 2016).We further validated our phylogeny by reconstructing the phylogenetic tree using plastome sequences.The structural comparison of the ML trees derived from sequences of mitogenomes and plastomes reveals a perfect consistency, validating the precision of our phylogenetic analyses.

Discussion
4.1 The moderate size of the L. japonicum mitogenome Previous studies that put forward the hypothesis that the typical size of the vascular plant mitogenome is around 400 kb (Guo et al., 2017).In this study, the mitogenome of L. japonicum is structured as a circular DNA molecule with a GC content of 44.4% and a size of 454,458 bp.The mitogenomes of three Lycopodiaceae species (L.japonicum, P. squarrosus, and H. crispata) are relatively conserved.Each of these mitogenomes exhibits a typical single circular structure and varies slightly in length, ranging from 412,594 bp (H.crispata) to 454,458 bp (L.japonicum).The mitogenomes of ferns, however, exhibit variability in structure, GC content, and length when compared with those of Lycopodiaceae mitogenomes.For instance, the mitogenome size of Dryopteris crassirhizoma is 313 kb in length with a GC content of 49.0% (Song et al., 2021).The Ophioglossum californicum (O.californicum) mitogenome exhibits a larger mitogenome of 372 kb and a higher GC content of 52.2% (Guo et al., 2017).The Psilotum nudum (P.nudum) mitogenome is composed of two circular molecules with respective sizes of 264 kb and 364 kb, having a total GC content of 51.2% (Guo et al., 2017;Song et al., 2021).
Studies conducted previously revealed that the mitogenomes of ferns contain abundant repeats, such as P. nudum (52.7-63.3%)and O. californicum (37.1-44.3%)(Guo et al., 2017).In comparison to P. nudum and O. californicum, the Lycopodiaceae mitogenomes exhibit a notable decrease in the proportion of repeats.Specifically, in the L. japonicum mitogenome, repeats make up 27.42% of the entire mitogenome, a proportion comparable to that of P. squarrosus (26.36%).Additionally, the results of MIPT analysis indicate that the L. japonicum mitogenome has limited foreign DNA from plastid or nucleus, which is a characteristic shared with the mitogenomes of S. moellendorffii and P. squarrosus (Liu et al., 2012;Wu et al., 2017).However, the foreign DNAs from plastid or nuclear are detected in the mitogenomes of H. crispata, I. engelmannii, Cycas taitungensis and numerous angiosperms (Troia et al., 2016;Cao et al., 2023).It seems that the mitogenome of L. japonicum maintains a stable size, probably because it has not experienced the substantial invasion of foreign DNA from plastid or nucleus, a phenomenon seen in certain angiosperm mitogenomes (Giegéet al., 2008;Babbitt et al., 2015).

FIGURE 4
Whole mitogenome collinearity analysis.Long strips of different colors represent different plant mitogenomes.Linear blocks represent collinear segments.Species names and IDs used in collinearity analysis are shown in Supplementary Table S7.Sun et al. 10.3389/fpls.2024.1446015Frontiers in Plant Science frontiersin.org4.2 The loss of ccm genes and abundant tRNAs in the of Lycopodiaceae Ccm genes are typically present in the mitogenomes of most gymnosperms and angiosperms; however, they are absent or present as pseudogenes in all lycophytes mitogenomes.In the L. japonicum mitogenome, three ccm genes (ccmB, ccmC, and ccmFn) are absent, with only ccmFc present as a pseudogene.The loss of ccm genes is also found in the mitogenomes of P. squarrosus and H. crispata (Liu et al., 2012;Cao et al., 2023).These ccm genes are absent in some fern mitogenomes like Dryopteris crassirhizoma and O. californicum, but present in other fern species such as P. nudum (Troia et al., 2016;Guo et al., 2017;Wu et al., 2017;Song et al., 2021) (Supplementary Figure S1A).The aforementioned results suggest that the pathway of cytochrome c maturation seems to have been inherited from the original proto-mitochondrion but it has been lost in Lycopodiaceae species and some ferns (Babbitt et al., 2015).
Protein synthesis in plant mitogenomes necessitates 21 different tRNAs, but some tRNAs are lost frequently during the evolutionary process of higher plant mitogenomes.For example, S. moellendorffii mitogenome lacks all tRNAs, and I. engelmannii only maintains 13 intact tRNAs (Troia et al., 2016;Wu et al., 2017).In this study, a total of 27 unique tRNAs were annotated in the mitogenome of L. japonicum.Similarly, the mitogenomes of P. squarrosus and H. crispata contain 26 and 29 unique tRNAs, respectively (Liu et al., 2012;Cao et al., 2023).The findings suggest that the tRNA content in Lycopodiaceae mitogenomes is much higher than that of other species.

The loss of introns during the evolution of land plants
Studies conducted previously showed that it is challenging to infer the ancestral intron contents of vascular plants due to the  S5.
significant differences in intron distribution observed among lycophytes, ferns, seed plants (Guo et al., 2017).In the mitogenome of L. japonicum, a total of 32 group II introns were detected, all of which are cis-splicing introns.Moreover, within the mitogenomes of Lycopodiaceae species, the number of introns is very abundant, and all these introns undergo cis-splicing (Supplementary Figure S1B).The results strongly support that the mitogenomes of Lycopodiaceae species exhibit low levels of recombination and are highly conservative.In the lineage of land plants, the mitogenomes contain a high number of group II introns, whereas group I introns are much less common and display an irregular pattern of distribution across different species (Mower, 2020).The evolutionary genesis of the five different group I introns (cox1i375, cox1i395, cox1i624, cox1i876, and cox1i1305) found in either ferns or lycophytes remains ambiguous, given their highly specific distribution patterns within the vascular plant lineage.To elucidate whether the group I introns have been acquired through vertical or horizontal gene transfer, it is essential to undertake a more comprehensive study that includes a wider array of taxa.These results would enhance our understanding of the evolutionary dynamics governing these genetic elements within the plant mitogenomes.

Low RNA editing level in the mitogenomes of Lycopodiaceae species
RNA editing event, which involves the post-transcriptional modification of an RNA sequence through substitutions, deletions, or insertions, contributes significantly to transcriptome diversity (Lukešet al., 2021).It is likely that RNA editing events take place in the L. japonicum mitogenome, given that 11 PCGs necessitate the reconstruction of start or stop codons for accurate annotation (Supplementary Table S1).In this study, 326 RNA editing sites (299 C-U sites and 27 U-C sites) were detected in the L. japonicum mitogenome based on 14.57Gb lncRNA sequencing data.Previous reports demonstrated that RNA editing events vary widely across lycophytes mitogenomes.Specifically, the amount of RNA editing sites varies among different species, with P. squarrosus exhibiting a relatively low count of 364 sites, while I. engelmannii and S. moellendorffii have significantly higher numbers, with 1,782 and 2,152 sites, respectively (Liu et al., 2012;Troia et al., 2016;Wu et al., 2017).It can be deduced that the mitogenomes of Lycopodiaceae species exhibit low RNA editing level according to the given information.Furthermore, we utilized HISAT2 to validate these RNA editing sites.The findings indicated that HISAT2 and bwa share identical 272 C-U RNA editing sites.It is noteworthy that no U-C editing sites were detected using HISAT2 (Supplementary Tables S5 and S6).Although all RNA editing sites were manually verified using IGV, it is imperative to conduct Sanger sequencing and PCR to acquire a more precise and accurate result.

Phylogenetic analysis of L. japonicum
Based on PPG I, two pteridophyte classes are recognized, including Lycopodiopsida (lycophytes) and Polypodiopsida (ferns) (PPG I, 2016).The fossil record indicates that lycophytes and ferns are the oldest vascular plants on Earth, dating back to approximately 400 million years ago (Qi et al., 2018;Shen et al., 2018).As a result, the mitogenomes of lycophytes represent the most ancient form of vascular plant mitogenomes.In this study, we reconstructed the phylogenetic trees of 23 plant species based on 15 conserved mitochondrial PCGs and whole plastid genome sequences, The ML trees of 23 plant species based on 15 mitochondrial PCGs (left) and complete plastome sequences (right).The Lycopodium japonicum mitogenome is bold and labeled with an asterisk.Saccharina japonica and Laminaria rodriguezii were selected as outgroup.Bootstrap values are marked on each branch.Colors indicate the groups for each species.Plant names and NCBI accession numbers used in the phylogenetic analysis are provided in Supplementary Table S7.
separately.Both phylogenetic trees support that L. japonicum is a group to the clade clustered with P. squarrosus and H. crispata, thereby corroborating the familial and generic relationships among the three species.Additionally, the results of phylogenetic analysis also suggest that lycophyte are the basal group of vascular plants, in accordance with the PPG I classification system.Our results conclusively demonstrate that lycophytes are the oldest vascular plants and serve as transitional evolutionary forms bridging the gap between bryophyte and higher vascular plants (Qiu et al., 2006;Guo et al., 2017).Furthermore, the bootstrap values of most nodes are greater than 95%, indicating the robustness and reliability of the recovered phylogeny based on well-conserved mitochondrial PCGs.

Conclusion
In this study, we successfully assembled the L. japonicum mitogenome into a circular molecule with a size of 454,458 bp.We annotated 64 unique genes in the L. japonicum mitogenome, including 34 PCGs, 27 tRNAs and 3 rRNAs.32 introns were identified in L. japonicum mitogenome and all of them are cissplicing.The MIPTs analysis revealed the absence of any complete plastid genes in the L. japonicum mitogenome.Collinearity results indicated that the Lycopodiaceae species mitogenomes are highly conserved in terms of structure and size.Additionally, we detected 326 RNA editing sites in the L. japonicum mitogenome, indicating that L. japonicum exhibits a low RNA editing level in lycophytes mitogenomes.The reconstructed ML trees of mitogenome and plastome strongly support that L. japonicum forms a sister group to the clade containing by P. squarrosus and H. crispata, while lycophytes are identified as the basal group of vascular plants.Moreover, we analyzed the repeats of the L. japonicum mitogenome and explored the pattern of intron evolution from bryophytes to angiosperms.As the inaugural reported mitogenome in the Lycopodioideae subfamily, this study provides invaluable information for the assembly of other species within the same subfamily and contributes insightful perspectives on intron evolution, RNA editing, phylogeny, and mitochondrial evolution in lycophytes.
(A) Raw assembly graph.(B) Disentangled assembly graph.(C) The genome map of the Lycopodium japonicum mitogenome.Different types of genes are depicted using different color blocks.
FIGURE 2 Repeat elements identified in the Lycopodium japonicum mitogenome.(A) The frequency of identified SSRs.(B) The frequency of identified tandem repeats.(C) The frequency of identified dispersed repeats.(D) Distribution of repeat elements in the mitogenome.The outer circle represents the tandem repeats, followed by the SSRs, and the inner lines represent the dispersed repeats.Detailed information of dispersed and tandem repeats is shown in Supplementary TablesS2 and S3, respectively.

FIGURE 3
FIGURE 3 Distribution of MIPTs in Lycopodium japonicum.The blue and green arcs represent the plastid and mitochondrial genomes of L. japonicum, respectively.The outer circle represents the GC content of the two genomes, followed by the bars representing the length of MIPTs.Lines between the two arcs represent the MIPTs.Detailed information of MIPTs is shown in Supplementary TableS4.
FIGURE 5 RNA editing events in the Lycopodium japonicum mitogenome.(A) Characteristic of RNA editing sites across all PCGs in the mitogenome.The green bars indicate C-U RNA editing sites while the green bars are U-C RNA editing sites.The quantity of both is marked above or inside the bar.(B) The synonymous and non-synonymous substitutions caused by RNA editing.The color and size of the circle represent the type of raw amino acid and the number of the substitution, separately.The asterisk on the vertical axis represents stop codons.Detailed information of RNA editing sites is shown in Supplementary TableS5.

FIGURE 6
FIGURE 6 Sun et al. 10.3389/fpls.2024.1446015Frontiers in Plant Science frontiersin.orgfrom bryophyte to seed plants.As shown in Supplementary Figure S1B, P. squarrosus and L. japonicum share 32 identical group II introns, with H. crispata containing 31 group II introns and only missing rps10i235.In the mitogenomes of five lycophytes, certain introns such as atp6i439, atp9i21, atp9i87, cobi693, and cobi787 are present, yet these are conspicuously absent in many other vascular plants.Among the 21 species examined, group I introns were exclusively found in the mitogenomes of I. engelmannii, S. moellendorffii, Psilotum nudum, and three bryophytes, but were notably absent from the mitogenomes of other vascular plants (Supplementary Figure