1 INTRODUCTION

The red algae (Rhodophyta) have been divided into two subphyla, the Rhodophytina and Cyanidiophytina, which are further classified into seven classes, the Bangiophyceae, Compsopogonophyceae, Florideophyceae, Porphyridiophyceae, Rhodellophyceae, Stylonematophyceae and Cyanidiophyceae (Yoon et al., 2006). The Bangiophyceae include only one order, Bangiales, and one family, Bangiaceae, which has been suggested to be monophyletic (Müller et al., 2001). Sutherland et al. (2011) revised the order Bangiales to include at least fifteen genus-level taxa (Bangia, 'Bangia' 1, 'Bangia' 2, 'Bangia' 3, Boreophyllum, Clymene, Dione, Fuscifolium, Lysithea, Minerva, Miuraea, Porphyra, Pseudobangia, Pyropia, and Wildemania). Bangiales have been classified as an ancient lineage, based on a 1 200 million-year-old (Ma) fossil of Bangiomorpha pubescens, a Bangialike taxon (Butterfield, 2000). However, Porphyralike fossils have been found in rocks dated at approximately 425 Ma (Blouin et al., 2011). Yoon et al. (2004) estimated the evolutionary timeline using a six-gene data set, and suggested ca. 800 Ma for the divergence of the Florideophyceae from the Bangiales. However, the multicellular holdfast of Bangiomorpha pubescens has greater similarity to Erythrotrichia than to Bangia. It is possible that Bangiomorpha evolved as an extinct lineage in parallel to the Bangiophyceae and Compsopogonophyceae (Yang et al., 2016). Although previous molecular studies have produced robust phylogenies of the Bangiales, divergence time estimation of this order has rarely been addressed.

Species of Bangiales include the most economically important marine crops grown worldwide from tropical to polar seas. Previous research suggested that the center of diversity of Bangiales was in the Northern Hemisphere, particularly the North Pacific (Conway et al., 1975). As more taxa have been described worldwide, New Zealand, in the Southern Hemisphere, has more recently been regarded as a center of Bangiales diversity (Broom et al., 2004). Aspects of the historical distributions of this group may be of interest, such as: the origin of the lineage, ancestral distribution and the factors that shaped the composition and distribution of this clade over time. Many studies have addressed biogeographic questions using phylogenetic analyses, molecular dating, and reconstruction of ancestral biogeography (Dumont et al., 2005; Deng et al., 2015; Berger et al., 2016). However, to date, no molecular systematic study has been performed to determine the historical biogeography of the Bangiales.

The Bangiales comprise approximately 190 species, only 0.43% of currently described species of Rhodophyta. Florideophyceae is the most taxon-rich class of the Rhodophyta (95% of described Rhodophyte species) (Guiry and Guiry, 2017). Species diversification represents the balance between speciation and extinction rates that gave rise to the extant diversity of a clade. Studies have suggested that diversification is somehow linked to rates molecular evolution over time (Barraclough and Savolainen, 2001; Webster et al., 2003; Davies et al., 2004). Various types of abrupt genomic disruption can play an important role in promoting speciation via reproductive isolation (Venditti and Pagel, 2010). Collinearity analysis of Bangiales plastid genomes indicated a highly conservative genome structure (Janouškovec et al., 2013; Lee et al., 2016), but the evolutionary rate at the nucleotide level in this group is unknown. In accordance with the central dogma of molecular biology, the protein-coding gene mutation rate has classically been regarded as a major unit of evolution. A previous study predicted that the evolutionary rate at neutral sites corresponds to the actual mutation rate if synonymous sites are assumed to be neutrally evolving (Kimura, 1984). Studies on the relative synonymous-site substitution rates (dS) of plastid DNA (ptDNA), mitochondrial DNA (mtDNA) and nuclear DNA (nDNA) from different lineages have provided insight into genome evolution (Buschiazzo et al., 2012; Leliaert et al., 2012; Smith, 2015). Moreover, in pairwise comparisons of orthologous genes, ratios of dN to dS (dN/dS) give an indication of the mode and strength of selection, where dN/dS > 1 suggests adaptive selection and dN/dS < 1 indicates purifying selection. The rate of molecular evolution may help us to better understand the speciation and environmental adaptation of Bangiales.

In this study, all available Bangiales organellar genomes (10 ptDNA and 11 mtDNA genomes) as well as the plastid RUBISCO LSU (rbcL) and nuclear SSU rRNA (nrSSU) gene sequences of 162 Bangiales species were downloaded from the GenBank database to 1) reconstruct phylogenetic relationships among the representative species, 2) estimate divergence times, 3) clarify the historical biogeography, and 4) identify the relative substitution rates of the Bangiales species.

2 MATERIAL AND METHOD 2.1 Phylogenetic analyses and divergence time estimation

Organellar genome data were downloaded from the GeneBank database (Table S1). To determine the phylogenetic relationships among Bangiales species, the concatenated protein-coding amino acid sequences from both plastid and mitochondrial genomes were used to construct phylogenetic trees (Table S2). We used the all-versus-all blast to select single copy homologous genes with the e-value 1e-5. These protein-coding genes are shared by all of the selected taxa. The Chondrus crispus genome was used as the outgroup. Genome sequences were aligned using the program MAFFT version 5 and adjusted manually (Katoh et al., 2005). The aligned sequences were further trimmed using trimAl with the option 'automated1' (Capella-Gutiérrez et al., 2009). Maximum likelihood (ML) analyses were conducted using RaxML-8.2.4 (Stamatakis, 2014). The best model and parameter settings were chosen according to the Akaike information criterion using ProtTest 3.0 for the ML analyses (Abascal et al., 2005). The ML searches relied on the cpREV+G+I substitution model for plastid sequences and the JTT+G+I substitution model for mitochondrial sequences (-f a, 1 000 bootstrap replicates). The MP trees were constructed with PAUP^*4.0b10 (Swofford, 2003). Heuristic searches were conducted using random taxon addition with branch swapping tree bisection and reconnection (TBR). Non-parametric bootstrap analysis was performed with 1 000 replicates with TBR branch swapping. Bayesian analyses were carried out using MrBayes3.2 with the best ProtTest model determined as noted above (Huelsenbeck and Ronquist, 2001). Four independent Markov Chain Monte-Carlo chains were run simultaneously and sampled every 100 generations for a total of 1 000 000 generations. The first 10% of the trees were discarded as "burn-in."

To calculate the divergence times of Bangiales species, a dataset of 66 protein-coding nucleotide sequences from 27 Rhodophyte plastid genomes were selected (Table S2). We used BEAST v1.8.3 with the uncorrelated lognormal relaxed clock model to estimate rate change (Drummond et al., 2012). The best model inferred by ModelTest3.7, GTR+G+I, was used (Posada and Crandall, 1998). The Yule tree prior and unlinked clock models were used across gene regions. We set the most recent common ancestor with a lognormal prior, an offset of 950 Ma, and a standard deviation of 25.0 based on the divergence of the Florideophyceae and Bangiophyceae (Herron et al., 2009; Yang et al., 2016).

2.2 Estimation of substitution rates

The bidirectional best hit algorithm (BDBH) was used to compile clusters of potential 1:1 orthologous sequences using BLAST with the -e threshold=1e-10 in GET_HOMOLOGUES v2.0.3 (Contreras-Moreira and Vinuesa, 2013). Selected core gene amino acid sequences were aligned with MAFFT and then converted to nucleotide sequences using a Perl script written in house to achieve the codon alignment pattern. Pairwise distances at non-synonymous (dN) and synonymous (dS) sites were estimated for individual genes (Pyropia endiviifolia vs. Pyropia kanakaensis, Porphyra purpurea vs. Porphyra umbilicalis, Cyanidiaceae sp. MX-AZ01 vs. Cyanidioschyzon merolae, Gelidium elegans vs. Gelidium vagum, Gracilaria chilensis vs. Gracilaria salicornia) using codeml (PAML 4.9b) (Yang, 2007), with the maximum likelihood method employed with settings seqtype=1, CodonFreq=2, and Runmode=-2. Rates of substitution at coding sites (μ) were calculated using the following formula:

where d represents the distance at synonymous (dS) or non-synonymous (dN) sites and T represents the divergence time between pairwise comparisons (Buschiazzo et al., 2012). Divergence times were estimated as described above.

2.3 Reconstruction of ancestral areas

For the reconstruction of ancestral areas, rbcL and nrSSU gene sequences of 162 Bangiales species were downloaded from GenBank (Table S3). The Bayesian analysis was conducted as described by Sutherland (Sutherland et al., 2011). Seven biogeographic areas were defined based on the collection information: (A) Southern South America and Antarctica; (B) New Zealand and Australia; (C) Southern Africa; (D) Eastern North America; (E) Western North America; (F) Eastern Asia; (G) Europe. The ancestral area of Bangiales was evaluated using dispersal-vicariance analysis (S-DIVA) by RASP software with two as the maximum number of ancestral areas at each node (Ronquist and Cannatella, 1997; Yu et al., 2010).

3 RESULT 3.1 Phylogenomic analysis

Trees were constructed using a dataset of 153 amino acid sequences from 11 plastid genomes and 19 amino acid sequences from 13 mitochondrial genomes to illustrate the evolutionary relationships within Bangiales, and most of the nodes were inferred with strong support using the MP, ML and BI methods. The ptDNA tree showed that the clade uniting Py. endiviifolia formed a sister group with Py. kanakaensis, and Py. yezoensis clustered with Py fucicola. The remaining Py. haitanensis and Py. perforate grouped together. Moreover, Bangia and Porphyra were assigned to the same clade supported by all three methods and they were added to the sister group with Pyropia, leaving Wildemania at the basal node of Bangiales in MP and BI analyses (Fig. 1a). However, the ptDNA ML tree showed that Wildemania formed a sister group with Pyropia first. The topology of mtDNA tree was similar with ptDNA tree (Fig. 1b). The added species Py. tenera was clustered with Py. yezoensis, and then assigned to Py. fucicola with high support. Py. nitida was clustered with Py. endiviifolia and then assigned to Py. kanakaensis. The mtDNA tree showed Wildemania was at the basal node of Bangiales supported by all three methods. On the whole, our results support the (((Bangia, Porphyra) Pyropia) Wildemania) topology.

3.2 Divergence time estimation and reconstruction of ancestral areas

Based on the analysis of aligned nucleotide sequences of Rhodophyte plastid genomes, we estimated divergence times. The results revealed that the multicellular red algae originated in the late Mesoproterozoic to early Neoproterozoic at approximately 1 001.87 Ma (95% HPD=910.27–1 281.12 Ma, PP=0.99). Bangiales diversified ca. 246.40 Ma (95% HPD=194.78–318.24 Ma, PP=0.99), in the Permian and Triassic (Fig. 2a). The Pyropia clade diversified at 121.37 Ma (95% HPD=111.31– 131.40 Ma, PP=1), in the Early Cretaceous, and Porphyra began to diverge from Bangia at 180.51 Ma (95% HPD=165.73–195.31 Ma, PP=1), in the Jurassic (Fig. 2b). In this study, the species Py. endiviifolia diverged from Py. kanakaensis at 68.07 Ma (95% HPD=62.01–73.91 Ma, PP=1), in the Late Cretaceous.

Ancestral area reconstruction within the Bangiales was conducted using S-DIVA (Ronquist and Cannatella, 1997). Because of the limited organellar genome data, which contained very few Bangiales species from the Southern Hemisphere, the biogeographic analysis was conducted using the plastid rbcL and nrSSU gene sequences of 162 widespread Bangiales species together with the corresponding distribution information (Fig. 3, Table S3). The analysis revealed that Bangiales had a center of origin in eastern Gondwanaland, in the current area encompassing New Zealand and Australia (Oceania) (S-DIVA=0.77). Moreover, the ancestral area of the subgroup Pyropia spread to Eastern Asia (S-DIVA=0.88). Porphyra dispersed to Oceania (S-DIVA=0.85) and "Bangia" dispersed to either Oceania (S-DIVA=0.85) or Western North America (S-DIVA=1.00). The details of the ancestral area distribution as determined based on S-DIVA analysis are shown in Fig. 4.

3.3 Substitution rates in Bangiales

To better understand the evolutionary rate of the order Bangiales, protein-coding genes from ptDNAs and mtDNAs were aligned, and the numbers of pairwise synonymous (dS) and non-synonymous (dN) substitutions per site were inferred (Table 1). The mean dS values were 0.528±0.242 (ptDNA) and 1.623±0.699 (mtDNA) within Pyropia, and 0.463±0.193 (ptDNA) and 1.738±0.581 (mtDNA) within Porphyra. The rates of substitution at dN were lower than those of dS. The dN/dS values were significantly smaller than 1, with averages ranging from 0.035 to 0.044, which suggest the existence of purifying selection acting on the organelle proteincoding genes of the Bangiales species.

Combining substitution rates with the divergence times (Table 2), we estimated the substitution rates per site per year at synonymous sites (μ_S). The mean μ_S in Bangiales ranged from 0.356×10^-8 to 0.388×10^-8 (ptDNA) and from 1.193×10^-8 to 1.337×10^-8 (mtDNA). The average synonymous mutation rates of ptDNA were all lower than those of other groups of red algae. However, compared with other lineages outside the red algae, Rhodophytes had dramatically higher μ_S, by more than 3-fold for ptDNA and 23-fold for mtDNA relative to land plants, and by more than 5-fold for ptDNA and 15-fold for mtDNA relative to green algae. Consequently, the mean dN/dS value was lower in algae than in land plants, 2.5-fold on average, which implies different levels of environmental adaptability.

4 DISCUSSION 4.1 Phylogenetic relationships within Bangiales taxa

Studies have shown that multigene phylogenies elucidate phylogenetic relationships more precisely when the different evolutionary rates of the different genes are considered (Yoon et al., 2006; Verbruggen et al., 2010). Several studies based on complete plastid or mitochondrial multigene sequences showed great potential to elucidate the phylogenetic relationships of red algae (Reith and Munholland, 1995; Ohta et al., 2003; Hagopian et al., 2004; Verbruggen et al., 2010; Wang et al., 2013; Campbell et al., 2014; Yang et al., 2015). We utilized the entire available set of ptDNA and mtDNA genes to explore the phylogenetic tree of the Bangiales. In this study, the species B. fuscopurpurea was ascribed to the clade 'Bangia', which was closer to Pyropia 2, based on the combined data sets of nrSSU and rbcL genes (Fig. 4). Nevertheless, phylogenomics based on information from the entire genome revealed that the phylogenetic relationship between Bangia and Porphyra was closer than that of either with Pyropia, with high bootstrap support. Regarding the position of the genus Wildemania, our results indicated with high support that this clade diverged at the base of the Bangiales, which was in accordance with the results of Yang et al. (2015) based on whole mitochondrial genomes. Based on plastid genes, Hughey (2016) proposed that Wildemania was more closely related to Pyropia. Such inconsistency may occur when taxon sampling is sparse and when substitution saturation is present (Zhao et al., 2016). Therefore, more taxon information should be obtained; further works that combine the nuclear, plastid and mitochondrial genomes could help to better understand the relationships of Bangiales groups.

4.2 Divergence times and biogeography

Estimating divergence times using molecular data and fossil constraints is an advantage of the evolutionary study of red algae. There have been several studies on the divergence time of the taxonomically diverse class Florideophyceae (Xiao et al., 1998; Saunders and Hommersand, 2004; Xiao et al., 2004). The multicellular filamentous red alga Bangiomorpha pubescens from the ca. 1 200 Ma Hunting Formation (Somerset Island, Arctic Canada) was considered the oldest known taxonomically resolved eukaryotic fossil (Butterfield, 2000). However, the multicellular holdfast of B. pubescens had greater similarity to Erythrotrichia than to Bangia. It is possible that Bangiomorpha evolved as an extinct lineage in parallel to the Bangiophyceae and Compsopogonophyceae (Yang et al., 2016). A recent study indicated that Rafatazmia and Ramathallus represent crown-group multicellular Rhodophytes, which would antedate the oldest previously accepted red algae in the fossil record by about 400 million years (Bengtson et al., 2017). In our work, the organellar genomes of multicellular Rhodophytes include only Florideophyceae and Bangiophyceae. The latest common ancestor dating of the Florideophyceae and Bangiophyceae may not be sufficient to represent the dating of multicellular Rhodophytes. The present study provides molecular evidence of divergence times and historical biogeographic information, which indicate that Bangiales diversified approximately 246.40 Ma in the late Permian and early Triassic. Based on the present results, Bangiales is not as ancient as other red algae groups. This finding implies that this group experienced a long period of selection or an extinction event. Further genomic information of long-branch species such as Dione and Minerva is needed to more accurately evaluate the divergence time. At least three major diversification events occurred, which determined the four largest genera of the Bangiales (Bangia, Porphyra, Pyropia and Wildemania). We traced the ancestral species of Bangiales of eastern Gondwanaland, in the present New Zealand and Australia (Oceania) area of the Southern Hemisphere. At the time of these ancestral species, the land masses were still bound together into the vast supercontinent Pangaea. This emerging supercontinent presented severe climatic and environmental extremes because of its vast size (Chumakov and Zharkov, 2002). The south was cold and arid, with much of the region frozen under ice caps. Although climatic cooling was supposedly a primary cause of Permian mass extinctions (Stanley, 1988), Bangiales and other coldadapted species may have survived and expanded in the extreme environments of this period. Pangaea began to break apart in the mid-Triassic, forming Gondwanaland (South America, Africa, India, Antarctica, and Australia) in the south and Laurasia (North America and Eurasia) in the north (Ogg et al., 2004). These changes made it possible for plant migration across these continents. The ancestral area of the Bangia subgroup was determined to be in the Oceania area or western North America around the Jurassic. We infer that there may have been geographic isolation events at that time. The filamentous members of Bangia were divided into four groups: Bangia, 'Bangia' 1, 'Bangia' 2 and 'Bangia' 3 (Sutherland et al., 2011). Bangia and 'Bangia' 1 were deduced to originate from the Oceania area, and the other two genera were determined to be from western North America. The genus Pyropia diversified at 121.37 Ma in the Cretaceous and was widely distributed. In the Cretaceous, the breakup of the supercontinent Pangaea, which began during the Jurassic, continued. This change led to increased regional differentiation in species between the northern and southern continents (Casey, 1964). Based on our results, we speculate that Pyropia originated in East Asia and subsequently spread worldwide.

4.3 Evolutionary rate of Bangiales compared with other groups

A previous study implied that the high collinearity of plastid genomes indicate that the organellar genomes of red algae evolve slowly (Janouškovec et al., 2013). Smith et al. (2012) studied the relative rates of evolution among the three genetic compartments of the red alga Porphyra (mtDNA, ptDNA, and nDNA) and concluded that the mutation rates do not correlate with genome architecture. The present study was based upon whole plastid and mitochondrial genomic loci rather than on a small set of genes, and is relevant to understanding the differences in evolutionary patterns between Bangiales and other groups. First, we found that evolutionary rates are dramatically faster in red algae than in land plants and green algae. Secondly, evolutionary rates are slower in Bangiales than in Florideophyceae. To date, Bangiales is the only order that has been identified in the class Bangiophyceae, and includes ca. 190 species. However, at least 27 orders, including approximately 6 799 species, have been classified in Florideophyceae (Guiry and Guiry, 2017). The low speciation rate of the Bangiophyceae may be due to the relatively low evolutionary rates at the nucleotide level. In seed plants, gymnosperms are more conserved than angiosperms and have experienced very low speciation rates (Levin and Wilson, 1976). In birds, diversification has also been shown to be positively correlated with mutation rate (Lanfear et al., 2010). The lower evolutionary rates of Bangiales relative to Florideophyceae may be the cause of the lower species diversity.

Interestingly, land plants possess more species diversity than algae, despite generally lower mutation rates. The number of currently known, described and accepted land plant species exceeds 300 000, whereas the corresponding number of algae species worldwide is approximately 44 000 (Christenhusz and Byng, 2016). Species arise after splitting from ancestral species when they acquire new adaptations to a changing environment. We found that the mean dN/dS was approximately five times lower in algae than in land plants, despite much higher substitution rates in red algae genes (Table 1). To some extent, the fixation rate of mutations depends on the purifying (background) selection of the environment (Buschiazzo et al., 2012). Lower dN/dS in algae may arise from high selective constraints on nonsynonymous mutations because of the purifying selection of relatively constant marine life. The more variable terrestrial environment of land plants induces more adaptive mutations, which leads to a staggeringly high rate of fixation for non-synonymous mutations.

Simultaneously, we found that the average number of substitutions per synonymous site was higher (1.75–3.75 times) in mtDNA than in ptDNA in red algae, which was consistent with previously reported results (Robba et al., 2006; Milstein et al., 2008; Smith et al., 2012). This trend is also observed in Phaeocystis algae (Smith et al., 2014). In contrast, synonymous substitution rates are estimated at lower rates in mtDNA than in ptDNA (1:3) for seed plants and at roughly equal rates for green algae (Drouin et al., 2008; Hua et al., 2012). Seed plant cpDNA had lower conservation than mtDNA because of constraints on synonymous codon choice in the mitochondria, which may reduce the rate of synonymous substitution (Wolfe et al., 1987). The inverse situation exists in Rhodophyte species; plastid DNA exhibits a slower evolutionary rate than mtDNA. This difference in rates may be largely because of the conservation of ribosomal RNA regions of red algae, whereas contraction or expansion of the IR regions often occurred in seed plant chloroplast genomes (Kim and Lee, 2004).

5 CONCLUSION

In summary, our results indicate that the lower evolutionary rates of Bangiales relative to Florideophyceae may be the cause of their lower species diversity. The Rhodophyta had a faster evolutionary rate than the Viridiplantae, however, this clade also had a lower dN/dS ratio, which corresponds to adaptation to marine life. Phylogenomic analyses were conducted using complete organelle genomes of Bangiales, and the results support the topology (((Bangia, Porphyra) Pyropia) Wildemania) with high supporting values. In addition, divergence time dating and ancestral area reconstruction within the Bangiales were conducted for the first time. The increasing number of sequenced organelle genomes may allow further improvement of the understanding of the phylogenetics of the order Bangiales.

6 DATA AVAILABILITY STATEMENT

The authors declare that all data supporting the findings of this study are available in the appendix.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S3) is available in the online version of this article at https://doi.org/10.1007/s00343-018-7054-4.