Putative Pseudo-nitzschia cells were isolated using micropipette and incubated in L1 seawater culture medium at temperature of 18–20°C, with an irradiance of 30 µmol photons m⁻² s⁻¹ and a photoperiod of 12/12 h light/dark. Nine Pseudo-nitzschia strains analyzed in this study were isolated from water samples collected in the Bohai Sea (strains CNS00141, CNS00142, and CNS00159) and the Yellow Sea (strain CNS00130) onboard the research vessel “Beidou” supported by the National Natural Science Foundation of China, Bohai and Yellow Sea Oceanography Expedition (NORC2019-01), the Jiaozhou Bay (strains CNS00133 and CNS00138) onboard the research vessel “Chuangxin” operated by the Jiaozhou Bay Marine Ecosystem Research Station, the East China Sea (strain CNS00150) onboard on the research vessel “Zheyu 2” supported by the Natural Science Foundation of China (NSFC), and the Western Pacific (strains CNS00090 and CNS00097) onboard the research vessel “Kexue” (Figure 1A; Table 1).

All strains isolated and studied in this project were deposited at the KLMEES of IOCAS (Nansheng Chen, chenn@qdio.ac.cn). Morphological features of cells were observed by a ZEISS IMAGER A2 microscope (Carl Zeiss AG, Oberkochen, Germany) equipped with differential interference contrast optics. Species were identified based on their morphological features and the similarity of molecular markers to reference molecular markers of known Pseudo-nitzschia species (Table 2).

DNA samples of nine candidate Pseudo-nitzschia strains were prepared using the modified CTAB method (Doyle and Doyle, 1987), which were used to generate paired-end sequencing libraries of 350 bp in size. Genomic DNAs were sequenced using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) at Novogene (Beijing, China). Raw data of 3.78 – 7.33 Gb were generated for each strain with 150 bp paired-end read lengths. Low-quality reads and adapters were removed from the raw data using Trimmomatic (Bolger et al., 2014). Genome size estimation was conducted using Jellyfish (Marcais and Kingsford, 2011) and GenomeScope (Vurture et al., 2017) with k-mer 17. Genome sizes of nine Pseudo-nitzschia strains were estimated to be ranging from 36.6 M (strain CNS00130 and CNS00133) to 252.8 M (strain CNS00159) (Table S1). 1,000,000 clean reads were randomly selected for each strain for Basic Local Alignment Search Tool (BLAST) (Camacho et al., 2009) search against the National Center for Biotechnology Information (NCBI) NT database for estimating bacterial contamination. Bacterial contamination was negligible (< 0.5%) in the DNA samples of all strains (including CNS00090, CNS00130, CNS00133, CNS00141, CNS00142, CNS00150, and CNS00159), except the strains CNS00097 and CNS00138, which contained 58.78% and 7.24% bacterial contamination, respectively (Table S1). Nuclear genome assemblies of nine strains were assembled using SPAdes (Bankevich et al., 2012) using clean data. Genome sequencing depth was estimated based on the base of clean data and nuclear genome size, considering bacterial contamination (Table S1).

Molecular markers including full-length ITS1-5.8S-ITS2, 18S rDNA, 28S rDNA D1-D3, and rbcL were assembled with SPAdes (Bankevich et al., 2012). Quality assessment was done by aligning paired-end reads against each assembled molecular marker using BWA v0.7.17 (Li and Durbin, 2010), and inspected using IGV v2.8.12 (Robinson et al., 2011). The ITS2 regions were identified according to the method described in a previous study (Ajani et al., 2018), using ITS2 sequences of Pseudo-nitzschia dolorosa strains BP3 and 300 (GenBank accession numbers DQ336151 and DQ336153 respectively) as references. The annotation of Pseudo-nitzschia strains were primarily based on ITS1-5.8S-ITS2 sequences and ITS2 sequences and structures (if necessary). The assembled ITS1-5.8S-ITS2 sequence for each strain was used as a query to search the NCBI NT database using BLAST for the target sequence with the highest bitscore (and percentage ≥ 99%). A Pseudo-nitzschia species was annotated as the species from which the reference ITS1-5.8S-ITS2 sequence was supported by publications (Table 2). This annotation was further validated by examining the ITS2 sequences and structures (if necessary), with the focus on the compensatory base changes (CBCs), which was an important indicator for species identification of Pseudo-nitzschia (Li et al., 2017b; Ajani et al., 2018; Chen et al., 2021). Furthermore, the annotation of a Pseudo-nitzschia strain was also checked by examining other molecular markers (including 18S rDNA, 28S rDNA D1-D3, and rbcL) of this strain for consistency (Table 2).

The maximum likelihood (ML) phylogenetic trees of molecular markers were constructed using MEGA7 with 1000 bootstrap replicates (Kumar et al., 2016). Bootstrap values were shown next to the branches (Felsenstein, 1985). The best-fit models were Hasegawa-Kishino-Yano model (HKY + G + I), Kimura 2-parameter model (K2 + G), and General Time Reversible model (GTR + G) for 18S rDNA, 28S rDNA D1-D3, and rbcL, respectively.

Each cpDNA was de novo assembled using GetOrganelle (Jin et al., 2020), which in turn used SPAdes (Bankevich et al., 2012) for assembly, Bowtie2 (Langmead and Salzberg, 2012) for alignment, and BLAST+ (Camacho et al., 2009) for searches. The paths of the cpDNA were viewed using Bandage version 0.8.1 (Wick et al., 2015). Subsequently, complete cpDNAs were examined by aligning sequencing reads against the cpDNAs using the MEM algorithm of BWA v0.7.17 (Li and Durbin, 2010). Alignments were visualized using IGV v2.8.12 (Robinson et al., 2011). Meanwhile, the sequencing depth of cpDNAs were also calculated. ORF finder (https://www.ncbi.nlm.nih.gov/orffinder) and MFannot (https://megasun.bch.umontreal.ca/RNAweasel/) were used to annotate the cpDNAs. The annotated cpDNA sequences were submitted to GenBank under accession numbers MW853965 (P. hainanensis CNS00090), MW853966 (Pseudo-nitzschia sp. CNS00097), MW715816 (P. delicatissima CNS00130), MW722940 (P. micropora CNS00133), MW722941 (P. americana CNS00138), MW722942 (P. pungens CNS00141), MW722943 (P. multistriata CNS00142), MW722944 (P. cuspidata CNS00150), MW722945 (P. multiseries CNS00159). Gene maps of the annotated Pseudo-nitzschia cpDNAs were drawn using the online program OGDRAW (Greiner et al., 2019).

Because psaE was not found in the cpDNAs of eight Pseudo-nitzschia species constructed in this study, alignment of bas1-ftsH regions of nine Pseudo-nitzschia strains was constructed using MEGA7 (Kumar et al., 2016) to examine the gene losses from the cpDNAs. Because the gene psaE could have been transferred from cpDNAs to their corresponding nuclear genomes via endosymbiotic gene transfer (EGT) (Lommer et al., 2010), to ascertain this possibility, we searched for psaE in the assembled genomes of all eight Pseudo-nitzschia strains whose psaE genes were missing using psaE protein sequence of P. americana (CNS00138) as the query using BLAST+ (Camacho et al., 2009). We also searched for potential psaE in Pseudo-nitzschia genomes and transcriptomes downloaded from NCBI (Table S2). This method was successfully applied previously to identify endosymbiotic gene transfer cases in other diatom species (Liu et al., 2021b). To further verify that the loss of the psaE gene was not due to miss assemblies of the genomes and transcriptomes, we PCR amplified an internal segment of psaE (150 bp) by designing the following PCR primers (F: ACTAATTCATCTAAAGCAA; R: TCGTATTCTTAGAAAAG) based on the alignment of psaE genes of P. americana and other diatom species including Nitzschia ovalis (OK505007), Skeletonema tropicum (MW679507), Thalassiosira nordenskioeldii (MW592698). PCR assays were carried out using genomic DNAs of Nitzschia ovalis, Skeletonema tropicum, Thalassiosira nordenskioeldii, and seven Pseudo-nitzschia strains including CNS00130, CNS00133, CNS00138, CNS00141, CNS00142, CNS00150, and CNS00159 as templates. PCR amplification conditions included an initial denaturation at 94°C for 4 min, followed by 34 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 15 s, elongation at 72°C for 15 s, and a final extension at 72°C for 5 min. To verify the quality of all extracted DNA samples, primers DPrbcL1 (AAGGAGAAATHAATGTCT) and DPrbcL7 (AARCAACCTTGTGTAAGTCTC) (Daugbjerg and Andersen, 1997) were used for the amplification of rbcL gene in all DNA samples. PCR amplification conditions for rbcL gene included an initial denaturation at 94°C for 4 min, followed by 34 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s, elongation at 72°C for 1.5 min, and a final extension at 72°C for 5 min.

A total of 95 PCGs including atpA; atpB; atpD; atpE; atpF; atpG; atpH; atpI; cbbx; ccs1; ccsA; chlI; clpC; dnaB; ftsH; groEL; lysR; petA; petB; petD; petG; petL; petM; petN; psaA; psaB; psaD; psaF; psaJ; psaL; psbB; psbC; psbD; psbE; psbF; psbH; psbI; psbJ; psbK; psbL; psbN; psbT; psbV; psbX; psbY; psbZ; rbcL; rbcS; rpl1, 2, 3, 4, 5, 6, 11, 12, 13, 14, 16, 18, 19, 20, 23, 24, 29, 31, 32, 34, 35; rpoA; rpoB; rpoC1; rpoC2; rps2, 3, 4, 5, 7, 9, 10, 11, 12, 13, 14, 16, 17, 18, 20; secA; secG; secY; sufB; sufC; tatC; ycf3, which were shared among 65 cpDNAs, including 55 previously published Bacillariophyta cpDNAs (Accession number were included in Table S3), nine Pseudo-nitzschia cpDNAs constructed in this study, and Triparma laevis (AP014625) (an Ochrophyta cpDNA used as an outgroup taxa), were used for phylogenetic analysis. The amino acid sequences of each of the 95 PCGs were individually aligned using MAFFT with default parameters (Katoh and Standley, 2013). The regions that were ambiguously aligned in each alignment were deleted using trimAl 1.2rev59 (Capella-Gutierrez et al., 2009) with the parameters gt = 1, and all amino acid sequences were concatenated using Phyutility (Smith and Dunn, 2008). Phylogenetic trees were constructed with IQ-TREE using default parameters (Trifinopoulos et al., 2016). Ultrafast bootstrap analysis with 1000 replicates of the dataset and approximate Bayes test was performed to estimate statistical reliability (Anisimova et al., 2011; Minh et al., 2013). Annotation information noted in the phylogenetic tree was based on Algaebase (Guiry and Guiry, 2021). In addition, a consensus phylogenetic tree was constructed refer to previous studies (Huang et al., 2015; Garrison et al., 2016), ASTRAL (Mirarab et al., 2014) was used for phylogenetic analysis under default settings based on ML trees of 95 shared PCGs constructed by RAxML (Stamatakis, 2014).

The program TOPD-FMTS version 4.6 (Puigbo et al., 2007) was used to compare the similarity of two trees constructed by IQ-TREE and ASTRAL using two different approaches: splits and disagree from the program with 100 repetitions. Phylogenetic trees of cpDNAs (constructed by IQ-TREE), 18S rDNA, 28S rDNA D1-D3, and rbcL were also analyzed by TOPD-FMTS version 4.6 (Puigbo et al., 2007).

Synteny analysis of 10 Pseudo-nitzschia cpDNAs was executed using Mauve v2.3.1 using progressive Mauve with default parameters (Darling et al., 2010). The comparative view of representative cpDNAs was performed using circos-0.69 (Krzywinski et al., 2009). The arrangements of genes in nine Pseudo-nitzschia cpDNAs inverted repeat (IR) region were displayed using OGDRAW (Greiner et al., 2019). IRscope (Amiryousefi et al., 2018) was used for the analyses of IR region contraction and expansion at the junctions of cpDNAs.

Ka/Ks rates were calculated using KaKs_Calculator2.0 (Wang et al., 2010) based on 120 protein-coding gene sequences from 10 Pseudo-nitzschia strains. The nucleotide diversity (Pi) values of Pseudo-nitzschia were evaluated by Perl script. Primer 5 was used to design molecular markers of ycf89 (F: ATGRGTTTARATGAWAA R: KRTCATTTGGAATWGGA) and the ML phylogenetic tree of target sequences of ycf89 was constructed by the method mentioned above.

MCMCTree in PAML (Yang, 1997) was used to perform Bayesian estimation of species divergence times, based on the 109 PCGs shared by Ectocarpus siliculosus (NC_013498), Proboscia sp. (MG755791), Coscinodiscus radiatus (KC509521), Rhizosolenia setigera (MG755793), Thalassiosira pseudonana (EF067921), Chaetoceros muellerii (MW004650), Attheya longicornis (MG755798), Phaeodactylum tricornutum (EF067920), Fragilariopsis kerguelensis (LR812620), and 10 Pseudo-nitzschia cpDNAs. Divergence times were calculated according to methods described previously (Matari and Blair, 2014) and fossil evidence was used to calibrate the molecular clock analyses (Medlin, 2015). Fossil evidence from Late Cretaceous (Turonian) provided a minimum age of 89.8 Mya on the divergence between Rhizosolenia setigera and Coscinodiscus radiatus (5-95% quantiles = 92–118 Mya), fossil evidence from Late Cretaceous (Campanian) pennate diatoms provided a minimum age of 72.1 Mya on the divergence between Thalassiosira and Bacillariophyceae (5-95% quantiles = 74–100 Mya), and Early Jurassic (Toarcian) diatom fossils provided a minimum age of 174 Mya on the divergence between diatoms and Ectocarpus (5-95% quantiles = 176–202 Mya).

Tree topology was constrained to reflect the ML tree, and a GTR substitution model was used. The Markov chain Monte Carlo (MCMC) process of PAML mcmctree was run to sample 1, 000, 000 times, with sample frequency set to 50, after a burn-in of 500, 000 iterations.

Nine putative Pseudo-nitzschia strains (CNS00141, CNS00142, CNS00159, CNS00130, CNS00133, CNS00138, CNS00150, CNS00090, and CNS00097) were first annotated based on their morphological characteristics (Hasle, 1994). Their cells were fusiform or lanceolate in shape and tapered at both ends (Figures 1B–G). In general, each cell contained two plastids symmetrically distributed on either side of the transapical axis. Because morphological features of these strains could not be used to adequately determine their taxonomical status, molecular markers constructed in this study were used to facilitate species identification. The Pseudo-nitzschia strains were first annotated using ITS (ITS1-5.8S-ITS2) sequences (Table 2), and the ITS2 regions of these strains differed by at most one base compared to their reference sequences (Table S4), suggesting that there were no compensatory base changes (CBCs), confirming the ITS-based annotation of the Pseudo-nitzschia strains. ITS-based annotation of the Pseudo-nitzschia strains was supported by all other molecular markers including 18S rDNA, 28S rDNA D1-D3, and rbcL, except the strain CNS00097, which was annotated as P. hallegraeffii based on ITS and ITS2 (Tables S4, S5). Based on 18S rDNA sequence (MZ267115) and 28S rDNA D1-D3 (MZ267146), this strain was annotated as P. simulans based on the high similarities to the reference 18S rDNA sequence (OM807226) and 28S rDNA D1-D3 (MF374776), respectively (Tables 2; Table S5), suggesting that the strain CNS00097 might actually represent an unidentified Pseudo-nitzschia species. Thus, we named it Pseudo-nitzschia sp. CNS00097 (Figure 1H; Figure S1). Phylogenetic analysis of these molecular markers including 18S rDNA (Figure 1H), 28S rDNA and rbcL (Figure S1), which were constructed primarily for strain annotation, supported the above annotation.

Complete cpDNAs were constructed for nine Pseudo-nitzschia strains characterized above. Together with one cpDNA constructed for P. multiseries (KR709240) (Cao et al., 2016), ten cpDNAs corresponding to nine Pseudo-nitzschia species have been constructed altogether (Table 3). Nine newly constructed Pseudo-nitzschia cpDNAs varied substantially, ranging from 116,546 bp (P. americana) to 158,840 bp (P. hainanensis) in length (Figure 2). Interestingly, the lengths of these newly constructed cpDNAs were all substantially longer than that of the recently published cpDNA of P. multiseries, which is 111,539 bp (Cao et al., 2016). Indeed, the length of cpDNA of the P. multiseries strain CNS00159 constructed in this study (123,195 bp) was much longer than the recently published cpDNA of P. multiseries (111,539 bp) (Table 3). Of the nine Pseudo-nitzschia cpDNAs constructed in this study, each had typical four conjoined structures with one long single copy (LSC) (59,316–64,301bp), one short single copy (SSC) (38,030–48,237 bp), and two inverted repeats (IRs) (7,188–23,151 bp). In contrast, a single IR region was present in the recently published P. multiseries cpDNA (Cao et al., 2016) (Table 3), which was the main reason for the shorter length of its cpDNA. LSC, SSC, and two IRs of nine cpDNAs accounted for 40.48–55.03%, 30.37-32.63%, and 12.34-29.15% of the total cpDNA lengths. GC contents of these cpDNAs were rather similar, ranging from 30.69% (P. hainanensis) to 35.67% (Pseudo-nitzschia sp. CNS00097). Coding sequences of these nine cpDNAs showed moderate variations, ranging from 99,196 to 116,087 in length. In contrast, non-coding sequences of these nine cpDNAs varied substantially, ranged from 17,350 to 42,753 in length (Table 3).

The lengths of intergenic regions of all cpDNAs analyzed in this study were short (Figure 3), which was consistent to previous studies (Yu et al., 2018), confirming that cpDNAs of Bacillariophyta are generally compact with short intergenic regions. In cpDNAs of Pseudo-nitzschia species, the average length of intergenic regions in P. hainanensis cpDNA was obviously larger than those in other species (Figure 3; Table S3). In general, besides P. hainanensis cpDNA, cpDNAs of all other Pseudo-nitzschia species had no significant difference in intergenic region length. Nevertheless, cpDNA of P. multistriata and Pseudo-nitzschia sp. CNS00097 had some large values in the intergenic region (Figure S2; Table S3). These large intergenic regions in the cpDNA of P. hainanensis were responsible for its large cpDNA size (158,840 bp).

Among the nine newly constructed cpDNAs, P. hainanensis cpDNA was the largest primarily due to its large IR regions. In contrast, the small P. americana cpDNA was primarily due to its shortened IR regions. The differences of gene numbers between Pseudo-nitzschia species were also caused by the numbers of orf genes. No introns were found in all Pseudo-nitzschia cpDNAs, which was not surprising because introns are generally rare in diatom cpDNAs (Ruck et al., 2014). Four pairs of genes overlapping with each other were found in nine cpDNAs, including rpl4-rpl23 (8 bp), psbC-psbD (53 bp), atpD-atpF (4 bp) and sufC-sufB (1 bp). Moreover, a unique pair of overlapping genes orf238-orf126 (7 bp) was found in P. hainanensis (Table 4).

Two gene loss events were identified, including psaE loss from the cpDNAs of all Pseudo-nitzschia species except for that of P. americana, and rpl36 loss from the cpDNA of P. hainanensis. To confirm the loss of psaE was not due to the misannotation of this gene, we aligned the genomic region of P. americana cpDNA containing psaE and its upstream gene bas1 and downstream gene ftsH against syntenic regions of other eight Pseudo-nitzschia strains. The genomic spaces between bas1 and ftsH in all eight Pseudo-nitzschia strains were much shorter than that of P. americana (Figure S3A) and no similarities were identified between pasE and the genomic sequences between bas1 and ftsH in the eight strains (Figure S3B), supporting the loss of psaE from this region. To further explore the possibility that the gene psaE transferred to the nuclear genomes of these eight Pseudo-nitzschia strains via EGT, P. americana psaE protein sequence was used as a query to search for potential targets in the assembled genomes based on Illumina reads of each strain. The searches did not find any candidate psaE genes. Further searches using other published Pseudo-nitzschia sequencing data, including nuclear genomes of P. multistriata and P. multiseries, assembled transcriptomes of P. delicatissima and P. pungens, also did not find candidate psaE genes (Table S2). To test the possibility that the genome and transcriptome assemblies might miss the regions containing psaE gene, we carried out PCR reactions using primers (as described in Materials and Methods) designed against an internal region of psaE. PCR experiments of rbcL gene demonstrated the quality of all DNA samples (Figure S3D), and experiments of psaE gene showed that PCR product around 150 bp was only present in P. americana of seven Pseudo-nitzschia species (Figure S3C). PCR amplification of this region was also successful for other diatom species including Nitzschia ovalis, Skeletonema tropicum, and Thalassiosira nordenskioeldii (Figure S3C), providing independent evidence that gene loss had occurred in eight Pseudo-nitzschia species. rpl36 was also not found in the cpDNAs or in the nuclear genome assemblies of P. hainanensis (CNS00090).

To explore the evolution relationship of 10 Pseudo-nitzschia strains and other diatom species, the amino acid sequences of 95 shared PCGs of Bacillariophyta and Ochrophyta were used for constructing a concatenated tree using the maximum likelihood method (Figure 3). In addition, we also constructed a coalescent tree (Figure S4). These two phylogenetic trees showed highly consistent topologies (Split Distance: 0.1290) with three disagreement taxa including Astrosyne radiata, Toxarium undulatum, and Cylindrotheca closterium (Figure S4). As expected, most diatoms species were well grouped into three main clades corresponding to three classes of Coscinodiscophyceae, Mediophyceae, and Bacillariophyceae, respectively. However, interestingly, Leptocylindrus was sister to all other diatoms, and Attheya plus Biddulphia were sister to Bacillariophyceae.

For nine Pseudo-nitzschia strains in this study, phylogenetic trees based on cpDNAs and different molecular markers showed some differences (Figure 1H; Figure S1; Figure 3; Table S6). However, these molecular markers were primarily used for our species identification, and the phylogenetic tree of cpDNAs was the focus of this study. Based on the phylogenetic tree of cpDNAs, ten Pseudo-nitzschia strains could be grouped into two clades based on their phylogenetic relationships (Figure 3; Figure S4), with clade 1 containing two cpDNAs of two P. multiseries strains, and cpDNAs of P. pungens, P. multistriata, and P. americana, and clade 2 containing cpDNAs of P. hainanensis, P. cuspidata, Pseudo-nizschia sp. CNS00097, P. delicatissima, and P. micropora. A previous study suggested a categorization that can separate Pseudo-nitzschia species into two groups by cell width: (1) seriata group (cell width > 3 μm) and (2) delicatissima group (cell width < 3 μm) (Hasle and Syvertsen, 1997). Based on statistics on the cell size of different Pseudo-nitzschia species (Lelong et al., 2012), species in clade 2 (including P. hainanensis, P. cuspidata, P. delicatissima, and P. micropora) were also known to belong to the delicatissima group (cell width < 3 μm). In contrast, P. multiseries and P. pungens belonged to the seriata group (cell width > 3 μm). Furthermore, P. multistriata and P. americana, whose cell widths span both groups, belonged to neither group.

Comparative analysis of cpDNAs of 10 Pseudo-nitzschia strains showed that these cpDNAs can be divided into four groups (Figure 4), compared with the two clades revealed by phylogenetic analysis (Figure 3), suggesting that full-length cpDNA synteny provide higher resolution in distinguishing cpDNAs of Pseudo-nitzschia species. The first group containing the cpDNAs of P. pungens (CNS00141), P. multistriata (CNS00142), and P. multiseries (CNS00159, KR709240), the second group containing the cpDNAs of P. delicatissima (CNS00130), Pseudo-nitzschia sp. (CNS00097), P. micropora (CNS00133), and P. cuspidata (CNS00150), and the third and fourth groups each containing a single strain.

Although within groups, cpDNAs showed high collinearity (Figure 4), such as P. pungens and P. multiseries (Figure 5A), cpDNAs of different groups showed substantial genome rearrangements (Figures 4, 5). For example, between the cpDNAs of P. americana and P. multiseries (CNS00159), multiple inversion and translocation events were identified (Figures 4, 5B). Similarly, multiple inversion and translocation events were also identified between the cpDNAs of P. hainanensis and P. multiseries (CNS00159) (Figures 4, 5C), and between the cpDNAs of P. delicatissima and P. multiseries (Figures 4, 5D).

The lengths of the IR regions of cpDNAs of the nine Pseudo-nitzschia species were quite different, ranging from 7,188 bp (P. americana) to 23,151 bp (P. hainanensis). Such large differences in the IR regions may cause differences in the gene content. To test this hypothesis, the arrangements of genes in IR region of nine Pseudo-nitzschia cpDNAs were analyzed (Figure 6A). The topology tree of Figure 6A on the left was constructed based on the phylogenetic tree of cpDNAs. In addition, the junctions JLB (LSC/IRb), JSB (IRb/SSC), JSA (SSC/IRa), and JLA (IRa/LSC) were examined to analyze the contraction and expansion of IR regions of the nine Pseudo-nitzschia species (Figure 6B). Although most IR regions contain nine genes including psaA, psaB, trnP(ugg), ycf89, rns, trnI(gau), trnA(ugc), rnl, and rrn5, many IR regions of these Pseudo-nitzschia species hosts rather different sets of genes (Figure 6; Table S7). The IRa and IRb of the P. americana cpDNA each contained seven genes. Interestingly, these two genes (psaA and psaB) missing from the IRa and IRb regions were located in the LSC region. Thus, compared with the cpDNAs of other eight species constructed in this project that each had two copies of psaA and psaB, the P. americana cpDNA contained a single copy of psaA and psaB. The loss of these two genes in the IRa and IRb regions of the P. americana cpDNA was the main reason for its small size. In contrast, the IRa and IRb regions of the P. hainanensis cpDNA each contained 15 genes (trnG(ucc), psbE, psbF, psbL, psbJ, psaA, psaB, trnP(ugg), ycf89, rns, trnI(gau), trnA(ugc), rnl, rrn5, and psbA) and five orfs (orf119, orf295, orf123, orf166, and orf104) (Figure 6A; Table S7). These 15 genes included all nine genes in other cpDNAs. The addition of six genes and five orfs made the IRa and IRb sizes substantially longer than that of other Pseudo-nitzschia species, which was the main reason for the large size of the cpDNA of P. hainanensis.

In addition to the changes that involves gene content contraction (in the P. americana cpDNA) or expansion (in P. hainanensis cpDNA), many other changes have also been observed in the cpDNAs of other Pseudo-nitzschia species (Figure 6A; Table S7). Two orfs (orf167 and orf125) were found to be added to the IRs of P. pungens cpDNA, while four orfs (orf181, orf191, orf173, and orf174) were found to be added to the IRs of Pseudo-nitzschia sp. CNS00097 cpDNA. Moreover, we have found many cases in which genes were found to overlap with junctions. The dnaK gene was found to overlap with the JSB junctions of the cpDNAs of P. multiseries, P. pungens, P. americana, and P. hainanensis, and the rps16 gene was found to overlap with the JSB junction of the P. cuspidata cpDNA (Figure 6B). Similarly, psbB was found to overlap with the JLB junction of the P. americana cpDNA, ycf4 was found to overlap with the JLB junction of the P. hainanaensis cpDNA, and psbJ was found to overlap with the JLB junctions of the cpDNAs of P. delicatissima and P. micropora (Figure 6B; Table S7).

The set of 120 shared protein-coding genes of 10 Pseudo-nitzschia cpDNAs were used to analyze Ka/Ks (Table S8). For these 120 genes, petM showed the highest average Ka/Ks of 0.2231, petN, psbH, and psbL had the lowest average Ka/Ks of 0.0010. All Ka/Ks values were found to be < 1, indicating that all common protein-coding genes in the cpDNAs had purifying selection.

We further examined sequence variability of 150 genes by computing nucleotide diversity (Pi) shared by 10 Pseudo-nitzschia cpDNAs (Figure S5). Among the 10 Pseudo-nitzschia cpDNAs, the Pi values were from 0.0027 (trnP(ugg) and trnR(acg)) to 0.2221 (petF), and the average value of Pi of 150 genes was 0.0847. There were 11 genes ccs1, clpC, dnaB, petF, rpoC2, rps16, secA, secG, secY, thiS, ycf33, ycf41, ycf89, and ycf90 exhibited high Pi values (>0.15). These mutational hotspots can be appropriate loci for developing molecular markers for population genetic studies. Among these 11 genes with high Pi value, the flanking regions of the gene ycf89 were appropriate for designing PCR primers. Phylogenetic trees based on the target sequences suggested that this region could be used as a potential molecular marker. (Figure S6). Primers targeting ycf89, which were described in methods, could be potentially applied to track Pseudo-nitzschia species.

To explore the speciation of Pseudo-nitzschia species, we constructed the time-scale of Pseudo-nitzschia phylogeny (Figure 7). Estimated divergence time result suggested that crown age of Bacillariophyta was dated at approximately 189 Mya. Within the genus Pseudo-nitzschia, all species were divided into two main clades at approximately 41 Mya. P. hainanensis diverged from other Pseudo-nitzschia species at approximately 35 Mya on one of the clades, after which P. cuspidata and Pseudo-nitzschia sp. CNS00097 diverged at about 27 and 19 Mya, and P. delicatissima and P. micropora diverged at about 12 Mya. Within the other clades, the estimated time of divergence between P. americana, P. multistriata, and P. pungens were 30, 21, and 12 Mya, respectively. Thus, comparative analysis of Pseudo-nitzschia cpDNAs suggested that most Pseudo-nitzschia species were generated within the last 40 Mya.