Three algal samples including Ulva flexuosa (Ufl1), Ulva prolifera (Upr1), and Ulva intestinalis (Uin) collected in China were transported to laboratory in coolers (5–8°C) after sampling (Supplementary Table 1). Fresh tissue from one individual thallus was used for DNA extraction using a Plant Genomic DNA Kit (Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. Four algal samples including Ulva torta (Uto), U. flexuosa (Ufl2), Ulva sp. TM637 (Usp1), and Ulva sp. TM708 (Usp2) from the United States were preserved in silica gel after collection and as herbarium vouchers, which were submitted to the University of Alabama Herbarium (UNA) (Supplementary Table 1). The algal thallus (UTEX LB 1423) of Percursaria percursa (Ppe) was from the UTEX Culture Collection of Algae¹. Dried tissue from one individual thallus for United States samples was used to extract DNA with a Qiagen Plant DNA Extraction Kit (QIAGEN, Valencia, CA, United States). Species identification was performed based on phylogenetic analyses of two common marker datasets (the nuclear ITS region and the chloroplast rbcL gene) (Liu et al., 2013).

The concentration and quality of isolated DNA were measured with a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). For the samples from China, the purified DNA was fragmented into 350 bp and used to construct short-insert libraries. The short fragments were sequenced using an Illumina HiSeq 4000 sequencing platform. For the samples from the United States, paired end reads (150 bp) were sequenced at Cold Spring Harbor Laboratory on an Illumina MiSeq platform. Poor quality sequences and sequencing adapters were removed using Trim Galore! v0.3.7. Eight ulvalean mitogenomes were constructed using a combination of de novo and reference-guided assemblies. The mitogenome of U. prolifera (KT428794) was used as the reference genome for assembly. Mitogenome assembly was done with both A5 (Tritt et al., 2012) and Geneious R7 (Kearse et al., 2012). The mtDNA assembly was examined using the MEM algorithm of BWA v0.7.17 (Li and Durbin, 2010), and the mutation sites were verified using VarScan v2.3.9 (Koboldt et al., 2009). Incomplete genomes were closed by iteratively mapping the trimmed reads on to the contig with the Geneious 7.1 software (Biomatters²).

Protein-coding genes (PCGs) were annotated by Open Reading Frame (ORF) Finder at the National Center for Biotechnology Information (NCBI) website³, DOGMA (Wyman et al., 2004) and ORF finder in the Geneious 7.1 software. Transfer RNA genes (tRNAs) were searched for by reconstructing their cloverleaf structures using the tRNA scan-SE 1.21 software with default parameters (Chan et al., 2021). The large and small subunit ribosomal RNA genes (rRNAs) were identified by comparing newly sequenced ulvalean mtDNAs with homologous rRNA genes from the known Ulva mtDNAs deposited in the GenBank database, respectively (Table 1). The 5S rRNA gene (rrn5) was found by the RNAweasel Tool⁴ and RNA Folding (Zuker, 2003). Intron insertion-sites were identified manually based on the alignments of nucleotide (nt) sequences for homologous genes with or without introns from these 33 ulvalean mitogenomes. The corresponding genes (atp1, cox1, cox2, nad3, nad5, rnl, and rns) in the U. compressa (KY626327) (Uco3) mitogenome were used as a reference (Liu F. et al., 2020). Intron name was defined as host gene plus insertion site. The free-standing and intronic orfs greater than 300 bp were found by Open Reading Frame Finder at the NCBI website. The class and core structure of all these introns were determined using the software RNAweasel and RNA Folding. We observed that some annotations were incomplete or incorrect in the Ulva mitogenome data deposited in the GenBank database. In order to ensure the accuracy of our comparative analysis, we re-annotated all of the deposited Ulva mitogenomes with the same method.

The aa sequences of DNA-dependent RNA polymerase genes (rpos) were searched based on three full-length rpo genes found in Ulva mitogenomes as a query dataset through the BLAST program at the NCBI website⁵. A total of 45 closest Rpo proteins were obtained from the GenBank database for phylogenetic analysis. The core Rpo domains in these 48 Rpo proteins were determined by significant Pfam-A matches (Punta et al., 2012). Multiple sequence alignments of core Rpo domains were conducted using ClustalX 1.83 with the default settings (Thompson et al., 1997). To avoid phylogenetic artifacts caused by convergent base composition induced by synonymous substitutions (Cox et al., 2014), the ML phylogenetic tree was constructed based on aa sequences of core Rpo domains with 1,000 bootstrap replicates using MEGA 7.0 (Kumar et al., 2016). The phylogenetic relationships were inferred based on the Jones et al. w/freq. model (Jones et al., 1992). Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using a JTT model. There was a total of 902 positions in the final dataset of Rpos.

In group IIA/IIB introns of Ulva organelle genomes, the intron-encoded protein (IEP or intronic orf) was one reverse transcriptase/maturase (RTM). Three novel group IIA introns including intron atp1-1095, cox1-312, and nad5-1057 have been found for the first time in Ulva mitogenomes. To further understand the relationships between the new RTMs and those previously reported (Liu F. et al., 2020; Liu and Melton, 2021), we performed the phylogenetic analysis of the conserved reverse transcriptase (RT) domains which were determined by significant Pfam-A matches (Punta et al., 2012). The amino acid (aa) sequences of 94 RT domains (54 in mtDNAs and 40 in cpDNAs) from group IIA/IIB introns were subjected to concatenated alignments using ClustalX 1.83 with the default settings (Thompson et al., 1997). Maximum Likelihood (ML) phylogenetic tree was constructed for the RT dataset based on the Jones et al. w/freq. model (Jones et al., 1992) with 1,000 bootstrap replicates using MEGA 7.0 (Kumar et al., 2016). Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using a JTT model. There was a total of 330 positions in the final dataset of RT domains.

Base composition of these 33 ulvalean mtDNAs was determined using MEGA 7.0 (Kumar et al., 2016). The nt sequences of 61 genes (including 29 PCGs, three rRNAs, 26 tRNAs, and three conserved orfs) and the aa sequences of 32 genes (including 29 PCGs and three orfs) were subjected to concatenated alignments using ClustalX 1.83 with the default settings, respectively (Thompson et al., 1997). For the nt sequence dataset of the 61 genes, the evolutionary history was inferred by using the ML method based on the Tamura-Nei model (Tamura and Nei, 1993). Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach. For the aa sequence dataset of the 32 genes, the evolutionary history was inferred by using the ML method based on the Jones et al. w/freq. model (Jones et al., 1992). Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using a JTT model. There was a total of 39,070 and 10,491 positions in the final nt and aa datasets, respectively. Phylogenomic analysis was conducted with 1,000 bootstrap replicates using MEGA 7.0 (Kumar et al., 2016).

Eight newly sequenced ulvalean mitochondrial genomes (seven Ulva mtDNAs and one Percursaria mtDNA) were acquired in this study and compared with 25 known Ulva mitogenomes to understand the evolution of Ulva mtDNAs. The Ulva mitogenomes ranged in size from the smallest one, 55,814 bp in Ulva sp. TM708 (Usp2), to the largest one, 88,416 bp in U. rigida (Uri) (Table 1). The Usp2 mtDNA is the second smallest mitochondrial genome in Ulvophyceae to date, which lies between the 45,971-bp mitogenome of Codium fragile (Bryopsidales) and the 56,761-bp mitogenome of Oltmannsiellopsis viridis (Oltmansiellopsidales) (Pombert et al., 2004, 2006). All ulvalean mitogenomes were biased toward A + T nucleotides and the A + T content ranged from 61.16% in U. compressa (Uco6) to 67.83% in Ulva sp. (Usp3).

These 33 ulvalean mitogenomes shared a total of 62 core genes including 29 PCGs, three rRNAs (rnl, rns and rrn5), 26 tRNAs, three conserved free-standing orfs, and one putative RNA subunit of RNase P (rnpB) (Table 2). The homologs of these three conserved orfs share identical positions in all sequenced ulvalean mtDNAs, but based on blastp and Pfam, we can not determine their function. In mitogenome of U. rigida (Uri), the conserved orf located between rps2 and trnL1(uaa) was split into two orfs, orf326 and orf199 (Figure 1). All these core genes are coded on the same strand, and have the identical gene order in these ulvalean mitogenomes with the exception of one U. compressa mitogenome (Uco1), in which a collinear block of eight genes (rps11-rps19-rps4-rpl16-trnR3-trnQ-trnE-trnS3) has been inverted (Liu F. et al., 2020).

It is worth noting that the mitochondrial 5S rRNA gene (rrn5) previously unrecognized in Ulva mitogenomes was situated between trnG(ucc) and trnW(cca) in all sequenced ulvalean mitogenomes, while it was flanked by rnl and trnI(gau) in ulotrichalean mtDNAs. The sequences of domain β were conserved among the mitochondrial 5S rRNAs in Ulvales, but were different from the counterparts in Ulotrichales. Similar to that in mtDNAs of brown algae and some Ochrophyta lineages, the mitochondrial rrn5 gene in Ulvales was folded into a secondary structure by adopting a permuted triskelion shape (Valach et al., 2014; Figure 2A). We observed that the 5S rRNAs in ulvalean mtDNAs displayed a much larger structural variability than these in other lineages. The interior loop B in stem C is very asymmetric (Figure 2B), which makes it difficult to recognize the rrn5 genes. The sequences (approximately 35 bp) of domain γ including helix V, loop E, helix IV and loop D were highly conserved among rrns genes in not only Ulvales but also Ulotrichales. However, the homologous sequence was not detected in the mitogenomes of Bryopsidales and Oltmannsiellopsidales, as might be due to the sequence divergence, compositional bias and/or structural deviation of mitochondrial 5S rRNA genes (Valach et al., 2014).

The putative RNA subunit of mitochondrial RNase P (rnpB) could be detected in ulvalean mitogenomes using RNAweasel, which was situated in the latter part of the conserved orf (e.g., orf500 in U. torta) between rps2 and trnL1(uaa). The rnpB gene is encoding the RNA component of RNase P which participates in the generation of mature 5′-ends of tRNAs by cleaving the 5′-leader elements of tRNA precursors (Daoud et al., 2012; Shaukat et al., 2021). The highly variable sequence of rnpB makes it difficult to be found in mtDNAs. Thus far, the mitochondrial rnpB gene was found only in mtDNAs of the prasinophyte green alga Nephroselmis olivacea (Turmel et al., 1999), some ascomycete fungi (Seif et al., 2005) and jakobid protists (Lang et al., 1997; Burger et al., 2013).

The number of tRNA genes in Ulva mitogenomes showed slight differences at the interspecific level, ranging from 27 to 30 (Table 2). A total of 26 core tRNA genes are highly conserved in terms of both composition and structure, and shared by all sequenced ulvalean mitogenomes. The secondary structure of most tRNAs showed typical clover structures with the exception of the trnS3(cga) gene located between trnE(uuc) and nad5. It seems that trnS3(cga) has undergone insertion or deletion mutations in Ulva mtDNAs, causing the change of their secondary structures, or even severe degradation in the DHU arm and DHU loop, which led to the loss of this gene in the mtDNA of U. rigida (Uri). The trnX1 flanked by cox3 and rpl14 displayed variation in the anticodon loop sequences, which caused a marked change in function among different Ulva mtDNAs (Liu F. et al., 2020). This tRNA gene has been remolded into trnK2(uuu) in U. flexuosa (Ufl1-3) and U. meridionalis (Ume); or trnP2(ugg) in U. prolifera (Upr1-3) and U. linza (Uli); or trnI2(uau) in U. compressa (Uco1-3, 5, and 6); or trnV2(uac) in Ulva sp. TM708 (Usp2), U. fenestrata (Ufe), U. expansa (Uex), and U. rigida (Uri). This gene was lost in the mtDNAs of U. compressa (Uco4) and U. australis (Uau1-3) (Liu et al., 2017; Liu F. et al., 2020). Loss and structural changes of tRNAs at interspecific and intraspecific level reflected their rapid evolution in Ulva mtDNAs (Noutahi et al., 2019). The functions of the lost tRNA genes could be replaced by those from the nucleus or other organelles (Adams and Palmer, 2003; Pino et al., 2010).

We found that the tRNA duplication events occurred frequently in the evolution of Ulva mitogenomes. The mitogenomes of U. ohnoi (Uoh), U. lacinulata (Ulc), Ulva sp. (Usp4), and U. rigida (Uri) have two perfect copies of trnT1(ugu), both of which are located between rnl and cob (Figure 3A). Considering that Uoh, Ulc, and Usp4 have a relatively distant relationship with Uri, it seems that the duplication of trnT1(ugu) happened at least twice independently in Ulva mitogenomes. The trnI1-1(gau) located between trnK1(uuu) and trnA1(ugc) was duplicated in U. fenestrata (Ufe), and the trnI1-2(gau) has been translocated at the upstream of trnV2 (Figure 3B). Our previous study found that trnI1(gau) in U. australis (Uau1-3) was translocated from the intergenic region of trnK1-trnA1 to that of cox3-rpl14 (Liu et al., 2017), which was similar to that in Ulva sp. TM708 (Usp2). Based on their phylogenetic relationships, it is more reasonable that the transposition of trnI1(gau) occurred in common ancestor of the U. australis-Ulva sp.-U. fenestrata clade after its duplication, and then the previous trnI1-1(gau) was subsequently lost in mtDNAs of Uau1-3 and Usp2. In the U. lactuca-U. ohnoi-U. lacinulata-Ulva sp. clade, there are two copies of trnM2(cau) situated in different locations (Figure 3C), i.e., the intergenic regions of rns-trnY1 and cox3-trnX1, respectively, indicating that trnM2(cau) experienced the process of duplication and transposition in this clade.

The frequency of some specific tRNAs differed significantly at the interspecific level (Table 2). For example, trnY2(gua) was observed to exist only in the mtDNAs of U. prolifera (Upr1-3) and U. linza (Uli), and trnX2 was present only in mtDNAs of U. lacinulata (Ulc) and Ulva sp. (Usp4). These specific tRNAs are most likely recruited into the mitogenomes through the integration of exogenous DNA fragments (e.g., mitochondrial plasmids) (Handa, 2008). Overall, the evolution of tRNAs in Ulva mitogenomes is similar to that in other eukaryotic lineages (e.g., insects, fungi and higher plants), which is related to different evolutionary processes, including duplication, transposition, remolding, degeneration, loss and recruitment of tRNAs (e.g., Lang, 2014; Zhang et al., 2018; Li et al., 2021).

Among these specific free-standing orfs, a total of 58 orfs were annotated as the full-length or split DNA-dependent RNA polymerase genes (rpos) in these 33 ulvalean mitogenomes. Only three rpos were complete, including orf903 in U. rigida (Uri), orf956 in Ulva sp. TM637 (Usp1) and orf973 in Ulva sp. (Usp4), and the left 55 orfs ranging in size from 98 to 371 aa seem to be the different remnants after the degradation of the intact rpos (Figure 4). Obviously, the rapid mutation accumulation (e.g., insertion and deletion) led to the frameshift and destruction of rpos in these ulvalean mtDNAs. These rpos were detected to be located in nine specific intergenic regions including nad6-trnS1, trnV1-trnS2, trnM2-1-trnY1, trnY1-trnG, rnl-trnT1-1, trnK1-trnI1-1, trnA1-trnR1, cox3-trnX1, and trnE-trnS3, either in a forward or reverse order (Figure 4).

The Ulva mitochondrial Rpos belong to T7-phage-type RNA polymerase which was a group of single-subunit Rpo (ss-Rpo) family. The ss-Rpo family is structurally and evolutionarily distinct from the multi-subunit family of Rpos (including bacterial and eukaryotic sub-families) (Cermakian et al., 1997). Ulva mitogenomes lack the ancestral a-proteobacterial polymerase genes, as was the same as almost all eukaryotic mitogenomes except jacobids (Yin et al., 2010; Peralta-Castro et al., 2020). Transcription of mitochondrial proteins involves nuclear or mitochondrial encoded single subunit T7-phage-type Rpos which probably has replaced the a-proteobacterial polymerase. In this regard, the mitogenomes in Ulva are completely different from chloroplast genomes which still retain the core subunits (rpoA, rpoB, rpoC1, and rpoC2) of the plastid-encoded Rpo derived from their cyanobacterial ancestor (Liu and Melton, 2021).

Based on the blastp searches, the entire T7-phage-type rpo genes closely related to those in Ulva mtDNAs could be found in the nuclear or mitochondrial genomes of eukaryotes (higher plants, brown algae and fungi), as well as genomes of prokaryotes (bacteria and viruses). A total of 45 Rpo proteins from the GenBank database were finally used for phylogenetic analysis. The ML tree inferring from the analysis of core Rpo domains in the 48 rpo genes showed that three Ulva full-length Rpos clustered together, representing a Chlorophyta-specific Rpo lineage (Figure 5). As found in higher plants and fungi (Clark-Walker, 1992; Handa, 2008; Warren et al., 2016), the observed Rpo genes in Ulva mtDNAs most likely come from mitochondrial linear plasmids which carry genes usually coding for an RNA/DNA polymerase or more often for both, as well as other ORFs (Cermakian et al., 1996). These linear plasmids were observed to be sometimes entirely integrated into the mitogenome. The mitogenome of brown alga Pylaiella littoralis contains a putative completely integrated linear plasmid which harbors an entire T7-phage-type Rpo gene (Rousvoal et al., 1998; Oudot-Le Secq et al., 2001).

At the interspecific level, Ulvalean mitogenomes sometimes shared some specific homologous orfs which were located in different positions in these genomes (Supplementary Table 3). For example, among the 17 homologous orfs in group 1, which were from four intergenic regions (Supplementary Table 3), 13 orfs showed sequence similarity to the putative bacterial transposase with alignment scores of approximately 50%, and the left four orfs are more like remnants of their degradation. The orf124 which was only found in mtDNA of Ulva sp. TM708 (Usp2) contained the transferase hexapeptide (six repeats) of putative bacterial origin. The left specific orfs had little sequence similarity to any PCGs in the GenBank database based on the search of blastp.

All sequenced Ulva mitogenomes shared the same set of core genes, but they showed great variations in the content of specific genes including rpo genes, specific orfs (more than 100 codons) and recruited tRNAs (Tables 1, 2). These variations are mainly caused by different acquisitions of foreign DNA fragments from diverse sources including nucleus, plastids, bacteria, viruses, and mitochondrial plasmids (e.g., Wang et al., 2007; Gandini and Sanchez-Puerta, 2017; Liu F. et al., 2020). The integrations of foreign fragments in Ulva mtDNAs lead to the generation of large specific intergenic regions. Nine intergenic regions flanked by core genes vary greatly in size and sequence even at intraspecific levels, indicating that these hot spot regions are undergoing drastic dynamic changes which involve the recent capture of exogenous DNA fragments (Liu and Melton, 2021) and the frequent mutation or loss of intergenic regions which may occur through recombination-excision processes or by slipped-strand mispairing for small regions (Clark-Walker, 1992). Based on these hot spot regions, specific DNA markers could be designed and developed to specifically identify species or populations in the genus Ulva.

Obviously, these sequenced ulvalean mitogenomes show great changes in intron content from four introns in Ulva sp. TM708 (Usp2), U. lactuca (Ula1-3), U. compressa (Uco3 and 4) and U. fenestrata (Ufe) to 14 in Ulva sp. (Usp4) (Table 3). Introns account for 8.77% of mitogenome in Usp2 to 27.80% in U. linza (Uli) in size. The reported ulvophycean mitogenomes contain a certain number of introns (Pombert et al., 2004, 2006; Melton et al., 2015; Turmel et al., 2016), which are not only ribozymes that catalyze their own splicing, but also retroelements that usually harbor an intron-encoded protein (IEP) and can insert themselves into new locations. Most of introns encoded a homing endonuclease or a reverse transcriptase/maturase (RTM) in Ulva mtDNAs, while IEPs in some introns degenerated and even completely lost (Liu et al., 2017; Liu F. et al., 2020).

From a comparative analysis of these intron locations, a total of 29 intron insertion sites were detected at seven mitochondrial genes (atp1, cox1, cox2, nad3, nad5, rnl, and rns) (Table 3). Three types of group I introns including group IB (complete, LHE), group IB (complete, GHE) and group ID, and three of group II introns including group IIA, group IIB and group II (LHE) were detected in these ulvalean mtDNAs. Five intron insertion sites were found for the first time in Ulva mitogenomes, two introns (intron cox1-214 and cox1-900) belonging to group IB and three (intron atp1-1095, cox1-312, and nad5-1057) belonging to group IIA. Intron DNA sequences at the same insertion site were homologous among these ulvalean mtDNAs, and these cognate introns shared the highly conserved RNA secondary structures of ribozyme components. Similar to that in chloroplast genomes, group IB intron was the most prevalent in ulvalean mitogenomes, and was found at ten insertion sites. The LHEs from different families of group IB introns displayed great genetic diversity (Liu and Melton, 2021), but their ribozyme components showed similarity in RNA secondary structures. One group ID intron (intron cox1-709) formerly detected to be present only in the U. prolifera-U. linza-U. flexuosa-Ulva sp. clade (Liu F. et al., 2020), were found to be scattered in mitogenomes of two Ulva lineages (I and II), indicating this group ID intron might be frequently homing or jumping.

In Ulva chloroplast genomes, intron-encoded homing endonucleases from three distinct families (LAGLIDADG, GIY-YIG, and H-N-H) have been found in group I introns (Wang et al., 2021), while all ORF-containing group I introns previously known in Ulva mitogenomes encode a LAGLIDADG homing endonuclease (LHE). In this study, we found that the IEPs from nine group IB (complete) introns were LAGLIDADG homing endonucleases (LHEs), while intron cox1-214 carried a GIY-YIG homing endonuclease (GHE). Group IB (complete,

GHE) intron was observed to be present only in the mtDNAs of two closely related species, U. lacinulata (Ulc) and Ulva sp. (Usp4), indicating that it is most probably the result of an independent insertion event occurring in their common progenitor. In addition, the GHE was also observed to be encoded in group I (derived, A) intron (e.g., intron psbB-489, psbB-772, psbC-882, and psbD-740) in Ulva chloroplast genomes.

Group IIA/IIB introns were present at eight and five insertion sites in these ulvalean mitogenomes, respectively, and both of them usually encoded an RTM (Dai et al., 2003; Liu F. et al., 2020). All of three newly discovered group II intron, i.e., intron atp1-1095, cox1-312, and nad5-1057, belonged to group IIA intron, indicating that group IIA intron is frequently involved in recent invasion or homing in Ulva mitogenomes. Based on phylogenetic analysis of reverse transcriptase (RT) domains in RTMs from Ulva mtDNAs and cpDNAs, RTs can be clearly clustered in two clades, representing group IIA and IIB lineages, respectively (Figure 6). RT family encoded in intron nad5-1057 had a close relationship with that in intron cox1-760, while both of RT families encoded in intron atp1-1095 and cox1-312 showed novel sequence characteristics, respectively, and represented two new intron families. The conserved YXRYADDXXXGXXG catalytic motif in RT segment 5 was shared by RT domains from group IIA introns (Supplementary Figure 1), which was much longer than those (RYADD) from group IIB introns (Supplementary Figure 2; Bonen and Vogel, 2001). All of group IIA introns (100%) contained an intact tripartite RTM gene in these ulvalean mitogenomes, while the RTMs in some group IIB introns have obviously degenerated or been completely lost, e.g., intron nad3-216 and rnl-1963 in U. gigantea (Ugi), and intron rns-780 in U. ohnoi (Uoh). Similar phenomena have been found in P. littoralis (Ikuta et al., 2008), but the mechanism by which the integrity of RTMs in group IIA introns is maintained successfully remains an open question. Group IIB introns with incomplete or missing RTMs should retain splicing competence to ensure that housekeeping genes function properly.

Interestingly, some group II introns in Ulva mitogenomes did not encode an RTM, but instead encoded an LHE. These mitochondrial LHEs in group II introns have close relationships with that in group IB introns (Liu and Melton, 2021). For a long time, it is considered that there is no genetic evolution relationship between group I and group II introns (Kelchner, 2002; Haugen et al., 2005). These findings indicated that there was a certain genetic relationship between group I and group II introns, at least between group IB (LHE) and group II (LHE) introns. Considering the great differences in ribozyme structure and splicing mechanism between group IB (LHE) and group II (LHE) introns (Seetharaman et al., 2006; Stoddard, 2011), the close genetic relationships of their IEPs are likely related to the evolution processes for both of introns. It appears that group I and II introns might recruit or employ the same IEP components (e.g., LHE) in their evolution, but these processes are still unclear.

Single DNA markers (e.g., ITS, rbcL, and 18S rDNA) or combined marker sequences contain limited phylogenetic signals to understand genetic relationships in Ulva species. Trees inferred from phylogenetic analysis based on these datasets have important implications for the phylogeny and classification of Ulva species (e.g., Hayden et al., 2003; Hayden and Waaland, 2004), but it is still restricted by the number of genetic difference signals in these DNA markers, especially in the differentiation of closely related species (e.g., U. linza vs. U. prolifera) or intraspecific relationships (Shimada et al., 2010; Liu et al., 2013). Phylogenomic trees based on mitogenome data and/or plastid genome data as well as organelle genome structure information could provide a more comprehensive understanding of evolutionary systematics and molecular species concepts in this morphologically simple group of macroalgae, due to their rich genetic information at the genome level (Liu F. et al., 2020; Liu and Melton, 2021).

The Ulva mitogenomes sequenced thus far represent a large portion of the genetic diversity at the intragenus level as they were sampled from each of the major clades in this genus. Phylogenomic analysis using maximum likelihood (ML) method based on the mitochondrial nt and aa sequence datasets showed that the 19 Ulva species were divided into two sister clades with strong support values (100%), representing Ulva lineage I and II, respectively (Figure 7). Many mutation sites in mtDNAs are completely consistent within each of these two lineages, but there are significant differences between these two lineages. For example, a three-base insertion occurred at the 5′ end of rps14 in all mtDNAs of Ulva lineage II, not in the mtDNAs of lineage I and P. percursa (Ppe) (Supplementary Figure 3), and a 15-base insertion was detected in rps2 of Ulva lineage II, not in that of lineage I (Supplementary Figure 4). In addition, some genomic features such as gene content, tRNA duplication, distribution of some introns, gene order, and genome rearrangement could partially support the current evolutionary relationships in lineage I and II. The trnY2(gua) was present only in the U. linza-U. prolifera (LP) clade, two copies of trnM2(cau) were present only in the U. lactuca-U. ohnoi-U. lacinulata-Ulva sp. clade, and the GHE-containing group IB introns were found only in the U. lacinulata-Ulva sp. clade.

Relationships of taxa in the present tree are generally congruent with those based on plastid genome dataset (Liu and Melton, 2021; Wang et al., 2021). Minor topological differences between trees are most likely due to the different evolutionary rates between two gene datasets from mitogenomes and plastid genomes. Although two monophyletic clades in Ulva, namely lineage I and II, were well supported by genomic data (Figure 7), few morphological synapomorphies for these two clades were identified in this group of green macroalgae, due to their high degree of phenotypic plasticity caused by environment conditions (Hayden et al., 2003). Both lineage I and II consist of green seaweeds with multiple morphotypes including tubular and blade morphologies. Considering that the genus Ulva contains more than 80 known species worldwide as well as many cryptic species, further phylogenetics studies including more taxa is needed to clarify the infrageneric taxonomy and to understand the evolutionary relationships in Ulva.