{'fraction_non_alphanumeric': 0.032531475436815976, 'fraction_numerical': 0.022931124338261513, 'mean_word_length': 5.370435086493098, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Diatom Biodiversity and Speciation Revealed by Comparative Analysis of Mitochondrial Genomes

Diatoms (Bacillariophyta) constitute one of the most diverse and ecologically significant groups of phytoplankton, comprising 100,000–200,000 species in three classes Bacillariophyceae, Mediophyceae, and Coscinodiscophyceae. However, due to the limited resolution of common molecular markers including 18S rDNA, 28S rDNA, ITS, rbcL, and cox1, diatom biodiversity has not been adequately ascertained. Organelle genomes including mitochondrial genomes (mtDNAs) have been proposed to be “super barcodes” for distinguishing diatom species because of their rich genomic content, and the rapid progress of DNA sequencing technologies that has made it possible to construct mtDNAs with increasing throughout and decreasing cost. Here, we constructed complete mtDNAs of 15 diatom species including five Coscinodiscophyceae species (Guinardia delicatula, Guinardia striata, Stephanopyxis turris, Paralia sulcata, and Actinocyclus sp.), four Mediophyceae species (Hemiaulus sinensis, Odontella aurita var. minima, Lithodesmioides sp., and Helicotheca tamesis), and six Bacillariophyceae species (Nitzschia ovalis, Nitzschia sp., Nitzschia traheaformis, Cylindrotheca closterium, Haslea tsukamotoi, and Pleurosigma sp.) to test the practicality of using mtDNAs as super barcodes. We found that mtDNAs have much higher resolution compared to common molecular markers as expected. Comparative analysis of mtDNAs also suggested that mtDNAs are valuable in evolutionary studies by revealing extensive genome rearrangement events with gene duplications, gene losses, and gains and losses of introns. Synteny analyses of mtDNAs uncovered high conservation among species within an order, but extensive rearrangements including translocations and/or inversions between species of different orders within a class. Duplication of cox1 was discovered for the first time in diatoms in Nitzschia traheaformis and Haslea tsukamotoi. Molecular dating analysis revealed that the three diatom classes split 100 Mya and many diatom species appeared since 50 Mya. In conclusion, more diatom mtDNAs representing different orders will play great dividends to explore biodiversity and speciation of diatoms in different ecological regions.


INTRODUCTION
Diatoms (Bacillariophyta) constitute a major group of phytoplankton as free living organisms in both marine and freshwater environments or as endosymbionts in dinoflagellates and foraminifers (Field et al., 1998;Pierella Karlusich et al., 2020). Diatoms are unicellular photosynthetic microalgae deemed to be of crucial significance in biogeochemical cycles and the functioning of aquatic food webs (Smetacek, 1998;Armbrust, 2009;Falciatore et al., 2020). Diatoms have the most effective RuBisCO enzyme within autotrophs and work as an essential part in the cycling of CO 2 , which is predicted to be comparable to that of all terrestrial rainforests combined in global carbon cycling (Field et al., 1998;Giordano et al., 2005) and constitute a large part of aquatic biomass contributing approximately 20% of the total primary production on Earth (Nelson et al., 1995;Falkowski et al., 1998). Diatoms also have great value in commercial applications. Combination of different diatom species provides better balanced nutrition and improves aquaculture animal growth better such as Chaetoceros, Skeletonema, and Thalassiosira (Spolaore et al., 2006). Odontella aurita is a marine diatom with rich eicosapentaenoic acid (EPA) content and displays beneficial effect in reducing risk factors for high-fat induced metabolic syndrome (Haimeur et al., 2012). The high neutral lipid content and growth rate of Fistulifera solaris are beneficial for biodiesel production (Tanaka et al., 2015). Diatom species can also pose substantial negative impact on environment. Extreme successions of some diatom species under certain situation can cause harmful algal blooms (HABs) with serious negative consequences on human health, as well as economic loss and social disruption (Jin et al., 2008;Dyson and Huppert, 2010;Lelong et al., 2012;Crosman et al., 2019;Moore et al., 2020). For example, Pseudo-nitzschia pungens has been found to cause HABs in coastal regions of New Zealand and the west coastal regions of America (Patrick et al., 2016).
The critical importance of diatoms is to some degree dictated by their biodiversity. Diatoms are divided into three classes Bacillariophyceae, Mediophyceae, and Coscinodiscophyceae (Medlin and Kaczmarska, 2019). Diatom species show high variations in cell sizes, ranging from a few micrometers to a few millimeters, with single cells forming chains of connected cells for many species (Foster et al., 2011). While morphological features have been effectively used as defining criteria for taxonomy of the diatom species, taxonomy for diatom species with small-cell sizes, and few discrete morphological characters under light microscope are challenging (Mann et al., 2021). Many molecular markers including 18S rDNA, 28S rDNA, ITS, rbcL, and cox1 have been successfully applied alone or in combination with microscopy to study diatom species (Evans et al., 2007;Hamsher et al., 2013;Park et al., 2017). While these common molecular markers are normally effective for distinguishing higher taxa, they are often inadequate for distinguishing congener species, let alone intra-species genetic diversity (Wilcox, 1998;Hollingsworth et al., 2011;Song et al., 2020). As such, although between 100,000 and 200,000 diatom species have been estimated to exist (Malviya et al., 2016), currently only about 12,000 species have been identified (Mann and Vanormelingen, 2013). Lack of molecular markers with sufficient resolution impedes in-depth analysis of diatom biodiversity and speciation, and particularly their differential contribution to ecological systems.
With the development of DNA sequencing technologies and bioinformatics analysis algorithms, mitochondrial genomes (mtDNAs) of thousands species have been fully constructed, which helped to uncover valuable knowledge about phylogenetic relationships and evolutionary trajectories (Smith, 2016). The research of diatom mtDNAs began with the construction of the mtDNA (a circular DNA of 43,827 bp) of Thalassiosira pseudonana (Armbrust et al., 2004). Subsequent construction of the Phaeodactylum tricornutum mtDNA revealed a 35,454-bp repeat sequence, providing a first glimpse of rich divergence among diatoms mtDNAs (Secq and Green, 2011). Further analysis of mtDNAs uncovered the existence of group I intron, encoding a LAGLIDADG endonuclease in the mtDNA of the endosymbiont diatom Kryptoperidinium foliaceum (Imanian et al., 2012). Comparative analysis of five mtDNAs (Toxarium undulatum, Psammoneis japonica, Eunotia naegelii, Cylindrotheca closterium, and Nitzschia sp.) with eleven publicly available diatom mtDNAs suggested that diatom mtDNAs shared strongly conserved genes and highly variable intron content caused by intron loss and horizontal transfer (Guillory et al., 2018). The mtDNAs have also been used as "super barcodes" with higher resolution than common molecular marker genes for comparative genomics analysis (Lv et al., 2020). Mitochondrial genomes have been used to track the genetic diversity of the harmful algal species Eucampia zodiacus (Zhang W. et al., 2021). Despite their critical importance, mtDNAs of only 53 diatom species of seventeen orders have been successfully constructed, with many orders underrepresented. In particular, only four mtDNAs have been constructed for species in the class Coscinodiscophyceae. Such limited availability of mtDNAs hinders in-depth research on diatom biodiversity and the speciation of diatoms.
The Jiaozhou Bay, which is connected to the Yellow Sea with a small opening, is an epitome of China's offshore ecosystem (Liu et al., 2004). Previous studies revealed rich diatom composition in the Jiaozhou Bay, which makes it suitable for studying the biodiversity and speciation of diatom . In this study, we constructed complete mtDNAs of 15 diatom species, which was the first step in constructing mtDNAs for all diatom species in the Jiaozhou Bay. All of these 15 strains were isolated from the Jiaozhou Bay, consisting of strains from 10 diatom orders including two orders Bacillariales (four species) and Naviculales (two species) of the class Bacillariophyceae, four orders Eupodiscales (one species), Hemiaulales (one species), Briggerales (one species) and Lithodesmiales (one species) of the class Mediophyceae, and four orders from Paraliales (one species), Stephanopyxales (one species), Rhizosoleniales (two species), and Coscinodiscales (one species) of the class Coscinodiscophyceae. For accurate comparison, we re-annotated 53 published mtDNAs publicly available at GenBank at NCBI and made several annotation corrections. This work enriched diatoms' mtDNA resources for comparative analysis of biodiversity and speciation, and represents the first large-scale construction of mtDNAs for diatom species from the Jiaozhou Bay.

DNA Library Preparation and Whole Genome Sequencing
Total DNA was extracted with DNAsecure Plant Kit (Tiangen Biotech, Beijing, China). The genomic DNA samples were fragmented by sonication (Covaris S220, Covaris, WBN, United States) to a size of 350 bp. DNA fragments were then end polished, A-tailed, and ligated with the fulllength adapters for Illumina sequencing, followed by PCR (MiniAmp Thermal Cycler, ThermoFisher, MA, United States) enrichment using generic adapter P5 and P7 oligos. The DNA libraries were sequenced using NovaSeq PE150 (Illumina, San Diego, CA, United States) and paired-end reads in size of 150 bp were generated.

Morphological and Molecular Identification of Diatom Strains
All diatom strains were observed using microscope (Zeiss Axio Imager Z2) (Figure 2A). Common molecular marker 18S rDNA for each strain were constructed using whole genome sequencing results. The phylogenetic trees of 18S rDNA were constructed using MEGA7 (Kumar et al., 2016). Phylogenetic relationships were inferred using the Neighbor-Joining method (Saitou and Nei, 1987). The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) was shown next to the branches (Felsenstein, 1985). The analyses involved 50 nucleotide sequences in 18S rDNA gene. All positions containing gaps and missing data were eliminated. There were a total of 1,518 positions in the final dataset of 18S rDNA gene.

Divergence Time Estimation
Phylogenetic analysis and molecular dating were analyzed by calculating the codon evolution rate of 27 mitochondrial proteincoding gene nucleotide sequences (atp6, 8, 9; cob; cox1, 2, 3; nad1-7, 4L, 9, 11; rpl2, 6, 14, 16; rps3, 4, 8, 11, 14; and tatC). The Ectocarpus siliculosus was used as outgroup with its fossil age (Supplementary Table 1). The nucleotide sequence of 27 protein-coding genes (PCGs) were aligned using MAFFT with Codon and concatenated in the software PhyloSuite (Zhang et al., 2020). The phylogenetic tree was performed with IQ-TREE. Molecular dating was conducted using PAML package (Yang, 2007). Briefly, rough estimation of substitution rate was conducted using baseml, and estimation of divergence times with the approximate likelihood method was performed using mcmctree. Five calibration points were used in this analysis (Supplementary Table 1). The phylogenetic tree was displayed in the Figtree and visualized with 95% highest posterior density interval (HPD) for each node.

Construction and Comparative Analysis of Mitochondrial Genomes
We successfully constructed the complete mtDNA for all 15 diatom strains, among which 14 mtDNAs representing 14 species were constructed for the first time (Figure 3). The mtDNA of Cylindrotheca closterium, whose mtDNA was previously published, was constructed for a second strain collected from Jiaozhou Bay for comparative analysis. The sizes of these 15 mtDNAs varied substantially, ranging from 34,323 bp (Guinardia delicatula, order Rhizosoleniales, class Coscinodiscophyceae) to 57,953 bp (Nitzschia ovalis, order Bacillariales, class Bacillariophyceae) ( Table 2). The AT contents also varied substantially, ranging from 68.0% (Nitzschia ovalis) to 77.8% (Actinocyclus sp.). Interestingly, of these mtDNAs, genome sizes and AT contents showed a significant negative correlation (Pearson correlation coefficient = −0.59, p < 0.05), suggesting that mtDNAs with large sizes tend to have relatively lower AT contents. Further study on more mtDNAs is needed to ascertain the relationship between AT contents and mtDNA genome sizes.
With the 15 mtDNAs constructed in this study, there are altogether 68 mtDNA constructed for all Bacillariophyta species. Of these 68 mtDNAs, 39 mtDNAs were for species in Bacillariophyceae, 20 mtDNAs were for species in Mediophyceae, while only nine mtDNAs were for species in Coscinodiscophyceae (Table 3). Most mtDNAs contained   35 PCGs and two non-coding rRNA genes (rns and rnl). The mtDNAs PCGs of Nitzschia traheaformis and Haslea tsukamotoi were exceptional because they possessed three and two copies of cox1, respectively (Figure 4). Different copies of cox1 genes had different lengths, suggesting prolonged evolution after the duplication events. The numbers of tRNA genes varied from 20 (Actinocyclus sp.) to 28 (Psammoneis japonica). Most instances of introns were found in cox1 and rnl (Table 3). However, introns were also found in other genes including nad5 and nad11 in the mtDNAs of Haslea tsukamotoi and Pleurosigma sp., in cox2, cox3, cob, nad7, and rns in Halamphora calidilacuna mtDNA, and nad2 in Eunotia naegelii mtDNA. Among these introns, only four group I introns were found, including one found in cox1 in the mtDNAs of Nitzschia ovalis, one found in rns in the mtDNA of Halamphora calidilacuna, and two found in cox1 and rnl in the mtDNA of Endosymbiont of Kryptoperidinium foliaceum (Table 3). Independent gene loss events have been identified in these mtDNAs. While the 5S rRNA gene (rrn5) was missing in the 52 of 68 mtDNAs, most other genes were missing in mtDNAs of only a few species. The gene rps2 was found missing in eight mtDNAs, including four mtDNAs of species in Mediophyceae and three mtDNAs of species in Coscinodiscophyceae (Figure 4). The gene tatA was found missing in two mtDNAs, both of which were species in Coscinodiscophyceae. Other genes, including rps7, rps10, rps12, rps13, rps19, and rpl5, were found to be missing in single instances of these 68 mtDNAs, suggesting that these genes were likely disposable from the mtDNAs. Indeed, similarity searches of genes rps2, rps13, rps19, and tatA that were lost in mtDNAs in their corresponding nuclear genomes did not identify close homologs. In the class Mediophyceae, mtDNAs of Lithodesmioides sp., Odontella regia, and Lithodesmium undulatum each lost a single PCG rps2. In the class Coscinodiscophyceae, Guinardia striata mtDNA lost a single PCG tatA, while Actinocyclus sp. mtDNA lost four PCGs rps2, rps13, rps19, and tatA. In addition, Coscinodiscus grani mtDNA lost two PCGs rps10 and rpl5, mtDNAs of Coscinodiscus wailesi and Melosira undulata each lost a single PCG rps2. In the class Bacillariophyceae, mtDNAs of three species each lost a single PCGs. Synedra acus mtDNA lost rps2, Pseudo-nitzschia multiseries mtDNA lost rps7, while Halamphora calidilacuna mtDNA lost rps12.

Phylogenetic Positions of the Mitochondrial Genomes of the 15 Species
To explore the evolutionary relationships between mtDNAs of 15 species constructed in this project and those of other diatom species reported previously, we constructed a phylogenetic tree using amino acid sequences of PCGs shared by these mtDNAs using maximum likelihood method. As expected, all diatom species generally fell into three clades, corresponding to three distinct classes Bacillariophyceae, Mediophyceae, and Coscinodiscophyceae ( Figure 5). Mitochondrial genomes representing 39 species in the class Bacillariophyceae were classified into nine orders. Of the 15 mtDNAs constructed in this study, six mtDNAs belonged to the orders Bacillariales and Naviculales in this class. In Bacillariales, Cylindrotheca closterium (MW849265 of the strain CNS00381) clustered with another strain of Cylindrotheca closterium (NC_037986) as expected (Guillory et al., 2018). The mtDNAs of Nitzschia ovalis, Nitzschia traheaformis, and Nitzschia sp. did not cluster with mtDNAs of other species in the genus Nitzschia, suggesting that that the genus Nitzschia may actually represent multiple genera. In the order Naviculales, Haslea tsukamotoi was found to group more closely with Navicula ramosissima than with Haslea nusantara suggesting that that the genus Navicula and Haslea had more complex relationships Prasetiya et al., 2019). And the mtDNA of Pleurosigma sp. was found to form an independent clade (Figure 5).
mtDNA was grouped with that of Odontella regia in the order Eupodiscales as expected (Wang et al., 2021). Hemiaulus sinensis mtDNA was clustered with that of Eucampia zodiacus (Zhang M. et al., 2021) in the order Hemiaulales. Lithodesmioides sp. mtDNA, which belonged to the order Lithodesmiales, was clustered more closely with the mtDNA of Helicotheca tamesis from another order Briggerales than with the mtDNA of Lithodesmium undulatum that belongs to the same order Lithodesmiales (Figure 5).
In the class Coscinodiscophyceae, nine mtDNAs (including five mtDNAs constructed in this study) representing five orders have been constructed. Mitochondrial genome of Paralia sulcata, which represented the first species of the order Paraliales, was clustered with the mtDNA of Stephanopyxis turris, which also represented the first species of the order Stephanopyxales. In the order Rhizosoleniales, mtDNAs of Guinardia delicatula and Guinardia striata from the same genus Guinardia formed a clade within the same order that also included Rhizosolenia setigera. Finally, mtDNA of Actinocyclus sp. was clustered with that of Coscinodiscus granii and Coscinodiscus wailesii (Figure 5).
Taken together, the 15 mtDNAs constructed in this study were grouped into 14 orders of diatom species. Notably, the mtDNAs of Helicotheca tamesis, Paralia sulcata, and Stephanopyxis turris, represented the first mtDNAs of the orders Briggerales, Paraliales, and Stephanopyxales, respectively, providing first reference mtDNAs for these orders.
We further explored the sizes of intergenic regions of mtDNAs in the context of the phylogenetic relationships ( Figure 6). As expected, the intergenic regions of diatom mtDNAs were generally small, with only minor exceptions. Coding sequence of the gene nad11 was split into two segments in the mtDNAs of many diatom species. The availability of 15 additional mtDNAs allowed us to systematically examine the evolutionary features of nad11 (Figure 6). Among 39 mtDNAs of class Bacillariophyceae, nad11 showed split model in 33 mtDNAs (Figure 6). The gene nad11 was not split in the mtDNAs of six species including Pleurosigma sp., Fistulifera solaris, Proschkinia sp., Psammoneis japonica, Asterionella Formosa, and Synedra acus. In contrast, the gene nad11 was not split in all mtDNAs of class Mediophyceae except Odontella aurita var.minima. In the class Coscinodiscophyceae, none of the nad11 genes were split in the mtDNAs of all species, suggesting that the split of nad11 was occurred multiple times in evolution.

Synteny Analysis of Full-Length Mitochondrial Genomes
In the class Bacillariophyceae, six mtDNAs constructed in this study exhibited highly diverse gene arrangements ( Figure 7A). Similar situation was also found in the classes Mediophyceae and Coscinodiscophyceae. Four mtDNAs constructed in this study that belonged to the class Mediophyceae and five mtDNAs constructed in this study that belonged to the class Coscinodiscophyceae also had a variety of rearrangements, respectively (Figures 7B,C).
For exploring the degree of conservation of the mtDNAs at the order level, we conducted a series of synteny analysis in six orders that contained mtDNAs constructed in this study (Figure 8). The mtDNAs of the same order showed more conservation with fewer genome rearrangement events than that of different orders. For example, the mtDNAs of two species Lithodesmioides sp. and Lithodesmium undulatum of the order Lithodesmiales showed nearly identical genome organization. Gene-level analysis revealed that gene translocations and/or inversion in the block of cox2-trnE-trnM-3 appeared in Lithodesmioides sp. and Lithodesmium undulatum, Helicotheca tamesis, and Lithodesmioides sp. (Figures 9D,E). The mtDNAs of three species Guinardia delicatula, Guinardia striata, and Rhizosolenia setigera of the order Rhizosoleniales were also highly syntonic with no rearrangements. Gene-level analysis of Guinardia striata and Guinardia delicatula mtDNAs showed that most genes arrangements were conserved except for the translocations of trnM-1, trnH, and nad3 ( Figure 9F).
In other orders, some genome rearrangements were identified. In the order Bacillariales, a major genome arrangement event FIGURE 5 | Maximum likelihood (ML) phylogenetic tree based on concatenated amino acid sequences encoded by 27 mitochondrial protein-coding genes (atp6, 8, 9; cob; cox1, 2, 3; nad1-7, 4L, 9, 11; rpl2, 6, 14, 16; rps3, 4, 8, 11, 14; and tatC). Mitochondrial protein-coding genes of Phytophthora ramorum were used as outgroup taxa. Numbers on the left and right side at the branches represent Bayesian posterior probabilities and bootstrap values. Branch lengths were proportional to the amount of sequence change, which are indicated by the scale bar below the trees.
(i.e., an inversion) was identified between the mtDNA of Nitzschia sp. and the mtDNAs of three other species in this order (Cylindrotheca closterium, Nitzschia ovalis, and Nitzschia traheaformis) (Figure 8). In the order Naviculales, a major genome rearrangement event (i.e., an inversion) was also identified between the mtDNA of Navicula ramosissima and mtDNAs of two other species (Haslea tsukamotoi and Pleurosigma sp.) (Figure 8). Further gene level analysis of Navicula ramosissima and Haslea tsukamotoi mtDNAs revealed multiple translocation events including trnI, nad11a, trnA, trnR, trnN, trnS, trnM-3, nad11b, trnC, nad4, and nad7 and the inversion of a gene block containing rnl, rns, and nad6 ( Figure 9A). In the order Eupodiscales, multiple translocation events were identified between the mtDNAs of two species of a same genus, Odontella aurita var.minima and Odontella regia, breaking the genes in these two genomes into five gene blocks (Figures 8, 9B). Notably, all genes on the mtDNAs of both Odontella aurita var. minima and Odontella regia were located on a single strand (Figure 8). In the order Hemiaulales, a major inversion event was identified between two species, Hemiaulus sinensis and Eucampia zodiacus. The gene arrangements of Hemiaulus sinensis and Eucampia zodiacus mtDNAs were almost consistent with seven gene rearrangements including trnL-1, trnL-2, rns, rnl, trnM-3, trnH, and nad3 ( Figure 9C).

Divergence Time Estimation for the Diatom Species
The availability of mtDNAs allowed us to explore the divergence time between diatom species, the speciation of which originated in Jurassic period (95% HPD: 173-192 Mya) (Lewitus et al., 2018) (Figure 10). The divergence between the classes Coscinodiscophyceae and Mediophyceae occurred 141 million years ago (Mya, 95% HPD: 159-130 Mya), while the divergence between the classes Mediophyceae and Bacillariophyceae occurred in 131 Mya (95% HPD: 146-120 Mya). The class Bacillariophyceae were represented by two major clades, including Clade A containing a single order Bacillariales, and Clade B containing three orders (e.g., Surirellale, Naviculales, and Cymbellales). The divergence between of the Clade A and Clade B occurred 109 Mya (95% HPD: 124-99 Mya). In Clade B, the order Naviculales was divided into two subclades, suggesting that this order should be divided into two independent orders.
The class Mediophyceae contained Clade C that consisted of the orders Thalassiosirales, Cymatosirales, and Eupodiscales, and Clade D that consisted of the orders Hemiaulales, Lithodesmiales, and Briggerales. In Clade C, the divergence between the order Thalassiosirales and the orders Cymatosirales and Eupodiscales occurred 108 Mya (95% HPD: 127-90 Mya), followed by the divergence of the orders Cymatosirales and Eupodiscales at FIGURE 6 | Boxplot of intergenic region length and the structure of nad11gene of 68 mtDNAs. On the left part, red boxes, green boxes, and purple boxes represented the three classes Bacillariophyceae, Mediophyceae, and Coscinodiscophyceae, respectively. The white point and black vertical line of every box were average and median of intergenic region length, respectively. The scale of intergenic region length is in a log10 units. On the right part, nad11a and nad11b were two unit caused by the break of the nad11 gene.

Diatom Biodiversity and Speciation
Although diatoms have been regarded to be enormously speciesrich with 100,000-200,000 species, only 12,000 species have been described based on their morphological features with some molecular support (Mann and Vanormelingen, 2013). As a consequence, the quantitative understanding of diatom biodiversity has been limited. The recent development of remarkable genomics technologies and bioinformatics methods have made it possible to systematically construct mtDNAs with high throughput and low cost. Thus, organelle genomes including mtDNAs have been proposed to be used as "super barcodes" that provide higher resolution for distinguishing. As a demonstration of this idea, we have constructed 15 mtDNAs for diatom species isolated from the Jiaozhou Bay, which is an epitome of China's offshore ecosystem (Liu et al., 2004). This research has demonstrated that systematic construction of diatom mtDNAs is an effective and fruitful approach for exploring diatom biodiversity and speciation.

Gene Loss and Gene Duplication
Ribosomal PCGs of the mtDNAs have been lost more often than other genes that were related to respiration including cytochrome b (cob) and cytochrome oxidase subunit 1 (cox1) (Adamsa and Palmer, 2003). Among 68 diatom mtDNAs, the absence of specific ribosomal PCGs occurred sporadically and independently in only a few mtDNAs instead of widespread loss, which indicated that the function of the missing genes was still needed and may be directly replaced by some preexisting genes just like the gene substitution in yeast and plant mtDNAs ( Karlberg et al., 2000;Adams et al., 2002). The loss of tatA from mtDNAs was first discovered in the diatom species Guinardia striata and Actinocyclus sp., and rrn5 was found to be absent from most diatom mtDNAs (Oudot-Le Secq et al., 2001).
Gene duplication and subsequent divergence played a vital role in evolution by providing genetic material for novel functions (Hoshino et al., 2001;Taylor et al., 2001). In this study, the duplication and subsequent divergence of cox1 gene FIGURE 8 | Synteny relationships using Mauve analysis for six orders of diatoms (Bacillariales, Naviculales, Eupodiscales, Hemiaulales, Lithodesmiales, and Rhizosoleniales) that included mtDNAs constructed in this study. Each colored block indicates a synteny block among different mitochondrial genomes. in the mtDNAs of Nitzschia traheaformis and Haslea tsukamotoi was found for the first time in diatom mtDNAs. Two cox2 genes were previously found in mtDNA of rapeseed (Brassica napus L.) (Handa, 2003). Further study of more mtDNAs will help us to learn the origin of gene duplication and explore the role gene duplication plays in the evolution of diatom mtDNAs, which will in turn facilitate research on diatom biodiversity and speciation.

Phylogenetic Relationships Between Diatoms
Many studies had been conducted to explore the evolutionary relationships among diatom mtDNAs. For instance, the class Coscinodiscophyceae was found to be non-monophyletic according to the comparative analysis of 17 mtDNAs of diatoms (Pogoda et al., 2019), while the class Mediophyceae was found to be non-monophyletic and was a sister group of the class Bacilariophyceae according to the analysis of 35 mtDNAs of diatoms (Wang et al., 2021). In this study, the phylogenetic analysis of 68 mtDNAs covering 20 orders indicated that Mediophyceae and Coscinodiscophyceae were non-monophyletic. In the order Bacillariales, the taxonomy of species-rich Nitzschia was found to be unstable (Trobajo et al., 2013), which may contribute to the complicated genetic evolution relationships in genus Nitzschia. The published study had proposed transferring Haslea tsukamotoi to Navicula according to morphological data and molecular marker analysis (Li et al., 2017). That Haslea tsukamotoi was closer to Navicula ramosissima than Haslea nusantara in our tree also underscored this point. Similar situation was found in the order Thalassiosirales in which Skeletonema species was found to be closer to species of the genus Thalassiosira rather than Thalassiosira pseudonana. T. pseudonana has been proposed to be more appropriately classified as Cyclotella nana on the basis of the phylogenetic analyses of morphological and molecular datasets (Alverson et al., 2011). More mtDNAs from different orders are needed to facilitate understanding the complex evolutionary relationships between pairs of orders including Naviculales and Cymbellales, Briggerales and Lithodesmiales, and Paraliales and Stephanopyxales.
The sizes of intergenic regions usually correlated linearly with the genome size (Watzlowik et al., 2021). The largest differences were due to repeat sequences. The longest intergenic region was found in Phaeodactylum tricornutum and Asterionella formosa which contained 35 and 25 kb-long repeats, respectively (Secq and Green, 2011;Villain et al., 2017). The nad11 split coding region made nad11 gene undergo fission into two separate submits with its own start and stop codons, respectively and corresponds to the iron-sulfur binding (nad11a) and molybdopterin-binding domains (nad11b) (Guillory et al., 2018). In distinct contrast to the balance of intergenic region average length, the presence of nad11 split coding region was obviously different in the three classes of diatom, only existing in the class of Bacillariophyceae and Mediophyceae, which pointed to that the break of nad11 gene was produced during the evolution of diatoms diversity instead of generating randomly. In the class Coscinodiscophyceae, none of the nad11 genes were split in the mtDNAs of all species suggesting that the split model of nad11 was the derived version, while the non-split model of nad11 was the original version.

Synteny Conservation
There were a variety of rearrangements in three classes, which pointed to that the mtDNA structure was diverse between species in the phylum of Bacillariophyta. Other comparative analyses of mtDNAs had also shown a number of gene rearrangements in the class Eustigmatophyceae and class Phaeophyceae of algae, respectively (Sevcikova et al., 2016;Liu et al., 2019b). Compared to mtDNAs of different orders in the phylum of Bacillariophyta, the mtDNAs within same orders were more conserved. Interestingly, the transcriptional direction of genes on mtDNAs showed substantial variations. While genes are found on both strands of most mtDNAs, genes of mtDNAs of the genus Odontella (Eupodiscales) were found to be located on the same strand. The functional implication of such genome arrangement was yet unknown, but may represent a more economic transcription mode, as suggested previously in Ulva pertusa (Liu et al., 2017). Obviously further studies are needed to explore how such genome arrangement was acquired by Odontella species during evolution, and the functional implication it has in the formation of diversity in Odontella. Aside from the structural similarity, mtDNAs of Lithodesmioides sp. was closer to Helicotheca tamesis than Lithodesmium undulatum in the phylogenetic tree, suggesting that the mtDNA sequences evolved faster than their structures within order Lithodesmiales.

The Origin and Speciation of Diatom Diversity
Fossil evidence suggested that diatoms originated in the late Jurassic period and became more common in the mid-Cretaceous (Round et al., 1990;Sims et al., 2006;Lewitus et al., 2018;Nakov et al., 2018). Despite the reliability of fossil records, many diatom deposits have undoubtedly been lost through diagenesis, making it difficult to track evolutionary trajectory of diatoms. Comparative analysis of diatom mtDNAs could provide valuable complementary genetic insights into diatom evolution. In this study, we first used diatom mtDNAs as a "super barcode" to infer the time frame within which the diatoms originated and diversified. Analysis of fossil records revealed that Mediophyceae and Coscinodiscophyceae diatoms evolved before Bacillariophyceae diatoms (Round et al., 1990). In our study, diatom species of the class Coscinodiscophyceae was first established, followed by diatom species of the class Mediophyceae, and diatom species of the class Bacillariophyceae were the last to appear, which was also generally consistent with the molecular phylogenetic dating analysis using nuclearencoded small-subunit ribosomal RNA (Sorhannus, 2007). Our study indicated that all three classes of diatom emerged 100 Mya during the Cretaceous period (Figure 10). During this time frame, diatoms played dominant role in the carbon cycle, when CO 2 levels of atmosphere were about five times higher than they are today and O 2 levels were rising (Falkowski et al., 2005;Armbrust, 2009). The evolution of the habitat of the class Bacillariophyceae had happened in approximately 120 Mya which reflected the propensity of diatoms to occupy and adapt to a variety of habitats, from marine to the freshwater (Falciatore et al., 2020). Our study also suggested that freshwater species of Bacillariophyceae had appeared in 121 Mya (95% HPD: 112-136 Mya) such as Asterionella formosa. Many species began to appear after the end of Cretaceous, especially those of the order Bacillariales and Thalassiosirales. The mass extinction about 65 Mya which caused loss of about 85% of all species except diatoms and provided ecological opportunities for the establishment of new species (Falkowski et al., 2004). Throughout the Cenozoic, the rising diversity of diatoms had been attributed to an expanded bioavailability of silica from increased weathering of silicate rock (Cermeno et al., 2015) and an influx of nutrientrich seawater into the South Atlantic brought by the Antarctic Circumpolar Current (Barker, 2001;Scher and Martin, 2006).

CONCLUSION
The availability of 15 mtDNAs grouped in ten orders of diatom species provided valuable reference sequences for further evolutionary studies. The 15 newly constructed mtDNAs display several new features of mtDNAs including the identification of first duplication events of the cox1 gene. Further study could help understand the role they play in the evolution and diversity formation of diatoms. The three classes of diatom are separate from each other. The evolutionary relationships for diatom species from different orders were complex, with mtDNAs within orders showing higher similarity. Diatoms of the class Coscinodiscophyceae appeared first in evolution. All three classes emerged 100 Mya during the Cretaceous period, the diversity of diatom began to rise obviously after the mass extinction about 65 Mya and many diatom species were generated 50 Mya. With the construction of more mtDNAs of more diatom species, more insight will be gained into the diatom biodiversity and speciation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers can be found below: https://www.

AUTHOR CONTRIBUTIONS
NC conceived of the project. YW, SL, JW, YY, YC, and ZZ organized the strain selection, cultivation, and genome sequencing. YW and QX organized the assembly, annotation, quality control, and analyzed the data with suggestions from others. YW and NC wrote the manuscript with inputs from others. All authors read and approved the final version of the manuscript.