PeDinoflagellate strains TGD and MGD established in Sarai et al. (2020) have been maintained in our laboratories and were used in this study. The culture conditions of the two strains were the same as described in Sarai et al. (2020). The algal cultures were observed by light microscopy regularly and the cells in confluent cultures were harvested by centrifugation at 1,720 g for 10 min. The genomic DNA was extracted from the harvested cells by the cetyl trimethyl ammonium bromide (CTAB) method. The cell pellet was dissolved in 500 μL of CTAB extraction buffer [per 100 mL, 2 g CTAB, 10 mL of 1 M Tris–HCl (pH 8.0), 4 mL of 0.5 M ethylenediaminetetraacetic acid (pH 8.0), and 35 mL of 4 M NaCl] at 65^°C for 1 h. After the cells were completely dissolved, 500 μL of Chloroform:Isoamyl Alcohol Solution [480 μL chloroform: 20 μL Isoamyl Alcohol] was added to the tube. After vigorous vortexing, the aqueous solution was saved in a fresh tube by centrifugation at 13,200 rpm for 5 min at 20^°C. Isopropanol precipitation was carried out by adding 500 μL of isopropanol to the aqueous solution, followed by centrifugation at 13,200 rpm for 20 min at 4^°C. The pellet was rinsed with 1 mL of 70% ethanol, and the supernatant was discarded after centrifugation at 13,200 rpm for 5 min at 4^°C. The DNA pellet was briefly dried up and then dissolved into 50 μL of sterile distilled water.

The genomic DNA sample of strain TGD was subjected to Genome-seq analysis using the Illumina Next-seq platform. Approximately 127 million of 150 base paired-end reads were generated (38.1 Gb in total). Initial reads were examined by fastQC to filter the reads containing low-quality bases (under 35) more than 20%. We trimmed the adapter sequence and excluded low-quality bases by FASTX toolkit, yielding approximately 82 million reads for the analyses described below.

We reconstructed a circular plastid genome of strain TGD by the two steps described below. As the de novo assembling of the 82 million reads was computationally intense, approximately 4 million reads were assembled into 116,551 contigs. TBLASTN search was carried out against the resultant contigs using the amino acid sequences of the plastid-encoded proteins of two green algae, P. minor and Chlorella vulgaris, as queries (GenBank accession numbers NC_016733 and AB001684). We retrieved 8 contigs, which were approximately 71 Kb in total length, as the tentative plastid genome of strain TGD. 37,646 reads were selected for the second assembling by mapping the 82 million reads on the tentative plastid genome contigs. We connected the resultant contigs by combining the paired-end information and PCR experiments. Finally, a circular DNA molecule of 71,225 bp was reconstructed as the plastid genome of strain TGD. SPAdes ver.3.7.1 (Prjibelski et al., 2020) and Bowtie2 (Langmead and Salzberg, 2012) were used for de novo assembling and mapping, respectively.

The genomic DNA sample of strain MGD was shipped to a biotech company (Hokkaido System Science Co., Ltd., Hokkaido, Japan) for Genome-seq analysis using the Illumina Hiseq 2000 platform. Approximately 224 million of 100 bp paired-end reads were generated (33.0 Gb in total). After the quality control (see above), approximately 60 million reads were assembled into 379,981 contigs using Ray (Boisvert et al., 2010). We repeated the TBLASTN search described above and identified a single contig of approximately 100 Kb in length as the plastid genome contig. The PCR experiment using a set of the primers, which were designed based on both edges of the 100-Kb contig (5′-GGGGAGAAATTTCAAGATACGG-3′ and 5′-GGGAGGCAAAGGATAAACTAAACG-3′), amplified a single DNA fragment of approximately 2 Kb in length. We failed to determine the complete nucleotide sequence of the amplicon which is largely composed of the repeats of 84 bp (5′-TTA TTTAATGTCACAAAGCCAAATATATAGGCTTTCTATATGTA GAAAGACCAGTTTATTTAATAAAAAGAATAAATTTTAT GT-3′). We conclude that strain MGD has a circular plastid genome of approximately 102 Kb in length, albeit the exact length of the plastid genome remains uncertain.

We annotated the plastid genomes of strains TGD and MGD as follows. The ORFs encoding polypeptides of equal to or more than 100 amino acid residues were identified by MFannot¹. Genes encoding transfer RNAs (tRNAs) were surveyed by tRNAscan-SE² (Chan and Lowe, 2019). Ribosomal RNAs, RNase P RNA, and introns were investigated by RNAweasel³.

The possibility of RNA editing was explored by comparing RNA-seq reads with the plastid genome sequences. The RNA-seq data of strain TGD (GenBank accession number DRR190720), strain MGD (DRR190721), and L. chlorophorum (DRR124369) were downloaded from DDBJ Sequence Read Archive (Kodama et al., 2012). We also downloaded the RNA-seq data of P. minor SAG 1965-3 (ERR2041093) generated by the 1000 Plant Transcriptomes Initiative⁴ (Carpenter et al., 2019; One Thousand Plant Transcriptomes Initiative, 2019). After quality control with FASTP v.0.12.4 (Chen et al., 2018), all reads were mapped with HISAT2 v.2.2.1 (Kim et al., 2019) to the corresponding plastid genomes. Besides the two plastid genomes of strains MGD and TGD determined in this study, we used those of L. chlorophorum and P. minor deposited in the GenBank database under the accession numbers NC_027093.1 and FJ968740.1, respectively. The numbers of RNA-seq short reads subjected to the mapping and those aligned with the corresponding plastid genomes are summarized in Supplementary Table 1.

The incongruities between the RNA-seq reads and plastid genome sequences were detected with mpileup and bcftools commands from bcftools v.1.9 (Li, 2011; Danecek et al., 2021). In this study, we set the two conditions to identify the positions that underwent RNA editing. First, quality-scores of variant calls were greater than 200. Second, the candidate positions are located in ORF/gene coding regions. Only the positions fulfilling both of the two conditions were regarded as post-transcriptional base conversions.

The phylogenetic relationship among green algal plastids and three peDinoflagellate plastids was examined by the maximum-likelihood (ML) phylogenetic analysis of 50 plastid-encoded proteins. The alignment generated by Kamikawa et al. (2015) was modified by adding the plastid-encoded proteins of strains TGD and MGD, as well as three free-living pedinophycean green algae, namely P. tuberculata, Marsupiomonas sp., and strain YPF701 (Lemieux et al., 2014; Jackson et al., 2018). The amino acid sequences were aligned by MAFFT v.7.490 with L-INS-I option (Katoh, 2002), following manual trimming of ambiguously aligned positions. The final “50-protein” alignment comprised 50 plastid-encoded proteins from 37 taxa (8,736 amino acid positions in total). The 50-protein alignment was subjected to the ML analysis with the LG + Γ + F + C60. Non-parametric bootstrap support values were calculated by 100-replicate ML bootstrap analysis with the LG + Γ + F + C60 + PMSF [posterior mean site frequencies; see Wang et al. (2018)] model. The ML and ML bootstrap analyses described above were repeated after excluding two out of the three peDinoflagellates considered in the original 50-protein alignment. We used IQTREE v.2.1.3 (Nguyen et al., 2015) for both ML and ML bootstrap analyses described above. The alignments generated in this study are available as a part of the Supplementary Data.

We also conducted Bayesian phylogenetic analysis with the CAT + GTR model by using Phylobayes v.1.8a (Lartillot and Philippe, 2004, 2006; Lartillot et al., 2007). Two Markov chain Monte Carlo (MCMC) runs ran for 10,000 cycles and the consensus tree with branch lengths and Bayesian posterior probabilities (BPPs) were calculated after the first one-fourth cycles were discarded as “burn-in.” Note that the maxdiff value stayed large (0.371194) but we found that the trees from the two MCMC chains agreed largely with the ML tree.

The 50-protein phylogeny indicated that the branches of the three peDinoflagellates were much longer than those of the four pedinophycean green algae (see section “Results”). We examined the magnitude of “long branch-ness” of the three peDinoflagellates in the individual single-protein alignments. We prepared a 5-taxon tree comprising the four pedinophycean green algae and one of the three peDinoflagellates. The 5-taxon tree was enforced to have the sister relationship between P. minor and P. tuberculata and that between Marsupiomonas sp. and Pedinophyceae sp. YPF-701. The branch lengths of the 5-taxon tree were optimized over each of the 50 single-protein alignments. The branch length optimization was performed by IQTREE v.2.1.3 (Nguyen et al., 2015) with the LG + Γ + F + C60 model.

For each tree, the length of the peDinoflagellate branch was subtracted by the sum of the lengths of the rest of branches to obtain “branch length-ratios.” We split the 50 plastid-encoded proteins into two functional categories, “photosynthetic” and “non-photosynthetic.” The former category contained 30 proteins involved in photosynthesis, while the latter was composed of 17 ribosomal proteins, translation elongation factor Tu, Ycf3, and ClpP. The branch length-ratios of the 30 photosynthetic proteins were compared with those of the 20 non-photosynthetic proteins by Wilcoxon rank-sum test.

We completely sequenced a single, circular DNA molecule of 71,225 bp in length as the plastid genome of the peDinoflagellate strain TGD (TGD plastid genome). The circular genome map is provided as Figure 2 and the general features are summarized in Table 1. The content of guanine plus cytosine (GC content) is 34.8%. 69 ORFs, the genes for small and large subunit rRNA (rns and rnl), and the genes for 27 tRNAs were identified in the genome. 67 out of the 69 ORFs were functionally assignable. Neither BLAST nor Pfam search provided any clue to the function of orf123 or orf156. By mapping the RNA-seq reads on the nucleotide sequence of the plastid genome, the incongruities of the nucleotide identity were detected between the genome and transcripts at 177 positions (0.0327% of the coding region; Table 2) in 40 ORFs and rnl (marked by stars in Figure 2). Thus, we concluded that base conversion editing occurred at the 177 positions in the plastid gene transcripts. The RpoB-coding region was found to be interrupted by a single stop codon in the genome. Likewise, a frameshift hinders the recovery of the continuous PsaA-coding region in the genome (a zigzag line in Figure 2). After referring to the RNA-seq data mapped on the two regions, we regard rpoB as a pseudogene while psaA is a functional gene. Although rpoB seems to be transcribed at a certain level, the stop codon remains in the corresponding transcripts. On the other hand, the single reading frame encoding the entire PsaA was recovered post-transcriptionally by the insertion of two consecutive nucleotides (Supplementary Figure 1). Thus, we conclude that the psaA transcript receives base insertion editing and is functional. The RpoC2-coding region was split into two separate ORFs (rpoC2_N and rpoC2_C; Figure 2) and we detected the distinct transcripts from the two ORFs by analyzing the RNA-seq data. No intron was found in any ORFs identified. Intergenic regions occupy 17.9% of the plastid genome. Three pairs of ORFs, namely (i) psaM and psbK, (ii) rps8 and rpl36, and (iii) petL and petG are fused to each other (colored in purple in Figure 2). We had no evidence for any post-transcriptional processing of the transcripts from the three fused ORFs, except for base conversion editing. Noteworthy, the coding region for PsbK-PsaM fused protein was found to encode the entire amino acid sequence of Ycf12 on a different reading frame (Supplementary Figure 2). Ten pairs of neighboring ORFs are on distinct reading frames but partially overlap each other. We found a single case of partial overlapping between a tRNA gene and an ORF. The ORF/gene overlappings described above are highlighted by red dots in Figure 2 (see also Supplementary Table 2 for the details). IRs were not found. One of the three Ile codons in the standard genetic code, AUA, most likely assigns Met (AUA = Met) in the TGD plastid genome, as reported in that of L. chlorophorum (Matsumoto et al., 2011).

We recovered a single, circular DNA molecule of approximately 102 Kb in length as the plastid genome of the peDinoflagellate strain MGD (MGD plastid genome). See Figure 3 and Table 1 for the circular genome map and general features, respectively. The precise length of the plastid genome remains uncertain, as the nucleotide sequence of the region composed of the 84 bp-repeats could not be determined completely. The GC content is 34.6%. Seventy-one ORFs, two intron-encoded ORFs, rns, rnl, and the genes for 28 tRNAs were annotated in the genome. By mapping RNA-seq reads on the genome, we detected base conversion editing at 18 positions (0.0281% of the coding region; Table 2) in the transcripts of 12 ORFs, rns, and rnl (marked by stars in Figure 3). No sign of base insertion editing was detected. The functions of three out of the 71 ORFs could not be assigned (i.e., orf158, orf155, and orf172). We regard ycf4 as a pseudogene (ψycf4), as the putative N- and C-termini were found to be coded in different reading frames and few transcripts were mapped on this region (Supplementary Figure 3). A single group II intron, which hosts two ORFs (orf107 and orf355), was found in psbB. The non-coding region occupies 34.6% of the plastid genome. Four pairs of ORFs, namely rpoA and rps9, rpl5 and rps8, rps19 and rps3, and rps12 and rps7 are fused to each other (colored in purple in Figure 3). No post-transcriptional processing was observed for the transcripts from the fused ORFs, except for base conversion editing that occurred on the rpoA-rps9 and rpl5-rps8 transcripts. Eight pairs of neighboring ORFs are on distinct reading frames but partially overlap each other. We found a single case of partial overlapping between a tRNA gene and an ORF. The ORF/gene overlappings described above are highlighted by red dots in Figure 3 (see also Supplementary Table 2 for the details). IRs were not found. The genetic code used in the MGD plastid genome appeared to be the same as those in the plastid genomes of L. chlorophorum and strain TGD (see above).

Until the RNA-seq data became available (Matsuo and Inagaki, 2018), we could not test whether RNA editing occurs to the plastid gene transcripts in L. chlorophorum. In this study, the mapping of the RNA-seq reads on the nucleotide sequence of the L. chlorophorum plastid genome revealed base conversion editing at 188 positions (0.327% of the coding region; Table 2) in 28 ORFs, two rRNA genes (rns and rnl), and isoleucine and asparagine tRNAs (trnI and trnN). The genes, of which transcripts receive base conversions, are marked by stars in Supplementary Figure 4. rpoC1 has been regarded as a pseudogene due to multiple in-frame stop codons (Kamikawa et al., 2015). We detected the rpoC1 transcripts as reported in Kamikawa et al. (2015), albeit the in-frame stop codons remained in the transcript reconstructed from the RNA-seq data.

We applied the above-mentioned analysis to P. minor using the plastid genome data and RNA-seq data. Although the amounts of the RNA-seq reads mapped on the plastid genome were compatible among P. minor and the three peDinoflagellates (Supplementary Table 1), no case of RNA editing was detected in the transcripts from the P. minor plastid genome.

The origin(s) of the plastids of peDinoflagellate strains TGD and MGD were examined by both ML and Bayesian phylogenetic analyses of an amino acid alignment comprising 50 plastid-encoded proteins. Regardless of the method for tree reconstruction, the monophyly of TGD, MGD, and L. chlorophorum was reconstructed (an MLBP of 100% and a BPP of 1.0; Figure 4). The clade of the peDinoflagellates was then nested within pedinophycean green algae with a specific affinity to Pedinomonas spp. The above-mentioned relationships received non-parametric ML bootstrap support values (MLBPs) of 100% and BPPs of 1.0 (Figure 4). Furthermore, the sampling of the peDinoflagellates in the ML analyses gave no impact on the affinity between the peDinoflagellates and Pedinomonas spp. (Figure 4).

The branches of the three peDinoflagellates in the 50-protein phylogeny were found to be longer than those of the pedinophycean green algae (Figure 4). The ratios of the branch length for a peDinoflagellate to those of the four pedinophycean green algae (branch length-ratios) were calculated for each protein considered in the alignment, aiming to identify the plastid-encoded proteins that contributed to the long branch-ness observed in the 50-protein phylogeny (Figure 5A). The peDinoflagellate branches were generally long in the vast majority of the single-protein trees, regardless of the peDinoflagellate included (see Supplementary Material). Nevertheless, the branch length-ratios calculated from the non-photosynthetic proteins tend to be larger than those from the photosynthetic proteins (Figure 5B). The Wilcoxon rank-sum test rejected the null hypothesis of no difference between the median values of the branch length-ratios calculated from the two categories at the 1% level (p = 7.50 × 10^–3, 3.49 × 10^–5, and 2.14 × 10^–3 in the comparisons considering L. chlorophorum, strain MGD, and strain TGD, respectively). These results suggest that the overall substitution rates in non-photosynthetic plastid-encoded proteins are higher than those in photosynthetic plastid-encoded proteins in the three peDinoflagellates. The significant difference in substitution rate between photosynthetic and non-photosynthetic plastid-encoded proteins may have been achieved separately in L. chlorophorum, strain MGD, and strain TGD.

Prior to the comparison between the plastid genomes of the three peDinoflagellates and those of pedinophycean green algae, we briefly review the four pedinophycean plastid genomes sequenced completely to date (Table 1). As do phylogenetically diverse algae, the four pedinophycean green algae have the plastid genomes with low GC content (33.4–40.3%) and IRs (Table 1). Although the sizes of the four plastid genomes vary from 91.8 to 126.7 Kb, they appeared to carry similar numbers of functionally assignable ORFs (73–75; henceforth here designated simply as “ORFs”). The plastid genome of Pedinomonas tuberculata (126.7 Kb) appeared to be larger than the rest of the three plastid genomes (91.8–98.3 Kb) compared here, partially due to 7 introns in the P. tuberculata plastid genome but none in others (Table 1).

We revealed 7 of the plastid genome features that were shared among the three peDinoflagellates but were not found or much less obvious in the pedinophycean green algae. (i) The peDinoflagellate plastid genomes commonly lack IRs (Figures 2, 3, and Supplementary Figure 3), although the four pedinophycean plastid genomes bear IRs (Table 1). (ii) We identified more than 10 cases of ORF/gene overlapping/fusion in the peDinoflagellate plastid genomes (Table 1). In contrast, only a few cases of ORF overlapping and no ORF fusion were detected in the four pedinophycean plastid genomes. (iii) The peDinoflagellate plastid genomes appeared to carry fewer numbers of ORFs (i.e., 61–68) than 73–75 ORFs found in those of the four pedinophycean green algae (Table 1). The ORF repertories of the three peDinoflagellate plastid genomes appeared to be similar to each other but the TGD and MGD plastid genomes retain 6–8 ORFs that are absent from the L. chlorophorum plastid genome (Figure 5). (iv) The peDinoflagellate plastid genomes commonly use a deviant genetic code, in which AUA is assigned as methionine, not as isoleucine as the standard genetic code (Table 1). (v) Certain levels of pseudogenization seemingly operated on the peDinoflagellate plastid genomes. A single non-functional ORF was found in the individual peDinoflagellate plastid genomes, albeit none has been reported in the pedinophycean plastid genomes. (vi) Base conversion type RNA editing appeared to occur to the plastid gene transcripts in the three peDinoflagellates (see below). Importantly, we found no evidence for RNA editing on the transcripts from the P. minor plastid genome. Finally, (vii) the significant difference in substitution rate between the plastid-encoded proteins in “photosynthetic” category and those in “non-photosynthetic” category appeared to be shared among the three peDinoflagellates (Figure 5B).

Base conversion editing on plastid gene transcripts was observed in the three peDinoflagellates (the details of the sites received base conversion editing are summarized in a spreadsheet in Supplementary Data). However, the pattern and frequency of RNA editing appeared to vary among the three peDinoflagellate plastids. We observed similar numbers of editing positions in the L. chlorophorum and TGD plastids (188 and 177). In the two peDinoflagellate plastids, more than 90% of the identified editings converted A to G or U to C (Table 2). In the L. chlorophorum plastid, three minor types of base conversion, namely conversion from G to A, that from G to C, and that from C to U, were observed. We observed a single case of the conversion from U to A and 6 cases of the conversion of C to U in the TGD plastid.

The trend of the editing in the MGD plastid appeared to be distinct from those in the L. chlorophorum and TGD plastids described above in two aspects. First, only 18 editing positions were detected in the MGD plastid (Table 2). Second, the conversion from A to G occupied 15 out of the 18 observed editing positions but no case of the conversion from U to C was found (Table 2). The rest of the editings converted A to U, G to A, and C to U.

The majority of the base conversion editing in the transcripts from protein-coding regions in the three peDinoflagellate plastids were found at codon first and second positions. In L. chlorophorum, 86, 60, and five cases of the editing were found at codon first, second, and third positions, respectively, after exclusion of those occurred in tRNAs, rRNAs, and the regions overlapped by two ORFs (Table 3). A similar strong bias of the editing toward codon first/second positions over third positions was observed in the TGD and MGD plastids (Table 3).

A single editing position in the petA transcript appeared to be shared between the L. chlorophorum and TGD plastids. A-to-G conversions occurred at the homologous positions in the L. chlorophorum and TGD petA transcripts (Supplementary Figure 5). These editings introduced the amino acid change from lysine to arginine and that from aspartic acid to glycine in the L. chlorophorum and TGD proteins, respectively. All editings except the case described above were found to occur at unique positions across the three peDinoflagellate plastids.

The plastid phylogeny inferred from plastid small subunit rRNA genes demonstrated the pedinophycean green algal origin of the L. chlorophorum, TGD, and MGD plastids (Sarai et al., 2020). In the current study, we strengthened the pedinophycean green algal origin of the three peDinoflagellate plastids by analyzing the 50-protein alignment (Figure 4). The L. chlorophorum, TGD, and MGD plastids were grouped together and connected to the P. minor and P. tuberculata plastids with full statistical support. At the face value, the 50-protein phylogeny indicates that the three dinoflagellates took up the same pedinophycean green alga as the endosymbionts. However, we had a suspicion of the three peDinoflagellate plastids being grouped together artifactually and misplaced in the tree of green algal plastids, as the proteins encoded in the three peDinoflagellate plastids appeared to evolve much more rapidly than the orthologs in the green algal plastids.

To examine the above-mentioned possibility, we reanalyzed the 50-protein alignment after excluding two out of the three peDinoflagellate plastids alternatively. Significantly, the 50-protein phylogeny constantly grouped the Pedinomonas plastids and one of the three peDinoflagellate plastids together, arguing against their intimate phylogenetic affinity being an artifact in the tree reconstruction. The conclusion deduced from the 50-protein analyses agrees well with the discussion in Sarai et al. (2020)—L. chlorophorum, strains MGD, and strain TGD separately transformed green algal endosymbionts belonging to the genus Pedinomonas or those closely related to Pedinomonas. Nevertheless, there is a large possibility for the endosymbionts in the three peDinoflagellates being closely related but different species belonging to the genus Pedinomonas, as the genuine diversity of Pedinomonas (and their close relatives) is underrepresented by P. minor and P. tuberculata, for which plastid genome data are currently available. The phylogenetic relationship among the three peDinoflagellate plastids should be reexamined after we obtain the plastid genome data from the species that sufficiently represent the diversity of pedinophycean green algae in the natural environment.

In order to address why the pedinophycean endosymbiosis was repeated in distantly related branches in the tree of dinoflagellates, we need to accumulate both environmental and physiological data of pedinophycean green algae and peDinoflagellates. First, the three peDinoflagellates are of marine and most likely took up marine pedinophycean algae as the endosymbionts. Nevertheless, P. minor and P. tuberculata were isolated originally from fresh water and soil environments, respectively. Thus, future studies may explore the diversity and distribution of Pedinomonas and their close relatives in the marine environment. Secondly, the repeated pedinophycean endosymbiosis imply the potential merits for the host cells to bear the pedinophycean-derived plastids in their natural habitats. Thus, it is necessary to examine the key factors for the emergence of peDinoflagellates in the marine environments and the physiological characteristics that differentiate peDinoflagellates from other dinoflagellate species.

The comparative plastid genomics identified the features that are common in the three peDinoflagellate plastids but absent in the pedinophycean plastids (see section “Results”). We here propose that some of the above-mentioned features were achieved in the three peDinoflagellate plastids in parallel.