Six candidate Thalassiosira strains (CNS00050, CNS00051, CNS00052, CNS00439, CNS00472, and CNS00561) were isolated from water samples collected in various China coastal regions including the Jiaozhou Bay (Vessel “Chuangxin,” March and June 2019), the Yellow Sea (Vessel “Beidou”, April 2019), the Changjiang Estuary (Vessel “Xiang Yang Hong 18”, August 2020), and the East China Sea (Vessel “Beidou”, April 2019) (Figure 1A and Table 1) using single cell capillary method. These strains were cultured using L1 medium (1 ‰ volume fraction Na₂SiO₃⋅9H₂O was added). The culture temperature was maintained at 19°C, and the illumination intensity was kept from 2000 Lx to 3000 Lx at the photoperiod of 12 h light: 12 h darkness. Morphological features of algal cells were observed using Zeiss microscope (Leica Microsystems). Full-length 18S rDNA sequences of the six strains were obtained by Polymerase Chain Reaction (PCR) amplification using forward primer 28F (5′-CGAATTCAACCTGGTTGATCCTGCCAGT-3′) and reverse primer 42R (5′-CCGGATCCTGATCCTTCT GCAGGTTCACCTAC-3′) (Gu, 2007). PCR amplification was programmed with thermal settings of an initial denaturation at 94°C (4 min), followed by 32 cycles of the following 3 steps: 94°C (1 min), 57°C (1 min 50 s), and 72°C (2 min) and finally, an elongation step of 72°C (10 min). Subsequently, the amplicons were sequenced by the methods of Sanger dideoxy. Full-length 18S rDNA sequences were used as molecular markers for identifying species based on their similarities to 18S rDNA sequences of known Thalassiosira species.

Algal cultures were collected by centrifugation, and the algal muds were stored in liquid nitrogen for subsequent DNA extraction. Total DNA of each sample was extracted using DNAsecure Plant Kit (Tiangen Biotech, Beijing, China). The DNA Library of each sample was prepared using NEB Next^® Ultra™ DNA Library Prep Kit for Illumina (NEB, United States). After qualification, each DNA library was sequenced using Illumina NovaSeq 6000 platform (Illumina, United States) at Novogene (Beijing, China), yielding about 5 Gb sequencing data of paired-end reads with 150 bp in length.

Illumina sequencing results of each strain were assembled using GetOrganelle v1.7.1 (Jin et al., 2020) with ORGANELLE_TYPE: embplant_pt assisted by SPAdes v3.13.2 (Bankevich et al., 2012) and Platanus-allee v2.2.2 (Kajitani et al., 2019). With the cpDNA of T. pseudonana (GenBank accession number: NC_008589) (Oudot-Le Secq et al., 2007) as a reference, we screened the scaffolds of the target cpDNA from the whole genome assembly for cpDNA assembly using BLASTN v2.10.0. The cpDNA assembly of each strain was examined by aligning Illumina reads using the MEM algorithm of BWA v0.7.17 (Li and Durbin, 2010). VarScan v2.3.9 (Koboldt et al., 2009) and IGV v2.8.12 (Robinson et al., 2011) were used to examine mutation sites and to verify the assembly results. Finally, a circular cpDNA was obtained for each strain. Using the same strategy, we also obtained DNA sequences of common molecular markers, including full-length 18S rDNAs, 28S rDNA D1-D2, ITS, 16S rDNA, rbcL, and cox1 for each strain (Supplementary Table 1).

Annotation of cpDNAs was performed using MFannot¹ (transl_table = 11) and annotation results were verified using multiple software including Open Reading Frame Finder (ORF finder)² with genetic code 11 for identifying open reading frames (orfs), tRNAscan-SE 2.0 (Chan and Lowe, 2019) with default setting for identifying tRNAs, MEGA X (Kumar et al., 2018) and BLASTN for identifying rRNAs and introns by sequence comparison with genes from closely related species, and RNAweasel³ for identifying the type of introns. Genome maps were drawn using Organellar Genome DRAW (OGDRAW) (Greiner et al., 2019). Annotations of 55 published diatom cpDNAs were downloaded from NCBI and were verified using the same strategy.

Six newly assembled complete cpDNAs of Thalassiosira species were compared against published cpDNAs of Thalassiosira species including T. weissflogii (NC_025314) (Sabir et al., 2014), T. pseudonana (NC_008589) (Oudot-Le Secq et al., 2007), and T. oceanica (NC_014808) (Lommer et al., 2010). Multiple sequence alignments of complete cpDNAs and synteny analysis were performed using Mauve v2.4.0 (Darling et al., 2010) with default setting. Pairwise gene rearrangements of cpDNAs were visualized using CIRCOS v0.69.9 (Krzywinski et al., 2009). The borders of the large single-copy (LSC), small single-copy (SSC), and inverted repeat (IR) regions were identified and analyzed to show the IR expansions and contractions.

Maximum likelihood (ML) phylogenetic tree was constructed using concatenated amino acid sequences of 95 PCGs (Supplementary Table 1) shared by 61 diatom cpDNAs, including six Thalassiosira and 55 other diatoms, with the Ochrophyta species Triparma laevis (AP014625) as the out-group taxon. Sequence alignment, editing and concatenation were completed using MAFFT v7.471 (Katoh and Standley, 2013), trimal v1.2 (Capella-Gutierrez et al., 2009) with default setting, and phyutility (Smith and Dunn, 2008), respectively. Phylogenetic trees were constructed using IQtree v1.6.12 (Trifinopoulos et al., 2016) with data_type: AA, model: edge-linked partition model with 1000 bootstrap alignments.

We extracted 123 PCGs shared by cpDNAs of nine Thalassiosira strains, which were aligned by MAFFT v7.471 (Katoh and Standley, 2013) and analyzed for nucleotide diversity (Pi) using DnaSP v5.10 (Rozas, 2009). Fragment length differences of the above conserved gene blocks in the nine Thalassiosira cpDNAs were identified using Mauve v2.4.0 (Darling et al., 2010) with a sliding window of 600 bp with the step size of 50 bp. To characterize regions with length differences, Primer Premier 6 (Premier Biosoft International, Palo Alto, CA, United States) was used to design PCR amplification primers (the forward primer: 5′- GTTCTRCATTCRGTACATTCTA -3′ and the reverse primer: 5′- ATCWGGTAGAGCAGTTRTR -3′) for amplifying in the conserved region, located on both sides of the regions with different length. PCR amplification was programmed with thermal settings of an initial denaturation at 94°C (4 min), followed by 32 cycles of the following 3 steps: 94°C (1 min), 50°C (1 min 20 s), and 72°C (2 min) and finally, an elongation step of 72°C (10 min).

We used 127 PCGs (Supplementary Table 2) shared by 23 cpDNAs (including nine Thalassiosira species) to estimate the divergence time at the nucleotide level using MCMCTree in PAML (Yang, 2007). Branch lengths, gradient (g), and Hessian (H) were estimated using maximum likelihood estimates (MLE) (dos Reis et al., 2013, p. 10) and GTR+G substitution model (model = 7) with independent rates clock model (clock = 2). We estimated divergence times by Bayesian analysis. After discarding the first 50,000 (10%) trees as burn-in, twice of 500 million generations were set for the MCMC chains, with sample frequency set to 50. Finally, the tree file was visualized with Figtree v1.4.3⁴.

Age calibration of divergence times was based on the fossil record. First, lower Jurassic diatom fossils were used to define a minimum divergence time of 174 MYA between diatoms and Ectocarpus (5–95% quantiles = 176–202 MYA) (Matari and Blair, 2014). The second calibration point, which corresponded to the prior divergence time of Rhizosolenia, was set at 91.5 MYA (5–95% quantiles = 90–93 MYA) (Sinninghe Damsté et al., 2004), according to the chemical fossil record of C₂₅ highly branched isoprenoid (HBI) alkenes from Upper Turonian. The crown node of Fragilariopsis was constrained with a uniform distribution from 5 to 25 MYA, which estimated by the change or absence of directional movement structures (e.g., terminal fissures) in the fossil (Sims et al., 2019).

Six candidate Thalassiosira strains (CNS00050, CNS00051, CNS00052, CNS00472, CNS00439, and CNS00561) were characterized based on their morphological features (Figure 1B). Similarity of their molecular marker sequences to known Thalassiosira species were used for species annotation (Figure 1C and Supplementary Figure 1). For example, the strain CNS00051 was determined to be T. rotula because its cells were 40–45 μm in diameter, which were similar to reported sizes (Hoppenrath et al., 2007; Li et al., 2013), and each cell contained several yellow-brown chloroplasts (Figure 1Bb). Its full-length 18S rDNA sequence (MW205690) clustered closely with those of two T. rotula strains CCMP1085-20 (Whittaker et al., 2012) and CCMP1812 (Alverson et al., 2007) in the 18S rDNA-based phylogenetic tree (Figure 1C) with high percentage identities (Supplementary Table 3). The annotation of the strain CNS00051 as T. rotula was also supported by the high similarities of its 28S rDNA D1–D2, ITS, 16S rDNA, and rbcL to their corresponding reference sequences (Alverson et al., 2007; Kooistra et al., 2008; Supplementary Figure 1).

Similarly, the strains CNS00050, CNS00052, CNS00472, and CNS00439 were annotated as T. profunda, T. nordenskioeldii, T. tenera, and T. weissflogii, respectively (Figure 1C and Supplementary Table 3). The strain CNS00561 could not be annotated with certainty to a specific Thalassiosira species although it had the typical morphology of Thalassiosira species (Figure 1Bf), suggesting that it may represent a new Thalassiosira species, or a known Thalassiosira species whose molecular markers have not be characterized. Thus, we referred this species as Thalassiosira sp. (CNS00561).

We successfully assembled complete cpDNAs for the above six Thalassiosira strains, each representing one Thalassiosira species (Table 2). The cpDNAs of five Thalassiosira species including T. profunda, T. rotula, T. nordenskioeldii, T. tenera, and Thalassiosira sp. (CNS00561) were assembled for the first time, while the cpDNA of a T. weissflogii strain has been previously published (Sabir et al., 2014). In this study, we assembled a second cpDNA of a different strain of T. weissflogii (CNS00439) that was isolated from the Jiaozhou Bay in China (Table 1). Each of the six cpDNAs was mapped as a single circular DNA, with a typical quadripartite structure consisting of one LSC, one SSC and a pair of inverted repeats (IRa and IRb) (Bendich, 2004; Daniell et al., 2016; Figure 2 and Table 2).

While GC contents of the nine cpDNAs, including the six cpDNAs assembled in this study, were quite similar, falling between 29.4% and 30.9%, sizes of these cpDNAs varied substantially, ranging from 127,601 bp in T. weissflogii (CNS00439) to 170,005 bp in Thalassiosira sp. (CNS00561). The larger cpDNA size of Thalassiosira sp. (CNS00561) was primarily due to the larger sizes of the pair of IRs, which accounted for 38.5% of its complete cpDNA size. Twenty-nine genes were found in each of the two IRs in the cpDNA of Thalassiosira sp. (CNS00561), which was substantially more than those found in the IRs of cpDNAs of other Thalassiosira strains. The cpDNA of T. oceanica was the second largest among the nine cpDNAs, with a length of 141,790 bp and a GC content of 30.4%. Its IR region was also relatively long (23,693 bp), accounting for 33.4% of the full length.

Intergenic spacers were small for all nine Thalassiosira cpDNAs described in this project, with the median sizes ranging from 68.5 to 82.0 bp and the average sizes ranging from 100.7 to 296.85 bp, and the largest intergenic spacers being 2676 bp, located between ycf33 and trnK(uuu) in Thalassiosira sp. (CNS00561). No intergenic spacer was found between a pair of genes rpl24-rpl14 in all nine Thalassiosira cpDNAs as expected because it is common in other reported diatoms, except for those in the cpDNA of Nitzschia palea (MH113811), which has a small interval of 1 bp (Crowell et al., 2019). Of all 61 published diatom cpDNAs, the median of intergenic spacers were less than 100 bp and the average value were less than 300 bp (Figure 3). The average length of intergenic spacers in diatom cpDNAs were much smaller than that of green algae, whose average values were roughly between 300 and 1000 bp (Turmel et al., 2017). This is most likely due to the fact that diatoms have larger IR regions (∼20 kbp for combined length of two IR regions), while green algae have smaller IR regions or no IR regions.

Nevertheless, a small number of relatively large intergenic spacers were identified in the cpDNAs of diatoms, including the intergenic spacer between rbcS and orf143 in the cpDNA of Plagiogramma staurophorum (MG755792), which was 3923 bp in size, the intergenic spacer between ycf3 and orf124 in the cpDNA of Halamphora americana (MK045450), which was 3522 bp in size, and the intergenic spacer between orf99 and atpB in the cpDNA of Halamphora calidilacuna (MK045451), which was 3545 bp in size.

Four pairs of overlapping genes were found in all nine Thalassiosira cpDNAs. The gene pairs atpD-atpF, sufC-sufB, rpl4-rpl23, and psbD-psbC overlapped by 4, 1, 8, and 53 bp, respectively. Among them, the gene pair psbD-psbC overlapped in all other reported diatoms, whose overlap length was between 44 bp in Entomoneis sp. and 57 bp in Roundia cardiophora. The gene pair rpl4-rpl23 also appeared in almost all diatom cpDNAs, except Pseudo-nitzschia multiseries (KR709240) and Fragilariopsis kerguelensis (LR812620). Similarly, the gene pairs atpD-atpF and sufC-sufB were widely found in most diatom cpDNAs, such as Rhizosolenia fallax (MG755802) of Coscinodiscophyceae, Odontella sinensis (Z67753) of Mediophyceae, and Nitzschia palea (MH113811) of Bacillariophyceae.

The nine Thalassiosira cpDNAs shared 155 genes, including 123 PCGs, 27 tRNA, three rRNA, and two additional RNA genes (Table 2). The proportions of coding sequences of these nine cpDNAs were rather different, ranging from 68.0% of Thalassiosira sp. (CNS00561) cpDNA to 86.1% of T. weissflogii cpDNA. The 27 tRNAs were sufficient to satisfy all the requirements for in organelle protein synthesis. All three shared rRNA genes, including 5S rDNA, 16S rDNA, and 23S rDNA, were present in IRs. Two other RNA genes, transfer-messenger RNA (ssra) and plastid signal recognition particle RNA (ffs), were found in all of these nine cpDNAs of Thalassiosira. In addition, the cpDNA of T. oceanica not only had these two RNA genes, but also encoded an additional flrn (ffs-like RNA) (Lommer et al., 2010).

A group II intron was identified in atpB of the T. tenera cpDNA. The intron was 2874 bp in size, contributing to the total length of the T. tenera cpDNA. In contrast, no introns were identified in the cpDNAs of all other Thalassiosira species analyzed in this project. This result was not surprising because introns were found in only five other reported diatom cpDNAs including atpB in Seminavis robusta and Proboscia sp., petD in Plagiogramma staurophorum and Halamphora calidilacuna, psaA in Toxarium undulatum, and petB in Halamphora calidilacuna (Table 3), most of which are group II introns in PCGs. A group I intron with 764 bp in size was also found in the rnl gene of Seminavis robusta. In this study, the intron found in the atpB of T. tenera cpDNA contains one orf (orf590) encoding a peptide of 590 amino acids, which was a reverse transcriptase (RT) and an HNH endonuclease with self-splicing activity. Phylogenetic analysis showed that orf590 had high similarity with orf586 within the intron of atpB of the Chlorophyta species Pseudoneochloris marina cpDNA (KY407657) (Supplementary Figure 2; Turmel et al., 2017).

Similar observations were also found in Seminavis robusta of diatom and Ulva of Chlorophyta (Brembu et al., 2014; Cai et al., 2017; Suzuki et al., 2018). orf582 in the intron of atpB gene of diatom Seminavis robusta has a closer phylogenetic relationship with orfs (orf581, orf572, and orf573) in the introns of atpB gene of three Ulva species (Ulva ohnoi (AP018696), Ulva ohnoi (KX579943), and Ulva compressa (MW548841)). Therefore, it is highly possible that this intron was transferred horizontally from the green algae Pseudoneochloris marina or related species of P. marina to T. tenera.

Compared with the cpDNAs of other Thalassiosira species analyzed in this study (e.g., T. weissflogii, T. pseudonana, T. tenera, T. profunda, T. rotula, and T. nordenskioeldii), the cpDNA of T. oceanica missed one PCG petF, and the cpDNA of Thalassiosira sp. (CNS00561) missed four PCGs including petF, psaE, rpl36, and ycf35. Among these genes, petF has been extensively studied for its role in iron uptake (McKay et al., 1997; Roy et al., 2020). The lack of petF in the cpDNA of T. oceanica has been found to be correlated with the transfer of this gene to the host genome (Lommer et al., 2010; Roy et al., 2020) via endosymbiotic gene transfer (EGT) (Timmis et al., 2004). We hypothesized that the lack of petF in the cpDNA of Thalassiosira sp. (CNS00561) was also correlated with the gain of this gene in the nuclear genome (Figure 4). To test this hypothesis, we searched for potential homologs of these four PCGs in the assembled whole genome sequences based on Illumina reads of Thalassiosira sp. (CNS00561), including both nuclear and organelle genomes. A candidate gene orf128 was found in the whole genome assembly of CNS00561 strain. It encoded a peptide with similarity to petF (63.8–76.0%) (Supplementary Table 4), suggesting that the nuclear gene orf128 corresponded to the cpDNA gene petF. Interestingly, a putative petF gene (named PETF) was also found in the nuclear genomes of T. profunda (CNS00050), T. rotula (CNS00051), and T. nordenskioeldii (CNS00052), respectively, suggesting similar EGT events in these Thalassiosira species. The exact locations of these putative PETF genes in the nuclear genomes remained be determined. The peptides encoded by these PETF genes showed high conservation with the peptides encoded by the cpDNA petF genes, suggesting that these nuclear PETF genes were obtained via EGTs (Figure 4B and Supplementary Table 4). All peptides encoded by the PETF genes had 60 aa at the N-termini that showed high similarity to chloroplast transport peptide, which can be found in the LI818 that codes for chlorophyll-binding light-harvesting protein (Armbrust et al., 2004). petF was also missing in the cpDNA of T. oceanica (Lommer et al., 2010).

The loss of three genes psaE, rpl36, and ycf35 from the cpDNA of Thalassiosira sp. (CNS00561) was not correlated with the gain of the corresponding genes in its nuclear genome, suggesting that these genes might have been lost in evolution. Alternatively, these genes could have been missed in the genome assemblies.

The cpDNA of Thalassiosira sp. (CNS00561) encoded three orfs (orf101, orf103, and orf134) that were not found in the cpDNAs of other Thalassiosira species.

To explore the evolutionary positions of these Thalassiosira strains in diatoms, a core set of 95 PCGs (Supplementary Table 1) shared by 61 diatom cpDNAs were used to construct Maximum likelihood (ML) phylogenetic tree (Figure 5). Most of the nodes in the tree had 100% bootstrap (BP) values, indicating strong support. The class Bacillariophyceae was a monophyletic group and sister to Clade 3c of the class Mediophyceae with 100% BP support. The class Coscinodiscophyceae contained two clades, Clade 2a and 2b with Clade 2a containing a single species Proboscia sp.

The class Mediophyceae contained four clades (Clade 1, 3a, 3b, and 3c). All Thalassiosira species (red) were contained in Clade 3a. In addition to Thalassiosira species, Clade 3a included some species from other genera, highlighting the complexity of their taxonomical relationships (Yu et al., 2018). T. pseudonana cpDNA clustered with those of species of the genus Cyclotella, which was consistent to previous studies that proposed renaming T. pseudonana (Alverson et al., 2011). Similarly, T. weissflogii might be renamed as well because it did not cluster with other Thalassiosira species (T. oceanica, T. tenera, T. profunda, T. rotula, and T. nordenskioeldii). Clade 3a also included the genus Skeletonema.

Comparative analysis of the cpDNAs of species in the genera Thalassiosira, Skeletonema, and Cyclotella with the cpDNA of Roundia cardiophora (KJ958483) (Sabir et al., 2014) uncovered high collinearity (Figure 6). In particular, the cpDNAs of T. weissflogii, T. pseudonana, T. rotula, T. nordenskioeldii, and C. pseudostelligera were perfectly collinear with that of R. cardiophora cpDNA (Figures 6, 7A), suggesting high conservation of these cpDNAs. Nevertheless, many intra-genus genome rearrangement events were identified among the cpDNAs of Thalassiosira species (Figure 6), suggesting the value of cpDNAs as super-barcodes (Li et al., 2015). An inversion event was found between the cpDNAs of two Cyclotella strains (Cyclotella sp. WC03_2 and Cyclotella sp. L04_2) and T. rotula (Figures 6, 7B). The gene block rpl20-rpl19 ranging from rpl20 to rpl19 contained 27 genes with a length of ∼16,000 bp (Supplementary Table 5). Between the cpDNAs of S. pseudocostatum and T. rotula, an inversion involving the complete SSC (rps6-rps16) involving 45 genes with a length of ∼26,500 bp was found (Figure 7C and Supplementary Table 5). This was due to the presence of two equimolar isomers in the cpDNAs, which differ from another by the relative orientation of the SSC (Linne and Kowallik, 1992). Another inversion event was identified between the T. rotula cpDNA and the cpDNAs of T. tenera and T. profunda (Figures 6, 7D). This inversion involved a gene block (clpC-dnaK) with 34 genes in partial SSC regions of T. tenera (CNS00472) and T. profunda (CNS00050) (Figure 7D and Supplementary Table 5). This block spanned ∼21,200 bp.

The cpDNA of Thalassiosira sp. (CNS00561) was exceptional with extensive rearrangements when it was compared with cpDNAs of other Thalassiosira species. For example, its cpDNA and T. rotula cpDNA showed extensive rearrangements (Figure 7E). Nevertheless, conserved multi-gene blocks could still be identified. For example, the large gene block named petA-trnW(cca) in LSC consisted of 16 genes and its gene order was highly conserved with cpDNAs of species of three genera (Thalassiosira, Skeletonema, and Cyclotella). Another gene block, atpA-sufB, containing eight genes, also showed conserved gene order. In addition, the gene order of gene block rps10-dnaK within SSC was shared by the three genera (Thalassiosira, Skeletonema, and Cyclotella) (Supplementary Table 5).

Comparative analysis of the quadripartite structures and corresponding boundaries of cpDNAs of species in the three genera (Thalassiosira, Skeletonema, and Cyclotella) revealed some IR regions experienced expansions or contractions (Figure 8).

In cpDNAs of species of the three genera Thalassiosira, Skeletonema, and Cyclotella, the psaA-trnP(ugg) genes were located in the LSC/IRb boundaries. The distances between psaA and the LSC/IRb boundaries ranged from 92 to 142 bp, while the distances between trnP(ugg) and the LSC/IR boundaries were from 69 to 116 bp. And ccsA-rps6 genes were located at the SSC/IRb boundaries. IRb region only extended into the PCG rps6 with 9 bp inside the IRb region of R. cardiophora cpDNA. However, rps6 of the three genera were completely included inside the SSC with 29 to 55 bp to the SSC/IRa boundaries. At the boundaries of LSC/IRa, the distances between ccsA were 229 to 288 bp to the boundaries. However, except for Cyclotella sp. WC03_2 (KJ958481), Cyclotella pseudostelligera (MG755804) and T. profunda, the SSC/IRb boundaries located in the PCG rps16 with a region ranged from 2 to 21 bp. And in all three genera, the psaJ were situated inside the LSC regions, which were 21 to 134 bp from the LSC/IRa boundaries.

While the IR regions of most cpDNAs were generally conserved with minor exceptions, IR regions of the cpDNAs of two species Thalassiosira sp. (CNS00561) and T. oceanica showed substantial cpDNA genome rearrangements relative to other Thalassiosira cpDNAs, which made some genes enter or leave the IR region, resulting in the difference of gene copy number. For example, the psaA and psaB genes were located in the IR region of Thalassiosira sp. (CNS00561), which were often located in the LSC region of other Thalassiosira species, resulted in one more copy of these genes in the cpDNAs of Thalassiosira sp. (CNS00561). For the same reason, the clpC of T. oceanica also had an additional copy (Figure 6). In contrast, the ccsA gene, which was usually located in the IR region, was located in the LSC region of Thalassiosira sp. (CNS00561), so it had one less copy than other Thalassiosira species.

Each of the 123 PCGs shared in the nine cpDNAs of the genus Thalassiosira was individually aligned and calculated for nucleotide diversity (Pi value) (Supplementary Figure 3). Pi ranged from 0.0162 (psbN) to 0.2237 (rpl29). Three genes, thiS, clpC and rpl29, have higher nucleotide diversity, with Pi greater than 0.2. The photosystem II protein gene (psb), except psbW, had lower Pi than other genes. Similarly, as an important light-regulated gene, psb often has low nucleotide diversity and high homology and conservation in the cpDNAs of higher plants (Xiong et al., 2020; Yang et al., 2020; Wen et al., 2021).

The analysis of Pi value of a sequence from multiple strains requires consistent collinearity. However, due to genome rearrangements among cpDNAs of the genus Thalassiosira, only three conserved gene blocks (atpA-sufB, petA-trnW(cca), and rps10-dnaK) were selected for mutation hotspots identification (Table 3). The block atpA-sufB was about 6700 bp in length and contained eight genes, corresponding to the region 122,247–128,970 bp of the cpDNAs of T. profunda (CNS00050), the block petA-trnW(cca) was about 19,300 bp in length and contained 16 genes, corresponding to the region 94,965–114,298 bp of the cpDNAs of T. profunda (CNS00050), and the block rps10-dnaK was about 17,800 bp in length and contained 31 genes, corresponding to the region 30,305–48,361 bp of the cpDNAs of T. profunda (CNS00050). With a sliding window of 600 bp with the step size of 50 bp, variation hotspots with length differences among the cpDNAs of nine Thalassiosira strains were found in the above-mentioned three gene blocks. A region about 1000 bp in block petA-trnW(cca) [location: 99,796–100,867 bp in the cpDNAs of T. profunda (CNS00050)], involving rpl33, rps20, and rpoB was identified to have high levels of variations (Figure 9A, red arrow). Figure 9A showed that in the region indicated by the arrow, cpDNAs of nine Thalassiosira strains have different sequence lengths, with the longest in the T. pseudonana cpDNA and the shortest in the Thalassiosira sp. (CNS00561) cpDNA. DNA sequences of this region from the cpDNAs had different lengths in these nine strains of Thalassiosira, ranging from 835 to 1087 bp (Figure 9B). The phylogenetic relationship revealed by this region of nine Thalassiosira strains was similar to but not identical to the phylogenetic relatationship of these strains revealed using cpDNA PCGs (Figure 5). We named this region Thalassiosira chloroplast 1, which was abbreviated as thcp1. Phylogenetic analysis based on ML method showed that eight Thalassiosira species could be adequately separated using thcp1 as a molecular marker (Figure 9B). This indicates that thcp1 can be used as molecular marker for distinguishing Thalassiosira species with high resolution. The primers were used to amplify DNAs extracted from 18 phytoplankton species including nine Bacillariophyta species [three Mediophyceae species (red), three Coscinodiscophyceae species (green), and three Bacillariophyceae species (yellow)] (Figure 9C), two Haptophyta species (Isochrysis galbana and Phaeocystis globosa), four Dinoflagellata species (Amphidinium carterae, Alexandrium tamarense, Karenia mikimotoi, and Prorocentrum donghaiense), and three Ochrophyta species (Heterosigma akashiwo, Aureococcus anophagefferens, and Chattonella marina). PCR amplifications were successful for all Thalassiosira species with amplicons with different sizes. However, PCR amplifications were also successful for five additional species of the classes Mediophyceae and Coscinodiscophyceae. With the abundance of cpDNAs information of more species, it is possible to improve the molecular markers to make up for the lack of specificity.

The 127 PCGs (Supplementary Table 2) shared in cpDNAs of the 23 species (15 of Mediophyceae, four of Bacillariophyceae, three of Coscinodiscophyceae and one of Phaeophyceae) were subjected to divergence times estimation (Figure 10).

The results suggested that the first event of diversification within the Phylum Bacillariophyta occurred 129 MYA in lower Cretaceous and the class Coscinodiscophyceae split from it preferentially at 111 MYA. The other two classes, Mediophyceae and Bacillariophyceae were estimated to have separated from each other about 103 MYA. And within Mediophyceae, age estimation for the crown node of Thalassiosira was 51 MYA in the early Eocene. The branching of T. weissflogii in the class Mediophyceae was estimated to have occurred 51 MYA and diverged into different strains at three MYA. Furthermore, the divergence time between the T. pseudonana and Cyclotella was estimated at 27 MYA in Oligocene. The crown node of the remaining Thalassiosira was estimated to have occurred 28 MYA and speciation of them concentrated between 17 and 28 MYA.

In this study, six complete cpDNAs of Thalassiosira strains were assembled, substantially expanded the total number of cpDNAs of Thalassiosira species from three to nine. Thus, cpDNAs of nine Thalassiosira strains corresponding to eight Thalassiosira species were compared and analyzed. The size of these cpDNAs varied significantly, from 127,601 to 170,005 bp, apparently due to a combination of gene density and IR expansions and contractions (Ravi et al., 2007). Among this, gene density was usually related to the length of intergenic spacers, the presence or absence of introns and gain and loss of genes (Ravi et al., 2007; Bedoya et al., 2019). In the nine cpDNAs, Thalassiosira sp. (CNS00561) had the longest intergenic spacers (median value was 82 bp, maximum value was 2676 bp and average value was 297 bp) and the longest IR region (32,702 bp). Similarly, T. oceanica and T. tenera also had longer intergenic spacers and IR region. In addition, the cpDNA of T. tenera had a group II intron with a length of 2874 bp in atpB. As a result of the above reasons, these three cpDNAs were incompact and had long copy length of IR region, which eventually lead to their larger genome size.

Gene rearrangement and the expansion and contraction of IR region are important characteristics for the analysis of genome evolution (Ravi et al., 2007). Although the organization of plastid is generally stable and conserved, gene rearrangement occurs frequently. The fragmented and large-scale rearrangement the cpDNAs of Thalassiosira sp. (CNS00561) and T. oceanica could be regarded as unique genomic characteristics of the two species, and species classification and identification by cpDNAs has been gradually accepted by taxonomists (Kumar et al., 2009; Wu et al., 2010; Yang et al., 2013; Li et al., 2015). As super-barcodes (Li et al., 2015), cpDNAs can show the differences between Thalassiosira species in terms of genome size, gene gain or loss and gene order, which can greatly improve the resolution of low classification levels (Parks et al., 2009) and help to describe species diversity more accurately.

The nine complete cpDNAs of Thalassiosira contained a set of 155 shared genes, including 123 PCGs, 27 tRNA, three rRNA, and two other RNA genes. Chloroplast is considered to be derived from free-living cyanobacteria, and cpDNA is much smaller than that of its predecessors, but it still retains many proteins that are not encoded by itself (Martin, 1998). It is widely accepted that endosymbiotic genes are exported to the nuclear genomes (nDNAs) through EGT (Timmis et al., 2004). Therefore, many primitive PCGs of cpDNAs have been or are being transferred to nDNAs. In addition to the reported missing petF in T. oceanica (Lommer et al., 2010), the cpDNA of Thalassiosira sp. (CNS00561) also missed four PCGs, including petF, ycf35, rpl36, and psaE. The loss of these four genes has also been recorded in the cpDNAs of other species. Among them, petF of Chlorophytes + Streptophytes, Euglenophytes and Alveolates was missing from cpDNAs and transferred to nDNAs (Cui et al., 2006), and the same event also occurred in T. oceanica (Lommer et al., 2010). Similarly, the conserved region of petF can also be found in the non-cpDNA of CNS00561 (PETF), and its front-end amino acid sequence was highly similar to the chloroplast transport peptide of T. oceanica PETF, which may provide the ability for PETF to transfer into chloroplast to play its role. Studies have proposed that the low iron in the sea can be the important factors causing the transfer of petF encoding ferredoxin from cpDNAs to nDNAs, which is an environment-driven evolutionary strategy (Strzepek and Harrison, 2004), which also confirms the remarkable plasticity of chloroplast and nuclear genomes (Lommer et al., 2010). In addition, in three of the six newly constructed cpDNAs of Thalassiosira, including CNS00050, CNS00051, and CNS00052, we found PETF in both cpDNA and non-cpDNA. Such functional PCGs were transferred to the nuclear DNAs, and degenerate copies were often present in organelle genomes, which was considered as an intermediate state of EGT (Brennicke et al., 1993). Until the nuclear copy gene can function normally, the chloroplast copy gene may accumulate some deleterious mutations and be deleted (Ravi et al., 2007).

Except Thalassiosira sp. (CNS00561), psaE was missing in the cpDNAs of five other diatoms, namely, Proboscia sp. (MG755791), Fragilariopsis kerguelensis (LR812620), Pseudo-nitzschia multiseries (KR709240), and Rhizosolenia imbricate (KJ958482), and Rhizosolenia fallax (MG755802). All classes of diatoms, including Mediophyceae, Coscinodiscophyceae, and Bacillariophyceae, have species with psaE loss. It’s worth noting that, the cpDNAs of another species of genus Fragilariopsis, Fragilariopsis cylindrus (NC_045244), and another species of genus Rhizosolenia, Rhizosolenia setigera (MG755793), had complete psaE, suggesting that the loss of psaE in Coscinodiscophyceae occurred at the branch node of Guinardia and the two newly evolved Rhizosolenia. And in class Bacillariophyceae, the loss of psaE occurred at the node where the branch including genus Fragilariopsis and Pseudo-nitzschia was generated, but it was obtained in F. cylindrus again. Similarly, rpl36 was also missing in the above-mentioned Proboscia sp. and R. fallax. However, the loss of ycf35 was more extensive, and it was missing in the cpDNAs of eight species including Proboscia sp., P. multiseries and F. kerguelensis in 61 diatoms. Gene deletion or transfer is a survival mode of adaptation to the environment, which may lead to the diversity of adaptive phenotypes (Roy et al., 2020). However, due to the lack of genome-wide information and the limited sampling of research species, the mechanism and evolutionary relationship of gene loss and transfer are still unclear. In addition, gene gain and loss could be used as a shortcut to distinguish between two species (Hebert et al., 2004), because super-barcodes are more effective than traditional barcodes in detecting genetic loss.

In this study, phylogenetic analysis of 61 diatoms was conducted based on 95 shared PCGs of cpDNAs. The results were consistent with 103 shared PCGs in the cpDNAs of 40 diatoms by Yu et al. (2018), and even expanded the analysis of representative species in each class compared with Yu’s analysis (Yu et al., 2018). This study focused on Clade 3a, including the genus Thalassiosira, which was generated roughly around 93 MYA. The results showed that T. weissflogii first differentiated from the root of all the branches of the genus Thalassiosira around 51 MYA, which supported by the conclusion that the genus Thalassiosira, as mentioned by Sims et al. (2019), evolved around the Eocene-Oligocene boundary. Then, T. pseudonana and two strains of Cyclotella were clustered into one branch with 100% BS, which also supports the view proposed by Alverson (Alverson et al., 2011) that T. pseudonana should be classified in the genus Cyclotella and named as Cyclotella nana. After that, the other six Thalassiosira strains were gradually differentiated at late Miocene to Pliocene around 17–28 MYA. During this period global changes such as the establishment of permanent Antarctic ice-sheets (Barron and Baldauf, 1995), along with changes in ambient temperature, increased the diversity of Thalassiosira species to adapt to the new marine environment (Sims et al., 2019). The genus Thalassiosira was separated by Skeletonema and Cyclotella in phylogeny, which was different from the classification of Thalassiosira based on morphology. The evolutionary status of the three genera may be more accurate when more extensive and balanced sampling analysis of each genus are performed, which may lead to a reclassification of some species.