Fresh thalli of S. ilicifolium var. conduplicatum (specimen number: 2015050009, Supplementary Figure 1A) from Qionghai City, Hainan Province, China (19°6′36″ N, 110°34′48″ E), S. graminifolium (specimen number: 2016030038, Supplementary Figure 1B) from Mirs Bay, Guangdong Province, China (22°32′58″ N, 114°34′7″ E), S. phyllocystum (specimen number: 2016050141, Supplementary Figure 1C) from Qionghai Base, Hainan Province, China (19°20′29″ N, 110°40′09″ E), S. muticum (specimen number: 2017020007, Supplementary Figure 1D) from Heishijiao, Dalian, China (38°52′7″ N, 121°33′30″ E), S. feldmannii (specimen number: 2017030072, Supplementary Figure 1E) from Huitou Xiaodonghai, Sanya City, Hainan Province, China (18°12′24″ N, 109°29′47″ E), S. mcclurei (specimen number: 2017030085, Supplementary Figure 1F) from Haiwei Port, Haiwei Town, Changjiang City, Hainan Province, China (19°26′52″ N, 108°50′9″ E), and S. henslowianum (specimen number: 2017050197, Supplementary Figure 1G) from Putian, Fujian Province, China (25°2′23″ N, 119°7′50″ E) were collected. The thalli were cultivated in sterilized filtered seawater under fluorescent light at 24–26°C (80–110 μmol photons m^–2 s^–1; 12 h light/dark cycles) and washed in autoclaved seawater before being dried with paper towels. All the samples were stored in the Culture Collection of Seaweed at the Ocean University of China.

Genomic DNA was extracted from approximately 1 g of fresh tissue using the reformative cetyl trimethyl ammonium bromide (CTAB) method (Sun et al., 2018). β-mercaptoethanol was used to suppress the oxidation of phenols, and two concentrations of CTAB (2% and 10%) were used for the removal of polysaccharides. The quality and quantity of DNA were determined using a NanoDrop ND1000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, United States).

Approximately 5 μg of purified DNA from each species was used to construct three short-insert libraries, according to the instructions of the manufacturer (Illumina Inc., San Diego, California, United States). Approximately 9 Gb of paired-end sequencing data was used as input in the SOAPdenovo (Luo et al., 2012), based on the default assembly parameters, to assemble the plastid genomes. The plastid genome of S. fusiforme (GenBank accession number: NC_048511) was used as the reference sequence to determine the proportion of the plastid-related contigs. Subsequently, all plastid-related contigs were aligned and ordered into a circular structure using the CodonCode Aligner (CodonCode Corporation, Centerville, Massachusetts, United States) (Liu N. et al., 2019).

The protein-encoding genes, ribosomal (r)RNA genes, and transfer (t)RNA genes of the seven species were annotated based on those of S. fusiforme using Geneious R10 (Biomatters Ltd., New Zealand; available on http://www.geneious.com/). The tRNA genes were predicted using tRNAscan-SE version 1.21² (Schattner et al., 2005). The physical mappings of the seven plastid genomes were obtained using Organellar Genome DRAW (OGDRAW) version 1.3.1³ (Greiner et al., 2019). The sequence alignment and base composition were conducted using MEGA 7.0, (Sudhir et al., 2016). The plastid genomes of the seven reported species from order Fucales were aligned using the ascensive Mauve genome aligner under the default settings in the Geneious software (Drummond, 2012).

The phylogenetic analysis of the class Phaeophyceae was conducted based on 73 shared plastid protein-encoding genes from 48 plastid genomes available in the GenBank database and the seven plastid genomes obtained in this study, with Eunotia naegelii (GenBank accession number: NC_024928) as the outgroup species. The protein sequences were aligned in MEGA 7.0, and the entire concatenated alignment was manually generated in BioEdit (Hall, 1999). The Gblocks server⁴ (Castresana, 2000) was used to remove the poorly conserved regions, and only 19,340 of the original 22,677 amino acids were retained in the concatenated protein alignments. To reconstruct the phylogenetic tree, the best-fitting model for the maximum-likelihood (ML) was selected using ProtTest 3.4.2 (Darriba et al., 2011), and subsequently, the phylogeny was conducted by ML bootstrap analysis using RAxML v8.2.12 (Stamatakis, 2006), with 1,000 replications under the JTT + G + F model. Bayesian inference (BI) was performed using MrBayes version 3.1.2 (Huelsenbeck and Ronquist, 2001), and the phylogenetic analysis was performed using two independent runs with four Markov chains, which ran for 1,000,000 generations until the average SD of split frequencies was below 0.01. The output trees were sampled every 100 generations. In addition, the first 25% of the samples were removed as residue, and the rest were used to build a 50% majority rule consensus tree accompanied by posterior probability values. The phylogenetic tree was plotted in FigTree version 1.3.1⁵ (Rambaut, 2009).

The MCMCTREE in the PAML4 (Yang, 2007) software package was used to estimate the species divergence time, based on the “global clock” molecular clock model and the “HKY85” nucleic acid replacement model. Additionally, a total of four calibrations were selected: Sargassaceae vs. Fucaceae, 49–66 million years ago (Mya); Laminariaceae vs. Lessoniaceae, 25.0–107.0 Mya; Scytosiphonaceae vs. Ectocarpaceae, 87.0–114.6 Mya; Laminariales vs. Ectocarpales, 129–220 Mya. The calibration points were obtained from the Timetree website⁶.

The seven plastid genomes were assembled as single circular molecules, with sizes ranging from 124,258 to 124,563 bp (Table 1). The average GC contents of the plastid genomes were 30.49% (S. ilicifolium var. conduplicatum), 30.43% (S. graminifolium), 30.42% (S. phyllocystum), 30.40% (S. muticum), 30.35% (S. feldmannii), 30.43% (S. mcclurei), and 30.44% (S. henslowianum), which were similar to those previously reported for Fucales species. All seven complete plastid genomes were submitted to GenBank under the following accession numbers: MW767830 (S. ilicifolium var. conduplicatum), MW784163 (S. graminifolium), MW784165 (S. phyllocystum), MW784166 (S. muticum), MW784167 (S. feldmannii), MW784168 (S. mcclurei), and MW784169 (S. henslowianum). The general features of the five previously reported complete plastid genomes from order Fucales, including three species from Sargassum (S. fusiforme, S. horneri, and S. thunbergii), one species from Coccophora (C. langsdorfii), and one species from Fucus (F. spiralis) are summarized in Table 1 for a comprehensive comparison.

The seven plastid genomes were predicted to encode a total of 173 genes, including 139 protein-encoding genes, 28 tRNA genes, and 6 rRNA genes, which were the same as those of the previously reported Fucales species (Table 1). The seven genomes were conserved in the gene distribution on both the heavy and light strands (Figure 1). Most of the tRNA genes ranged from 71 to 90 bp in size, but one tRNA gene of trnL was 294 bp in size in all seven plastid genomes. For the seven plastid genomes, all tRNA and rRNA genes retained a single copy and no introns were found. However, S. thunbergii and C. langsdorfii held two IR regions (IRa and IRb) and one intron, and S. horneri also had one intron. In addition, similar to the majority of the Fucales species, tmRNA and miscRNA were not identified in our seven plastid genomes. In comparison with the previously published plastid genomes of the Phaeophyceae species, we confirmed that there was high conservation in the gene content and genome organization.

The coding sequence regions in the seven plastid genomes ranged from 107,070 to 107,154 bp, accounting for 86.02–86.14% of the corresponding plastid genomes. The total spacer sizes in the seven plastid genomes were 17,344, 17,490, 17,387, 17,247, 17,179, 17,257, and 17,428 bp, which were smaller than those in S. thunbergii (17848 bp), C. langsdorfii (17,846 bp), and F. spiralis (17,913 bp). This indicated that the seven plastid genomes from our study, together with those of S. fusiforme and S. horneri from the family Sargassum, held a relatively more compact organization. As insertions and deletions were rare, the size variations of the seven plastid genomes were mainly due to the differences in the intergenic region or non-coding region.

The single nucleotide polymorphism (SNP) sites were statistically analyzed, and the Ka/Ks values (Supplementary Table 1) were calculated to test the gene selection. Based on the analysis of 139 protein-encoding genes in the seven species, there were almost no non-synonymous changes in most of the genes, as with the value of synonymous changes of less than 1, it resulted in a final Ka/Ks value of much less than 1. Therefore, it was concluded that the genes of the seven Sargassaceae species were subject to purification selection, which suggested that in most cases, selection eliminates harmful mutations and keeps the protein unchanged. It shows that the plastid genome sequences of the seven species were relatively conservative.

The genomes of seven species in this study were compact, and there were five pairs of genes overlapping each other, where the common regions of ftrB-ycf12, sufC-sufB, rps1-thiS, psbC-psbD, and rpl123-rpl4 were 6, 4, 4, 53, and 8 bp, respectively. These five overlaps were also found in the plastid genomes of the other Fucales species, suggesting that the five overlaps were conserved in the order Fucales. The overlapping regions of rpl23-rpl4 were also observed in some plastid genomes of red algae and diatoms, indicating its high conservation in algae.

The seven plastid genomes in this study were mapped as a canonical quadripartite structure with two IRs, which divide the circular molecule into an SSC (39,892 to 40,012 bp) and an LSC (73,500 to 73,693 bp) (Figure 1). Among them, S. phyllocystum, S. mcclurei, and S. henslowianum had the same length of IRs (5,421 bp), which was different from those of S. graminifolium (5,429 bp), S. feldmannii (5433 bp), S. muticum (5,436 bp), and S. ilicifolium var. conduplicatum (5,419 bp) (Table 1).

In terms of the predicted protein-encoding gene content, Venn diagrams of the gene contents revealed the shared and/or unique genes in the plastid genomes of the Fucales (Supplementary Figure 2). In total, the seven plastid genomes in our study encoded 139 genes, and 136 protein-encoding genes were shared by the 12 Fucales plastid genomes, accounting for 97.84% of all plastid genes. The other three genes, orf76, orf501, and atpI, were shared between the seven species in our study and C. langsdorfii/S. horneri, while the seven plastid genomes had lost four genes, orf470, orf467, orf219, and atpl. The atpl gene was found only in the C. langsdorfii plastid genome. The plastid genomes of S. fusiforme and S. horneri contained the unique genes orf470 and orf467, respectively, and both covered the orf219 gene. Several ancient genes, such as glnB, were absent in the plastid genomes of 12 species in Fucales.

Nearly all protein-encoding genes used ATG as the start codon in the plastid genomes of S. ilicifolium var. conduplicatum (135 genes, 97.12%), S. graminifolium (136 genes, 97.84%), S. phyllocystum (135 genes, 97.12%), S. muticum (136 genes, 97.84%), S. feldmannii (135 genes, 97.12%), S. mcclurei (135 genes, 97.12%), and S. henslowianum (135 genes, 97.12%). GTG, ATT, TTA, and TTG were also used as start codons in addition to ATG. GTG was used as the start codon for psbF and rpl3 genes in the plastid genomes of S. ilicifolium var. conduplicatum, S. graminifolium, S. phyllocystum, S. muticum, S. feldmannii, S. mcclurei, S. henslowianum, S. fusiform, and S. thunbergii, and gene psbF in S. horneri and F. spiralis. ATT was used as the start codon for genes rps8 and rpl3 in C. langsdorfii and gene thiS in the other six genomes. The gene rps18 used TTA as the start codon in S. ilicifolium var. conduplicatum, S. phyllocystum, S. muticum, S. feldmannii, S. mcclurei, S. henslowianum, and used TTG in S. thunbergii (Table 1).

All three typical stop codons (TAA, TAG, and TGA) were identified with an obvious preference for TAA, based on the high percentage values in the plastid genomes of S. ilicifolium var. conduplicatum (109, 78.42%), S. graminifolium (108, 77.70%), S. phyllocystum (108, 77.70%), S. muticum (106, 76.26%), S. feldmannii (110, 79.14%), S. mcclurei (110, 79.14%), and S. henslowianum (111, 79.86%). The rbcR gene used CAA as a stop codon in S. ilicifolium var. conduplicatum and S. phyllocystum.

We further conducted a detailed co-linear analysis for the gene order in the 12 complete plastid genomes of order Fucales and confirmed strong collinearity with similar gene orders and genomic organizations (Figure 2). In Figure 2, the homologous regions are marked with similar colors and are linked by a vertical line. The white blocks indicate genomic regions with high variations, for example, the right end of the yellow blocks. Two IRs (IRa and IRb) are highlighted in red and boxed in all 12 genomes. There were no differences in gene sequences among the IR region in the Sargassum species. Two inversions covering two tRNA genes, trnC(gca) and trnN(gtt), were identified in 10 Sargassum species, against the other two Fucales species (C. langsdorfii and S. horneri). The trnC (gca) gene is on the forward strand in S. horneri and C. langsdorfii, but is on the reverse strand in the other species. The trnN (gtt) gene is forward in S. horneri, but is reversed in the others. The region is the only structural variation in the 12 plastid genomes, as the plastid genomes of Fucales species are highly conserved.

In addition to the genus Sargassum, we downloaded the complete plastid genome sequences of the class Phaeophyceae together with the other selected algae that were available in the National Center for Biotechnology Information (NCBI), to perform phylogenetic analysis and elucidate the global phylogeny. We collected protein-encoding sequences from a total of 55 algal plastid genomes (Supplementary Table 3). Finally, we obtained 73 common genes to construct the ML and BI phylogenetic trees, with Eunotia naegelii as the outgroup (Figure 3A). The results showed that the topological structure of the phylogenetic trees inferred by the ML and BI methods were highly concordant and divided all species into five clades corresponding to their orders: Desmarestiales, Dictyotales, Ectocarpales, Fucales, and Laminariales. In the ML and BI trees, high bootstrap support and posterior probability values were obtained, except for the clade containing Hapterophycus canaliculatus and Scytosiphon lomentaria, and the clade containing Egregia menziesii and Ecklonia arborea.

We further performed a phylogenetic analysis based on 125 common genes of 18 species of the order Fucales, and the phylogenetic tree divided them into three clades, one for 15 Sargassum species, one for Fucus, and another for C. langsdorfii. Fifteen Sargassum species were divided into two subclades, and this phylogenetic tree was strictly stable with high bootstrap support and posterior probability values (Figure 3B).

The brown algae chronogram was inferred based on nucleotide sequences (Figure 4). Firstly, the Phaeophyceae node was estimated to be formed about 201.42 Mya, [with a 95% highest probability density (zHPD): 135.29–237.93 Mya]. The Fucales node was estimated to be formed about 92.52 Mya (with a 95% HPD: 62.05–109.16 Mya), and the divergence of the Sargassum node from Coccophora langsdorfii was 28.25 Mya (with a 95% HPD: 18.93–33.38 Mya). The divergence time of S. ilicifolium var. conduplicatum and S. graminifolium was 4.55 Mya (with a 95% HPD: 3.05–5.45 Mya), and the node of S. feldmannii was 9.51 Mya (with a 95% HPD: 6.37–11.27 Mya). The node of S. phyllocystum and S. mcclurei was.07 Mya (with a 95% HPD:0.03–0.12 Mya), and the neighbor node of S. henslowianum was 2.14 Mya (with a 95% HPD: 1.42–2.62 Mya).

Generally, the plastid genome contains a typical pair of IR regions (Henry, 2005), which are often observed in the plastid genomes of advanced plants and green algae (Cao, 2015). But the IR loss has occurred in some algae, for example, Rhodomonas salina (Khan et al., 2007), Chlorella vulgaris (Wakasugi et al., 1997), Ulva fasciata (Melton and Lopez-Bautista, 2015), and U. linza (Wang et al., 2017). It has been reported that the existence of IRs increases the capacity of ribosomal genes (Bock and Knoop, 2012). The expansion and contraction of IR regions in the plant plastid genome is a common phenomenon that contributes to the plant plastid genome size (Huang et al., 2014). For example, green algae of both Dunaliella salina (Smith et al., 2010) and Chlamydomonas reinhardtii (Dron et al., 1982) contained IRs with sizes of 14.4 kb and 22.2 kb, respectively, which resulted in larger plastid genomes than the algae without IRs. To our best knowledge, the IR regions were found in the chloroplast genomes of all the brown algae. However, compared with green algae, the IR regions of brown algae are relatively small. For example, the IR sizes in Vaucheria bursata (Linne and Kowallik, 1992), Phaeodactylum tricornutum, Odontella sinensis (Secq et al., 2007), Sargassum fusiforme (Zhang et al., 2019), and Sargassum thunbergii (Yang et al., 2016) are 5.9, 6.9, 7.7, 5.3, and 5.4 kb, respectively (Supplementary Table 2). We revealed the high conservation of plastid genome structure and gene content and confirmed two IR regions (IRa and IRb) in Sargassum species (Table 1). Notably, IRa covers a part of the bases of the gene rpl21, while IRb covers a part of the gene cbbx. The protein-encoding genes of these seven species are under the pressure of purification selection, which indicates that most variations of the genes have been eliminated during the evolution process, and that the conservation is relatively high. Some studies believe that the gene conversion ability of the plastid genome is high, which can ensure the consistency and stability of the two IR region sequences and enhance the conservation of the plastid genome (Khakhlova and Bock, 2010). In comparison with the other published plastid genomes of red algae and green algae, our analysis pipeline revealed that most of the red algae have lost the IR region, and that the IR length of green algae is longer than that of brown algae (Table 1). It could be concluded that the seven species in this study experienced identical IR evolution events, which had occurred before the formation of order Fucales. This provides an updated understanding of the evolution of brown algae.

Some studies suggested that the evolution of the non-coding region duplication has obvious lineage specificity, while the coding region is conservative (Ku et al., 2013). The proportion of the repetitive sequences in the genome increases with the increase in the evolutionary level of an organism (Shu-Fen et al., 2017). These repeats contain a great deal of genetic information and play an important role in plastid genes (Cavalier-Smith, 2002). We found that the brown algal plastid genomes are relatively small in comparison with the plastid genomes of red algae (Liu et al., 2016; Yang et al., 2016; Zhang et al., 2018). After analyzing the plastid genomes of the Sargassum species, we found that there were few repeats in the coding regions, and most of the duplicates existed in the non-coding regions. Among them, three rRNAs were repeated in pairs, and 50% of the rRNAs are replicated. The plastid genomes of the seven species in this study consisted of 173 genes, of which 170 genes (6 rRNAs, 28 tRNAs, and 136 protein-encoding genes) were shared in the plastid genomes of the 12 species of order Fucales. It is worth noting that only the C. langsdorfii plastid genome lost the atpI gene, which is present in the plastid genome of most marine brown algae. This gene deletion might be due to gene loss or horizontal transfer to the nucleus (Feng and Shaojun, 2016). In addition, we found that some ancient genes were missing from the plastid genomes of the Sargassum.

Among all the available plastid genomes from the Fucales species, the majority of the protein-encoding plastid genes used ATG as the start codon, followed by GTG, ATT, TTA, and TTG. In addition to the standard start codon ATG, GTG was another commonly used start codon, which was mainly present in the bacteria, and the analysis of the organellar genomes revealed that several plastid and mitochondrial genes used GTG as the start codon in both red and brown algae (Corguillé et al., 2009; Wang L. et al., 2013; Wang X. et al., 2013; Zhang et al., 2013; Liu N. et al., 2019). The start codon GTG was utilized in the species of Fucales, showing a conserved evolutionary process of these genes. Usually known as an unusual start codon in eubacteria and archaea, TTG was located in a frame upstream from the originally proposed ATG start codon (Golderer et al., 1995). TTG as start codon had been found in S. thunbergii.

Despite the strong collinear conversation in the 12 reported plastid genomes from the order Fucales, we still identified an inversion of trnC(gca) in S. horneri and C. langsdorfii, and trnN(gtt) in S. horneri. This indicates the evolutionary positions of S. horneri and C. langsdorfii in the order Fucales. In plastid genomes, gene inversions have been reported among the various algae orders (Liu N. et al., 2019). Inversions and rearrangements of large fragments have been reported in red algae (Zhang et al., 2020), 30 kb inversion has been found in all living vascular plants (Raubeson and Jansen, 1992), and two transient reversals were shared in all ferns (Roper et al., 2007; Gao et al., 2009).

The Phaeophyceae species are a group of multicellular marine eukaryotes that belong to the stramenopiles, and their plastids had originated from secondary endosymbiosis, in which a red alga was engulfed by a non-photosynthetic protist (Keeling, 2010). It was supposed that plastid genomes could provide more evolutionary information than mitochondrial genomes (Kuipeng et al., 2018), and brown algal plastid genomes appear to be highly rearranged in genome architectures (Feng and Shaojun, 2016).

According to our time-calibrated phylogeny, most brown algal orders had diversified in the timespan ranging through the Triassic period (251–119 Mya). For instance, the diversification of Fucales, Laminariales, and Ectocarpales was inferred to occur around 17.60 Mya. Brown algae have a particularly scarce fossil record (Draisma et al., 2003), and many of the fossils cannot be confidently assigned to brown algae (Clayton, 1988). This resulted in difficulty in identifying suitable algal fossils to elaborate a time-calibrated phylogeny (Verbruggen et al., 2009). Therefore, we decided to include four calibration scales that were based on fossils with clear evidence of a relationship with extant brown algal lineages. The origin of Sargassum was estimated to be more recent in the present study, i.e., 28.25 Mya (with a 95% HPD: 18.93–33.38 Mya), which was similar to the estimate of 22 Mya as the lower bound of the 95% HPD by Silberfeld et al. (2010), but was more ancient compared with the estimate of 6.7 Mya as the lower bound of the 95% HPD by Zhi et al. (2020).

Some previous studies have divided Ochrophytes algae into three branches: SI, SII, and SIII (Yang et al., 2012; Wetherbee et al., 2018; Bringloe et al., 2020), which is consistent with the results in a previous study (Bringloe et al., 2020). These phylogenetic analyses were mainly based on a few genes, such as SSU and LSU rDNA plus five proteins (Riisberg et al., 2009), and SSU rDNA plus four proteins (Yang et al., 2012). More genes are needed to construct reliable phylogenetic trees. In the Ochrophytes species, 34 plastid proteins were used to produce reliable phylogenetics (Ševèíková et al., 2015), and a collection of eukaryotic phylogenetic markers were used to assemble a phylogenomic matrix (Yang et al., 2012). The phylogenetic relationship of the Sargassum in our previous study was not well resolved based on only 13 concatenated chloroplast protein sequences (Jia et al., 2014). In our phylogenetic tree of this study, combined with the results of the co-linear analysis, we deduced that the genus Sargassum in the order Fucales had relatively significant differentiation in two subclades, which implied the possible non-monophyletic relationship of Sargassum. It also suggested that the species in the order Fucales form independent branches compared with other types of brown algae (Peters and Burkhardt, 1998). As there are still low bootstrap and posterior probability values for the individual branches with a lack of nuclear genome information, the phylogenetic relationships among the genus Sargassum could not be fully resolved based on the existing molecular data of the chloroplast genomes. In the future, the molecular phylogeny based on nuclear-encoded protein genes will provide more evidence for the evolutionary divergence of the Sargassum species, as their phenotypic differences were observed. To date, our results provide a complete and reliable maternal phylogenetic tree, and more whole genome sequences will possibly be included in further studies on phylogenetic trees and the evolution of Sargassum in brown algae.