Seven candidate Chaetoceros species studied in this project were isolated from water samples collected during multiple expeditions in Chinese coastal waters, among which CNS00389, CNS00390, and CNS00516 were isolated from the Jiaozhou Bay (July and August, 2019) on the research vehicle “Chuangxin,” CNS00386 and CNS00396 were isolated from the Changjiang Estuary (July, 2019) on the research vehicle “Zheyu 2,” and CNS00394 was isolated from the East China Sea (May, 2019) on the research vehicle “Xiang Yang Hong 18” (Figure 1A and Supplementary Table 1). The Chaetoceros cells were individually isolated using a micropipette, followed by multiple washes before transferring each single cell to 24-well culture dishes for growth and characterization. These Chaetoceros strains were cultured in L1 medium with 1 ‰ volume fraction Na₂SiO₃ with H₂O added (Guillard and Hargreaves, 1994). The culture temperature was set at 19 ± 1°C, and the illumination intensity was from 2000Lx to 3000Lx at the photoperiod of 12 h light-12 h dark.

Identification of the cultured Chaetoceros strains was done according to both microscopic morphological characters and phylogenetic analysis using universal markers, including full-length 18S rDNA, rbcL, and 28S rDNA D1-D3. The morphological features of the Chaetoceros species were observed using ZEISS Axio Imager 2 (ZEISS, Germany). Molecular marker sequences were assembled using Illumina reads with SPAdes (Bankevich et al., 2012) and GetOrganelle (Jin et al., 2020), with publicly available molecular marker sequences of Chaetoceros species as reference sequences. The assembled sequences were validated by the following steps. (1) Reads were aligned to the assembled sequences using BWA (0.7.17) (Li and Durbin, 2009). (2) Alignment results were extracted using SAMtools (1.10) (Li et al., 2009). (3) Resulting alignments were inspected for validation and error correction using IGV (Thorvaldsdottir et al., 2013). Phylogenetic trees based on molecular markers were constructed using MEGAX (Kumar et al., 2018). Phylogenetic relationships were inferred using the Maximum Likelihood (ML) (Tamura et al., 2004). The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) was shown next to the branches (Felsenstein, 1985).

For DNA library preparation for whole genome sequencing, Chaetoceros cells were collected by centrifugation, and algae mud samples were stored in liquid nitrogen for subsequent DNA extraction. Total DNA was extracted for each sample by using DNAsecure Plant Kit (Tiangen Biotech, Beijing, China). The integrity and purity of DNA were examined by 1% agarose gel electrophoresis and DNA concentration was accurately quantified by Qubit 2.0 Flurometer (Life Technologies, CA, United States). DNA libraries were prepared using NEB Next^TM Ultra^® DNA Library Prep Kit for Illumina (NEB, United States). PCR products were purified using AMPure XP system (Beckman Coulter, Beverly, MA, United States), libraries were analyzed for size distribution using NGS3K/Caliper and quantified using real-time PCR (Qubit^®3.0 Flurometer, Invitrogen, United States). Qualified libraries were sequenced using NovaSeq PE150 (Illumina, San Diego, CA, United States) at Novogene (Beijing, China).

We obtained an average of 5.85 Gb of Illumina paired-end clean sequencing data from genomic DNA of seven Chaetoceros strains. An average of 19,511,552 paired-end reads were retrieved from each sample, with a sequence length of 150 bp. Raw reads in FASTQ format were first processed through a series of quality control (QC) procedures to obtain clean reads, according to method described previously (Song et al., 2020). Complete cpDNAs were assembled using the GetOrganelle (Jin et al., 2020) with clean reads. Chloroplast genome sequences were verified using the same method used for verifying molecular markers described above in 2.1. Annotation of cpDNA was made using MFannot¹ using genetic code of Bacterial, Archaeal, and Plant chloroplast. For genes whose lengths were different from expected, whose start and stop codons were non-canonical, or open reading frames (orfs) that did not show similarity to known genes, Open Reading Frame Finder (ORF finder)² was applied to examine and edit gene models. Annotated results were further validated and formatted using NCBI’s Sequin15.10³. cpDNAs in the Genbank format of the cpDNAs were converted into genome maps by using Organellar Genome DRAW (OGDRAW) online software (Greiner et al., 2019).

Phylogenetic tree based on protein-coding genes (PCGs) was constructed using extracting 95 PCGs (Supplementary Table 2) shared by published Bacillariophyta cpDNAs (Supplementary Table 3) including seven Chaetoceros cpDNAs constructed in this study. The amino acid sequences of each of the 95 PCGs from different diatom cpDNAs were individually aligned using MAFFT v7.310 (–auto) (Katoh and Standley, 2013). Regions that were ambiguously aligned in each alignment were deleted and all amino acid sequences were concatenated using PhyloSuite v1.2.2 (Zhang et al., 2020). Phylogenetic tree was constructed using IQ-TREE v1.6.1 with SH-aLRT support (%)/aBayes support/ultrafast bootstrap support (%) (parameters: -st AA -m TEST -bb 1000 -alrt 1000 -abayes) (Trifinopoulos et al., 2016). Triparma laevis (AP014625) in Ochrophyta was included as out-group taxa.

Alignment of Chaetoceros cpDNAs were performed by using Mauve v2.4.0 (Darling et al., 2010) with default parameters. The cpDNAs borders were analyzed to show the IR expansions and contractions using irscope_pack.3.1 (modified from IRscope) (Amiryousefi et al., 2018).

Chaetoceros cpDNAs were aligned using MAFFT v7.310 (Katoh and Standley, 2013). Nucleotide diversity (Pi), which could be used to estimate the degree of nucleotide sequence variations, which could be used as potential molecular markers, was calculated using the software DnaSP v6.12.03 (Rozas et al., 2017) and cpDNA alignment as input. The window size was set to 600 bp and the step size was 50 bp.

Divergence time estimation was performed by 95 PCGs (Supplementary Table 2) shared by the published 55 Bacillariophyta cpDNAs (Supplementary Table 3) and seven Chaetoceros cpDNAs constructed in this study using MCMCTree in PAML v4.8a (Yang, 2007). Branch lengths, gradient (g) and Hessian (H) were estimated using maximum likelihood estimates (MLE) and GTR + G substitution model (model = 7) with independent rates clock model (clock = 1). Three calibration points⁴ were included in this analysis (Supplementary Table 4), including the calibration point between Ectocarpus siliculosus and diatoms [176.0–202.0 Million years ago (Mya)], the calibration point between Rhizosolenia setigera and Skeletonema pseudocostatum (90.5–91.5 Mya), and the calibration point between Pseudo-nitzschia multiseries and Fragilariopsis cylindrus (10.0–35.3 Mya). The phylogenetic tree was displayed using FigTree v1.4.3 and visualized with 95% highest posterior density interval (HPD) for each node.

All seven Chaetoceros species studied in this project formed chains in which cells were separated by apertures, with long setae protruding from each of the four corners of the cells. These Chaetoceros strains all displayed substantial morphological variations (Li et al., 2017; Xu et al., 2020). The strain CNS00047 was annotated as C. muelleri because these cells were rectangular with long setae, with valve diameter varying from 4.5 to 20.0 μm (Figure 1B), similar to previous description of C. muelleri (Reinke, 1984). Phylogenetics analysis of full-length 18S rDNA sequences of these candidate Chaetoceros species and reference sequences of known Chaetoceros species confirmed that the strain CNS00047 was C. muelleri because its 18S rDNA sequence (MW832831) clustered well with that (AY485453) of C. muelleri (Damste et al., 2004) with high percentage identity (PID) 99.43% (Figure 1I and Table 1). This annotation was also supported by phylogenetic analysis using another molecular marker rbcL (Supplementary Figure 1A and Table 1). The strain CNS00386 was annotated as Chaetoceros costatus, which contained a single, lobed plastid, formed straight chains (Figure 1C; Kooistra et al., 2010). The strains CNS00389 was annotated as C. socialis, which was fan-shaped, with one of the four setae longer than the others and the long setae of adjacent cells joining together (Figure 1D; Degerlund et al., 2012; Pelusi et al., 2019). The strain CNS00390 was annotated as C. pseudo-curvisetus, whose chains were curved, with a large aperture between adjacent cells, and the aperture was large in the middle and small on the sides (Figure 1E; Oku and Kamatani, 1990). The strain CNS00394 was annotated as C. tenuissimus, whose cells were very small, being square to rectangular, with setae being narrow, arising from the two poles of the valve at an angle of 45° to its apical axis (Figure 1F; Sar et al., 2002). The strain CNS00396 was annotated as C. laevisporus, whose cells contained multiple plastids, were rectangular in broad girdle view and formed straight chains (Figure 1G; Chen et al., 2019). The strain CNS00516 was annotated as C. curvisetus, whose chains were helical, with a large elliptical aperture between adjacent cells. Each cell contained only a single plastid, and all setae curve toward the convex side of the chain (Figure 1H). The strains CNS00386, CNS00389, CNS00390, CNS00394, CNS00396, and CNS00516 were further confirmed as C. costatus (Gaonkar et al., 2018), C. socialis (Gaonkar et al., 2017), C. pseudo-curvisetus (Gaonkar et al., 2018), C. tenuissimus (Gaonkar et al., 2018), C. laevisporus (Li et al., 2017), and C. curvisetus (Gaonkar et al., 2018), respectively, according to their molecular features (Figure 1, Supplementary Figure 1, and Table 1).

We constructed full-length cpDNAs for seven Chaetoceros species, among which cpDNAs of six Chaetoceros species (C. costatus, C. curvisetus, C. laevisporus, C. pseudo-curvisetus, C. socialis, and C. tenuissimus) were constructed for the first time. Together with two cpDNAs of C. muelleri and C. simplex that have been previously published (Sabir et al., 2014; Li and Deng, 2021), altogether nine cpDNAs representing eight Chaetoceros species (Table 1) were analyzed in this project. The sizes of these nine Chaetoceros cpDNAs were rather similar, ranging from 116,284 bp (C. muelleri; NC_053621) to 119,034 bp (C. laevisporus; MW845779) (Figure 2A, Supplementary Figure 2, and Table 1). The GC contents of these cpDNAs were also similar (30.26–32.10%). These Chaetoceros cpDNAs all formed typical quadripartite structure with two inverted repeats regions (IRa, IRb), a large single copy (LSC) region, and a small single copy (SSC) region (Figure 2 and Supplementary Figure 2). The lengths of LSC regions of these cpDNAs were similar (ranging from 61,902 to 63,586 bp), so were their SSC regions (ranging from 39,367 to 39,785 bp). In contrast, the lengths of IR regions showed larger variations among these cpDNAs, with the shortest being 6995 bp (C. costatus), while the longest being 8039 bp (C. laevisporus) (Table 1).

Annotation of these cpDNAs revealed that the cpDNAs of five species, including C. curvisetus, C. muelleri, C. pseudo-curvisetus, C. socialis, and C. tenuissimus, each contained 131 PCGs. In contrast, the cpDNAs of C. costatus and C. laevisporus contained different numbers of PCGs, and C. costatus and C. laevisporus contained 128 and 133 genes, respectively (Figure 2A, Supplementary Figure 2, and Table 1). All Chaetoceros cpDNAs contained 30 tRNA and six non-coding rRNA genes (rns, rnl, and rrn5 in IRs) (Figure 2A, Supplementary Figure 2, and Table 1). No introns were found in any of the cpDNAs of these seven Chaetoceros species, which was consistent to previous findings that no introns were identified in cpDNAs of C. simplex (NC_025310) and C. muelleri (NC_053621).

A comparison of 46 genes variably present among the diatom cpDNAs between these nine Chaetoceros strains and Thalassiosira pseudonana revealed many instances of gene gains and losses (Figure 2B). Except for acpP, acpP2, and psaI, the presence or absence of genes in Chaetoceros was generally consistent. These events included peroxiredoxin gene (bas1), three genes encoding subunits of protochlorophyllide reductase (chlB/L/N), the large and small subunits of acetolactate synthase (ilvB/H), cytochrome C6 gene (petJ), two putative serine recombinase genes (serC1 and serC2), putative tyrosine recombinase gene (tyrC), florigen genes (tsf), and hypothetical protein ycf91 (Figure 2B).

These cpDNAs were rather compact, with small intergenic regions, and the median lengths of intergenic regions ranged from 63 to 66 bp. The total lengths of intergenic regions of these nine cpDNAs were similar, ranging from 14,425 bp (12.4% of the total cpDNA of C. muelleri; NC_053621) to 16,621 bp (14.0% of the total cpDNA of C. laevisporus).

To explore the evolutionary relationship between Chaetoceros species and other diatom species, we constructed a phylogenetic tree using 95 PCGs (Supplementary Table 2) that were shared by 62 cpDNAs constructed for Bacillariophyta species (including six sp.). The cpDNA of T. laevis (AP014625), which belonged to the class Bolidophyceae in phylum Ochrophyta, was included as an out-group taxon (Figure 3). All Bacillariophyta species were clustered into three major clades in the phylogenetic tree, corresponding to three classes including Mediophyceae, Bacillariophyceae, and Coscinodiscophyceae, respectively (Figure 3). As expected, all nine cpDNAs of the Chaetoceros species clustered into a single clade. In particular, the cpDNA of C. muelleri (MW845774) and that of C. muelleri (NC_053621) clustered together, and the cpDNA of C. tenuissimus (MW845778) clustered closely with that of C. simplex (NC_025310) (Gaonkar et al., 2018; Figure 3). Additionally, cpDNAs of C. curvisetus and C. pseudo-curvisetus clustered closely, which was consistent with previous report (Gaonkar et al., 2018). In contrast, cpDNAs of C. laevisporus formed an independent clade (Li et al., 2017). The Chaetoceros clade clustered closely with Acanthoceras zachariasii (Chaetocerotaceae), which was consistent with previous study (Yu et al., 2018; Li and Deng, 2021).

Comparative analysis of Chaetoceros cpDNAs revealed near perfect synteny among nine cpDNAs of eight Chaetoceros species (Figure 4A). All genes in the nine cpDNAs exhibited nearly identical gene order (Figure 4B), with only four minor differences identified. First, while a single gene acpP (234–246 bp) was found between rpl20 and trnP in the IRb in the cpDNAs of C. curvisetus, C. muelleri (both MW845774 and NC_053621), C. pseudo-curvisetus, C. simplex, C. socialis, and C. tenuissimus, two genes acpP (240 bp) and acpP2 (267 bp) were found in the corresponding region in the cpDNA of C. laevisporus, and no genes were found in the same region in the cpDNA of C. costatus. Second, the same difference was also found in the IRa. Third, while a single protein-associated ncRNA gene ffs (109 bp) was found between the two genes psbX and trnF in the cpDNAs of C. costatus, C. laevisporus, C. muelleri (MW845774), C. simplex, and C. tenuissimus, it was not found in the corresponding region in the cpDNAs of C. curvisetus, C. muelleri (NC_053621), C. pseudo-curvisetus, and C. socialis. Notably, cpDNAs of two C. muelleri strains were different at this site as well, while ffs was found in C. muelleri (MW845774), it was not found in the cpDNA of C. muelleri (NC_053621), suggesting that this site was highly polymorphic. Lastly, the gene psaI was found in all cpDNAs studied in this project, except the C. costatus cpDNA.

Comparative analysis of the nine Chaetoceros cpDNAs revealed that the IR regions were generally similar, but with important differences (Figure 5). In particular, the IRa/LSC and IRb/SSC boundaries were different among these cpDNAs. Except for the C. costatus cpDNA, the distances between rpl20 and the LSC/IRb boundaries ranged from 38 to 956 bp, while the distances between acpP and the LSC/IRb boundaries ranged from 52 to 263 bp. The distances between psbY and the SSC/IRb boundaries ranged from 54 to 176 bp, while the distances between ycf35 and the SSC/IRb boundaries ranged from 0 to 117 bp. Except for the C. simplex cpDNA, the C. socialis cpDNA, and the C. tenuissimus cpDNA, rpl32 was located at the SSC/IRa boundaries. Except for the C. costatus cpDNA and the C. socialis cpDNA, ycf45 was located at the SSC/IRa boundaries (Figure 5). Because of the loss of the acpP gene in the C. costatus cpDNA, the boundary between the LSC/IRb shifted, causing a contraction of both IRa and IRb regions. The IRa and IRb of the C. costatus cpDNA each contained eight genes (trnP-UGG, ycf89, rns, trnI-GAU, trnA-UGC, rnl, rrn5, and psbY), compared to nine genes in the cpDNAs of C. curvisetus, C. muelleri, C. pseudo-curvisetus, C. simplex, C. socialis, and C. tenuissimus. In contrast, the IRa and IRb of the C. laevisporus cpDNA each contained 10 genes.

Although these nine Chaetoceros cpDNAs showed generally high collinearity, local regions of these cpDNA sequences showed substantial variations at the DNA level. To quantify sequence divergence in these Chaetoceros cpDNAs, we calculated and compared the nucleotide diversity (Pi) values of the Chaetoceros cpDNAs with a window size of 600 bp and a step size of 50 bp. Pi values ranged from 0.0031 to 0.3442. In this analysis, 62 windows were found to have high nucleotide diversity, with Pi greater than 0.25. The main regions contained syfB, ropC2, rpoB genes, and acpP-ycf45 (Figure 6A). Based on the total number of single nucleotide variations (SNVs) and gaps in each window and the ability to distinguish different Chaetoceros species, we identified a hotspots region with 354 SNVs and no gaps (position: 19,025–19,624 bp in C. muelleri cpDNA), which contained a gene syfB (Figure 6B).

To evaluate the resolution power of the hotspots region as a molecular marker, we carried out phylogenetic analysis using it and the result showed that the hotspots region could be used as molecular markers to distinguish different Chaetoceros species (Supplementary Figures 3A, 4A).

Furthermore, based on the sliding window analysis, we also showed distribution of variations of the Chaetoceros cpDNAs. We calculated the actual sequence length of each Chaetoceros cpDNA in each window. We identified a region with high presence and absence variations (ranging from 806 to 961 bp), which corresponded to the region spanning 106,895–107,700 bp in the C. muelleri cpDNA. Phylogenetic analysis of this region showed that it represented a mutation hotspot, which could be used as a potential molecular marker to distinguish different Chaetoceros species (Supplementary Figures 3B, 4B).

Divergence time estimation of the Chaetoceros species was achieved by analyzing DNA sequences of 95 PCGs shared by 62 cpDNAs (Figure 7). The branching of the class Coscinodiscophyceae was estimated to have occurred 131 Million years ago (Mya). The two classes Mediophyceae and Bacillariophyceae were estimated to have separated from their common ancestor 101 Mya. Furthermore, divergence time estimation revealed that the common ancestor of Chaetoceros species, which formed a monophyletic clade at approximately 58 Mya, split from A. zachariasii (Chaetocerotaceae) at about 70 Mya. Among the Chaetoceros species, the age estimate for C. laevisporus was 58 Mya. The divergence time between C. costatus and C. socialis was inferred to have occurred at 33 Mya. The branching of C. muelleri was estimated to have occurred 37 Mya and diverged into different strains at 14 Mya. And the divergence time between C. curvisetus and C. pseudo-curvisetus was inferred to have occurred at about 18 Mya. While the divergence time between C. simplex and C. tenuissimus was estimated at 5 Mya. Taken together, the majority of these Chaetoceros species arose within 50 Mya.