1 INTRODUCTION

Ulva species are widely distributed worldwide, thriving in intertidal, brackish, estuaries, and even freshwater environments (Mantri et al., 2020), with more than 80 identified species documented in the Algaebase (Guiry and Guiry, 2021). The thallus of Ulva is composed of distromatic blade or monostromatic tube. As their morphological features are very limited and unstable, which are sensitive to various factors such as salinity (Blomster et al., 1998), temperature (Blomster et al., 2002), and associated bacteria (Kessler et al., 2018), the morphological identification for Ulva are always very difficult (Blomster et al., 1999). The development of molecular approaches has significantly improved this dilemma, resulting in the reconstruction of genera Ulva (Hayden et al., 2003), and identification of some cryptic species (Hofmann et al., 2010), but some of related species still lack appropriate molecular markers to distinguish them from each other (Kang et al., 2019; Steinhagen et al., 2019).

The type locations for U. linza Linnaeus 1753 and U. prolifera O. F. Müller 1778 were at Kent of England, and Danish island of Lorand respectively. Due to the lack of holotype of U. prolifera, the representative sequences of genetic markers for this species were originally derived from the samples collected in the British Isles that have been identified as "U. prolifera" on the basis of morphological characteristics (Blomster et al., 1998; Tan et al., 1999). However, Shimada et al. (2008) suggested that in the ITS-based phylogenetic tree, those U. prolifera collected in Japan were separated from European "U. prolifera", and they were almost completely indistinguishable with U. linza and U. procera (synonym of U. prolifera), forming a cluster named LPP. Cui et al. (2018) confirmed that the epitype of U. prolifera collected from the type location were also located in the LPP cluster, and suggested to revise the previous "U. prolifera" to U. splitiana. In addition to the genetic similarities between U. prolifera and U. linza, their distribution areas are often overlapped as well, including the Baltic Sea in Europe (Cui et al., 2018; Steinhagen et al., 2019), the Atlantic coast of North America (Guidone et al., 2013), and many parts of the Northwest Pacific (Shimada et al., 2008; Zhao et al., 2018). In general, U. linza is mainly spread in marine habitats, while strains of U. prolifera are found commonly in estuaries and brackish waters (Shimada et al., 2008; Ogawa et al., 2013). However, in the Southern Yellow Sea area, the two species grew intermixed and the biomass of both species is very high (Han et al., 2013). In particular, U. prolifera, as the dominant species, has caused the largest green tide in the world for consecutive years (Zhao et al., 2013). Thus, the accurate discrimination between these two related species have become necessary for investigations of Ulva species composition and green tide monitoring especially in this sea area.

Because molecular markers commonly used in Ulva, including ITS, rbcL and tufA, always failed to distinguish between U. prolifera and U. linza (Leliaert et al., 2009; Zhang et al., 2011; Xie et al., 2019), the 5S rDNA spacer region, which was polymorphic in individual, was developed (Shimada et al., 2008), and has been used widely to discriminate these two species (Hiraoka et al., 2011; Duan et al., 2012; Zhang et al., 2015; Song et al., 2019). From each species, this marker could generate multiple amplified products of different sequences and lengths, of which the smallest fragment of about 300 bp was considered to be specific to U. linza, and was not available in U. prolifera. However, this likely U. linza-specific genotype was later found in the epitype of U. prolifera as well (Cui et al., 2018), suggesting that the 5S rDNA spacer region was probably not a substantial interspecific marker (Melton III and Lopez-Bautista, 2021). New efforts focused more on the organelle genomes, since they have much richer polymorphic sites which are usually used for phylogenetic analysis among populations or related species (Yang et al., 2013; Zhang et al., 2021). Liu et al. (2020b) reported that a newly developed mitochondrial marker rps2-trnL can well distinguish four Ulva species including U. linza and the drifting U. prolifera causing the Yellow Sea green tide. However, the drifting U. prolifera has been revealed as a unique floating ecotype, which was clearly different from the widely-distributed attached populations, in terms of both genetics and the performances of reproductive isolation with U. linza (Hiraoka et al., 2011; Zhao et al., 2015), whether the usage of rps2-trnL can be extended to distinguish these two species still needs further verification with the attached populations of U. prolifera.

In this study, the chloroplast genome of a representative strain for attached U. prolifera was sequenced, and two existing chloroplast genomes which were from U. linza and the floating ecotype of U. prolifera respectively, were combined for a comparative analysis. The identified interspecific variations were used to develop new markers for the discrimination between these two related species.

2 MATERIAL AND METHOD 2.1 Seaweeds and molecular identification

Each Ulva strain used in this study was unialgal culture maintained in our laboratory, the collection information were shown in Supplementary Table S1. All the samples were cultured in Von Stosch's Enriched (VSE) medium renewed once a week, at 20 ℃ with a 12-h꞉12-h light (L)꞉dark (D) photoperiod and a photosynthetic irradiance of about 80 μmol photons/(m²·s).

Genomic DNA of each sample was extracted using a Plant Genomic DNA Extraction Kit (Tiangen Biotech Co. Ltd., Beijing, China) according to the manufacturer's instruction. The molecular identification for all samples were performed using ITS, 5S rDNA spacer, and a sequence characterized amplified region (SCAR) marker which was specific to the floating ecotype of U. prolifera dominating the green tide in the Yellow Sea. The primers and PCR procedures for ITS, 5S rDNA spacer, and SCAR markers referred to Leskinen and Pamilo (1997), Shimada et al. (2008), and Zhao et al. (2015) respectively. PCR products were sequenced in Ruibo Bio Tech Co. Ltd, Qingdao, China by a Genetic Analyzer (ABI3730XL, USA). Phylogenetic analysis were performed according to previous descriptions from Xie et al. (2020).

2.2 Chloroplast genome sequencing, assembly, annotation, and phylogenetic analysis

An attached U. prolifera sample U161 was selected as a representative for chloroplast genome sequencing. A single thallus was cut into segments for vegetative growth, then the algal tissue was sent to HengChuang Gene Co. Ltd. (Shenzhen, China) for high-throughput sequencing. Total genomic DNA was extracted using a Plant Genomic DNA Extraction Kit (Tiangen Biotech Co. Ltd., Beijing, China). The DNA library with an insert size of 350 bp was constructed using a library preparation kit (New England Biolabs Co. Ltd., USA) and sequenced using the Hiseq 4000 platform (Illumina Co. Ltd., USA) to obtain 150 bp ×2 paired-end reads. The low-quality sequences which are those with over 50% bases having quality values of Q < 19 or over 5% bases being 'N' were removed. The filtered reads were assembled into contigs by SOAPdenovo v2.04 (Luo et al., 2015), then aligned and ordered according to the reference genome. Last, raw reads were again mapped to the assembled draft chloroplast genome and the majority of gaps were filled through local assembly.

The chloroplast genome was annotated using program PGA (Qu et al., 2019). Ribosomal RNA genes (rRNAs) were identified by RNAmmer v1.2 (Lagesen et al., 2007), and transfer RNA genes (tRNAs) were searched using the tRNAscan-SE v2.0 (Chan and Lowe, 2019). The OGDRAW v1.3.1 was applied to draw the genome map (Greiner et al., 2019). The whole chloroplast genome sequence with annotation information was submitted to GenBank of NCBI using Bankit.

For phylogenetic analysis with whole chloroplast genomes, a total of 48 shared protein-coding genes among all available 26 chloroplast genomes of Ulva, including our data from U161 and other 25 which were obtained from NCBI as references, were selected for alignment by MAFFT v7.475 (Kuraku et al., 2013). After alignment and concatenating of the shared genes, the full length of 48 gene sequences were about 36 kb. The maximum likelihood (ML) phylogenetic tree with alignment sequences from 26 chloroplast genomes of Ulva was constructed using a GTR + G + I model and the sequence divergences were calculated with MEGA 6.0 (Tamura et al., 2013).

2.3 Comparative genomic analysis between U. linza and two ecotypes of U. prolifera

The complete chloroplast genomes of U. linza (NC030312), the floating U. prolifera (NC036137) collected from the Yellow Sea green tide, and the attached U. prolifera U161 (MZ571508), were used for comparative genomic analysis. The codon usage biases was analyzed using PhyloSiute v1.2.2 (Zhang et al., 2020) and codonW v1.4.4 (Meade et al., 1997). The collinearity analysis with these three chloroplast genomes was carried out to check the genome rearrangement by Mauve v2.4 with the ProgressiveMauve algorithm (Darling et al., 2010). Single nucleotide polymorphism (SNP) sites were searched by Mauve v2.4, and indel (insertion-deletion) sites were identified by Dnasp v5.1 (Librado and Rozas, 2009). In order to visualize structure variations across the genomes, the chloroplast genomic sequence comparative analysis were conducted using the mVISTA following a global pairwise alignment of the sequences with the LAGAN program (Frazer et al., 2004).

2.4 Development of new species-specific markers from chloroplast genomes

From the identified SNP, indels or structural variations, some of those regions that were homologous between the two ecotypes of U. prolifera but had obvious divergences between U. prolifera and U. linza were selected as molecular marker targets, and the flanking sequences at both ends, which were completely identical among the three chloroplast genomes, were used for design of species-specific primers using Primer 3.0. All primers were synthesized by Sangon Biotech (Shanghai) Co. Ltd. (Shanghai, China). The effects of species distinguish for designed primers were evaluated with each of twelve Ulva samples by PCR reactions. The profile of the PCR reactions consisted of one initial denaturation of 10 min at 94 ℃, then 35 cycles of denaturation of 45 s at 94 ℃, primer annealing of 45 s at 55 ℃ and extensions of 2 min at 72 ℃, and a final extension of 10 min at 72 ℃. Following the cycles, there was a final hold at 4 ℃. PCR products were detected using gel electrophoresis in a 1.5% agarose gel stained with Super GelRed (US Everbright Inc., Suzhou, China). The sequencing and phylogenetic analysis were performed following the previous descriptions for ITS.

3 RESULT 3.1 Molecular identification

The phylogenetic tree for ITS showed that all 12 samples fell into the U. prolifera - U. linza complex (Supplementary Fig.S1), and the tree for 5S rRNA spacer showed that they were clearly resolved into two clades, i.e., U. prolifera and U. linza. After that, eight samples of U. prolifera were detected by SCAR marker further, and the results showed that all four floating samples belonged to the floating ecotype of U. prolifera (Fig. 1).

3.2 Chloroplast genome of U. prolifera U161 with phylogenetic analysis

To develop interspecific genetic markers for U. prolifera and U. linza based on the intraspecific variations within U. prolifera, an attached U. prolifera strain U161 was selected for sequencing of chloroplast genome since both references of floating U. prolifera and U. linza are readily available. After genome sequencing, assembly, and annotation, it was shown that the complete chloroplast genome of U161 is 99 724 bp in size (Fig. 2) (GenBank accession No. MZ571508), encoding 95 genes including 67 protein-coding genes, 26 tRNAs, and 2 rRNAs. There are five genes (psbB, psbD, atpA, atpB, and psaB) containing one intron and there is one gene (petB) containing two introns. The overall base composition was A (37.7%), T (37.0%), C (12.6%), and G (12.7%). The voucher (assigned number MBM 287038) was deposited in the Marine Biological Museum of Chinese Academy of Sciences (MBMCAS) at the Institute of Oceanology, Chinese Academy of Sciences, China.

The ML phylogenetic tree of chloroplast genomes of Ulva was shown in Fig. 3. It was shown that the attached and floating U. prolifera gathered into a cluster which was separated from U. linza. The chloroplast genome sequence divergence was 0.3% between U. linza and the attached U. prolifera, and 0.4% between U. linza and the floating U. prolifera. This result suggested that U. linza and U. prolifera can be distinguished as two species by the whole chloroplast genome despite the intraspecific divergences within U. proifera.

3.3 Comparative analysis of chloroplast genomes among U. linza and two ecotypes of U. prolifera

As shown in Fig. 4, the relative synonymous codon usage (RSCU) values were calculated and summarized with chloroplast genomes of U. linza and two ecotypes of U. prolifera respectively. In general, it was clearly indicated that the codon selection strategies in the three chloroplast genomes were extremely similar. Except for methionine and tryptophan (RSCU=1), most amino acids were exhibited to have codon bias. A total of 26 high frequency codons (RSCU > 1), including a stop codon, were identified with A/T ending as usual in Ulva (Cai et al., 2017), while the codons with negative bias (RSCU < 1) were prone to end with G/C. The results showed that the codon usage of the three genomes are extremely conservative without potential to provide resources for interspecific discrimination. Furthermore, we analyzed the genetic variations within the non-coding regions, including introns and gene spacer regions.

The collinearity analysis was conducted with these three chloroplast genomes. It was obviously shown that none of structural rearrangements such as inversions or translocations were detected among three genomes, and the orders of similarity sequences in the chloroplast genomes of U. linza and two ecotypes of U. prolifera were almost identical except for some slight variations such as insertions and deletions mainly located in the regions of introns or gene spacers (Fig. 5). Therefore, results of both the codon bias and collinearity analysis showed that U. linza had a very close genetic relationship with U. prolifera.

To investigate the interspecific variations between U. linza and U. prolifera, these three chloroplast genomic sequences were compared using mVISTA. As shown in Fig. 6, plenty of variations were detected which were distributed in both the conserved non-coding sequences (CNS) and exon regions. By ignoring the intraspecific variations between the two ecotypes of U. prolifera, such as the psbB-psbC spacer region, only those signals which were identical between two ecotypes of U. prolifera but divergent between U. linza and U. prolifera, were further searched out to represent the interspecific variations between these two related species. A total of 454 SNPs, 131 indels and six structural variations were identified. In particular, three of the six structural variations were found to be longer than 1 000 bp. According to the position displayed on the X axis which was based on the chloroplast genome sequence of the attached ecotype of U. prolifera (MZ571508), these three regions of large structural variations were found to be located at psaB (3 kb–4 kb), petB (70.5 kb–71.8 kb), and psbB (91 kb–92 kb) respectively. Upon further analysis, each region was determined as an intron in the chloroplast genomes of U. prolifera, while it was a complete deletion in that of U. linza.

3.4 Development of new species-specific markers from chloroplast genomes

Dozens of pairs of primers were designed to target those interspecific variations between U. prolifera and U. linza which were located in either CNS or exon regions. After validation with PCR amplifications, those primers generating no products, polymorphic products, different products between two ecotypes of U. prolifera, or identical products between two related species, were all abandoned. Finally, two pairs of primers, which were designed to match the coding regions within psaB gene and petB gene respectively (Supplementary Table S2), were proved to be capable of generating species-specific signals to distinguish U. prolifera and U. linza (Fig. 7), and the sequences data have been uploaded to NCBI (Supplementary Table S3). The primers for psaB marker could amplify approximately 2 100-bp bands of the same size from both ecotypes of U. prolifera, while only about 1 000-bp bands can be amplified from U. linza. Similarly, the primers for petB marker could amplify an about 1 800-bp band from each of U. prolifera samples, whereas about 150-bp bands in U. linza samples. The ML phylogenetic trees for petB and psaB markers showed that all samples were clearly resolved into two clades, i.e., U. prolifera and U. linza, without significant genetic divergency in each clade (Supplementary Figs.S2–S3).

4 DISCUSSION

The phenotypic differentiations between the two ecotypes of U. prolifera have long been concerned, in terms of the morphology (Wang et al., 2010; Hiraoka et al., 2011; Gao et al., 2016; Ma et al., 2020), habitats (Ding et al., 2009), and transcription level of some key metabolism-related genes (He et al., 2019). Their significant differences in performances of reproductive isolation with U. linza were also described (Hiraoka et al., 2011). In particular, the genetic variations have also been revealed, by using inter-simple sequence repeat (ISSR) markers which were located throughout whole genomes (Zhao et al., 2011). A SCAR marker specific to the floating ecotype has been developed to find that this unique ecotype almost never formed a colonization population in the intertidal zone (Zhao et al., 2018). These findings implied the genetic differentiation between the two ecotypes of U. prolifera, which was confirmed to some extent by the comparative chloroplast genomic analysis in this study. In contrast, the results of all four tested molecular markers, especially 5S spacer, showed that all U. linza samples from different geographic populations were almost genetically identical, suggesting that the intraspecific genetic differences in U. linza were not significant. Therefore, in order to develop species-specific molecular markers for U. linza and U. prolifera, the influence of intraspecific differences, especially for U. prolifera that consisting of different ecotypes, should be fully considered.

In this study, the chloroplast genome from one representative attached population of U. prolifera was completely sequenced, and comparative analysis was performed with other chloroplast genomes from U. linza and the floating ecotype of U. prolifera. A strategy was proposed that only those signals of variation which were identical between two ecotypes of U. prolifera but divergent between U. linza and U. prolifera, were selected to develop the interspecific markers for U. linza and U. prolifera. Two candidate markers, i.e. psaB and petB, were validated to be capable of distinguishing these two related species. These new markers are expected to be used in surveys for Ulva species composition and green tide monitoring especially in the southern Yellow Sea. This sea area has experienced severe green tides for more than a decade (Yu et al., 2018). It was proposed that, the fouling green seaweeds on the nori rafts at Subei, in which both U. linza and U. prolifera were major members (Fan et al., 2015; Huo et al., 2015), provided the origin of biomass for the green tides in the Yellow Sea (Liu et al., 2009), and only the floating ecotype of U. prolifera finally succeed to be extremely dominant (Zhao et al., 2015). In addtion, it was suggested that the Ulva species composition and biomass in samples, including those fouling green seaweeds and the Ulva micro-propagules distributed in seawaters or surface sediments in this area, might contribute greatly to the interannual characteristics of the green tides (Song et al., 2015). Therefore, the new interspecific markers developed in this study, in combination with the existing floating ecotypes-specific marker, are expected to be able to characterise the detailed dynamic characteristics of the Yellow Sea green tide, and provided important data for effective risk mornitoring and management.

Organelle genomes contain abundant genetic resources, the mitochondrial genome sizes in Ulva vary between 55 kb to 88 kb, and the chloroplast genome sizes are 86 kb–119 kb. At present, genome sequences including 33 mitochondria and 26 chloroplasts from Ulva have been available in the GenBank database (https://www.ncbi.nlm.nih.gov), which were conducive to the development of molecular markers and used for inter- or intra-specific phylogenetic analysis. Recent studies showed that, for some widely distributed Ulva species, such as U. pertusa (synonym of U. australis) and U. compressa, a certain degree of intraspecific variations in organelle genomes have been detected among different geographic populations (Liu et al., 2017, 2020a; Cai et al., 2021). Since U. linza and U. prolifera also occurred worldwide, the organelle genome resources we provided in this study could contribute to validation or development of interspecific markers in future.

Moreover, in addition to the interspecific markers for U. linza and U. prolifera, the candidates of intraspecific markers specific to the floating ecotype of U. prolifera were also noted in this study. Novel organelle genome-derived markers could be developed from the chloroplast genomes and could be used together with the nuclear genome-derived SCAR marker for the ecological investigation of the green tide of Yellow Sea (Zhao et al., 2015).

5 CONCLUSION

In this study, a chloroplast genome from one attached population of U. prolifera was completely sequenced, and comparative genome analysis was performed with other existing chloroplast genomes from U. linza and the floating ecotype of U. prolifera. The results showed that in spite of the high level of collinearity among three genomes, there were plenty of interspecific and intraspecific genetic variations. Two developed markers, psaB and petB, were shown to be able to distinguish these two closely related species and were applicable to more attached populations of U. prolifera from a wide range of geographical sources.

6 DATA AVAILABILITY STATEMENT

The genome sequence data of this study are openly available in GenBank of NCBI at https://www.ncbi.nlm.nih.gov under the accession No. MZ571508. The datasets analyzed during the current study were available from the corresponding author on reasonable request.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S3 and Figs.S1–S3) is available in the online version of this article at https://doi.org/10.1007/s00343-022-2045-x.