A laboratory-cultured genetically pure line of P. yezoensis RZ-58 was used in the experiments. The original thallus was collected from the nori farm in Putian, Fujian Province, China. At the laboratory, individual somatic cells were enzymatically isolated from the thallus and the allelically homozygous sporophytes (conchocelis) were obtained after spontaneous haploid doubling. The gametophytes from the homozygous sporophytes are homogenous and haploid. The gametophytes were then cultured in bubbling sterilized seawater with Provasoli’s enrichment solution (PES) medium at 8 ± 1°C with a light concentration of 50 μmol photons m^–2s^–1 and a 12:12 light: dark (L: D) photoperiod. The wild P. yezoensis were collected from three geographical populations: Qingdao (N36°03′29″, E120°18′56″), Yantai (N37°49′35″, E120°44′40″), and Weihai (N37°27′53″, E122°08′27″) of the Shandong Province, China, during March 2016. The samples were randomly collected from the intertidal zone at intervals of at least 5 m to prevent the collection of asexual reproductive ramets. A total of 53 haploid thalli with morphological diversity were selected, 19 from Qingdao, 19 from Yantai, and 15 from Weihai, for the genetic analysis of the different groups. The collected samples were thoroughly cleaned with sterile seawater and frozen at −20°C.

Total DNA was extracted from the frozen thallus material in accordance with previously reported methods (Wattier et al., 2000). The SMRT sequencing library was constructed using the PacBio^® Template Prep Kit with an insert size of 20 kb according to the manufacturer’s instructions. The genome was sequenced using the PacBio RS II platform (Pacific Biosciences, Menlo Park, CA, United States). The PacBio raw reads (polymerase reads) were filtered using the following parameters: removal of the adaptor sequences, minimum polymerase read quality = 0.80, and minimum polymerase read length = 1,000 bp. BLASR was used to extract the organellar sequences using the reference genome (NC_017837 and KC517072) (Chaisson and Tesler, 2012). The extracted sequences were assembled using the RS_HGAP_Assembly 0.3 protocol, and Quiver was run to polish the accurate consensus in the SMRT Analysis v2.3.0 using default parameters (Chin et al., 2013). The circular genomes were then validated by mapping the Illumina HiSeq data and the PCR-based sequencing. We constructed 500-bp paired-end libraries according to the instructions of the NEBNext^® Ultra^TM II DNA Library Prep Kit (NEB, Beijing, China), and the libraries were sequenced on an Illumina Genome Analyzer. The adapters and low-quality reads were removed using the NGS QC Toolkit (Patel and Jain, 2012). Organellar genomes were assembled using the NOVOPlasty software, which utilized a seed-and-extend algorithm that starts from a related or distant single seed sequence (Dierckxsens et al., 2017). We selected the rbcL and atp8 genes as the seeds to assemble the plastid DNA (ptDNA) and mitochondrial DNA (mtDNA), respectively. The genomes were annotated using the online tool ORFfinder² and aligned using BLASTX and BLASTN searches at the National Center for Biotechnology Information (NCBI) website³. Multiple genomes were aligned and visualized using the mVISTA tool (Mayor et al., 2000). For the molecular identification, the rbcL genes from the wild collections of P. yezoensis and the 82 species of Pyropia (Supplementary Table S1) were downloaded from GenBank and used to construct the phylogenetic tree. The sequences were aligned using the MAFFT version 5 program (Katoh et al., 2005) and trimmed using trimAl with the option “automated1” (Capella-Gutiérrez et al., 2009). The neighbor-joining (NJ) phylogenetic tree was created using MEGA7 with a 1,000-bootstrap replicate (Kumar et al., 2016).

After the quality controls were applied, the clean reads were mapped to the reference genome using BWA (Li and Durbin, 2009). The reference genomes were indexed, and the command “aln” was used to find the suffix array coordinates of good matches for each read. The best alignments were then converted into the SAM format using the “sampe” command. The alignment results were further corrected using the Picard package with the following two commands: “AddOrReplaceReadGroups” was used to replace all read groups in the INPUT file with a new read group and to assign all of the reads to this read group in the OUTPUT BAM; and “FixMateInformation” was used to ensure that all of the mate-pair information was synchronized between each read and its mate pair⁴. Potential PCR duplicates were removed using the SAMtools command “rmdup” (Li et al., 2009).

After alignment, we performed SNP calling on a population scale using a Bayesian approach based on the SAMtools program. The command “mpileup” was used with the parameters “- q 1 - C 50 - S - D - m 2 - F 0.002 - u.” The subsequent analysis filtered the raw SNPs using VCFtools (Danecek et al., 2011) and GATK (McKenna et al., 2010) with the following criteria: (1) coverage depth ≥8; (2) mapping quality ≥30; (3) the distance of adjacent SNPs ≥5 bp; (4) only SNPs that occurred in more than 50% of the individuals were retained; and (5) all SNPs with a minor allele frequency (MAF) ≥0.1 were removed. We then used the high-quality SNPs to infer genetic analysis.

Nucleotide diversity (pi), Tajima’s D, transition/transversion ratio (Ts/Tv) and population divergence (F_ST value) were calculated using VCFtools with a sliding window of 1,000 bp in length and a 500-bp step size (Danecek et al., 2011). Tajima’s D test statistic was calculated using DnaSP based on the differences between the number of segregating sites and the average number of nucleotide differences (Librado and Rozas, 2009). The phylogenetic tree was inferred using TASSEL under the p-distances model using SNPs in a population scale (Bradbury et al., 2007). Principal component analysis (PCA) was performed using the EIGENSOFT software (Price et al., 2006). The organelle genetic groups was further verified with structural classes identified from population structure analysis algorithms using Structure 2.3.4 (Pritchard et al., 2000), which used a Bayesian approach to identify subpopulations with distinct allelic frequencies and places the individuals into K clusters. The number of populations (K) was set from 2 to 7 and repeated 10 times with a burn-in period of 100,000 iterations and 200,000 Markov chain Monte Carlo (MCMC) repeats. The K optimum was evaluated using Structure Harvester (Earl and Vonholdt, 2012).

To assemble the high-quality reference organellar genomes, 712 Mb of PacBio^® subread bases (a single pass of the template without the adapter sequences) and a total of 92,569 reads were generated with a mean read length of 7,690 bp and an N50 size of 9,898 bp. The finished mitochondrial genome of P. yezoensis is 41,692 bp in size with an overall GC content of 32.72% (Figure 1A). The plastid genome is 191,977 bp with a GC content of 33.09% (Figure 1B). We confirmed that all of the differences identified in our assembly were strongly supported by mapping with the highly accurate Illumina short reads.

We then resequenced 53 accessions of wild P. yezoensis that were collected from Qingdao, Weihai and the Yantai coast of China. After removing low quality reads and enrichment, we obtained ∼4 Gb of organelle sequence bases with coverage that ranged from 260× to 370×. For the molecular identification, a phylogenetic tree that used the plastid rbcL gene showed that all the samples formed a well-supported clade together with the previously identified P. yezoensis (Supplementary Figure S1). The plastid genomes of P. yezoensis were similar in length (191,885 to 192,002 bp) and GC content (33.07 to 33.09%) and had highly conserved gene synteny (Figure 2B and Supplementary Table S2). Sequence differences mainly resulted from tandem repeat polymorphisms that were located in intergenic regions. The repeat motifs ranged in length from 1 to 20 base pairs and the repetitions ranged from 1 to 2. For example, there was a (TTTA)_n SSR polymorphism between the chlN and ycf46 genes. A (TAGTTAGTAGTTTAAATAGA)_n minisatellite polymorphism was also found between the ORF287 and ycf24 genes.

The mitochondrial genome lengths within the populations of P. yezoensis were significantly different, ranging from 34,031 to 40,758 bp (Figure 2A). The GC content ranged from 32.45 to 32.79% (Supplementary Table S2). The observed differences were mainly caused by different numbers of mitochondrial group II introns in the large subunit ribosomal RNA (rnl) gene and the cytochrome c oxidase subunit 1 (cox1) gene. Recent studies have shown that in the Pyropia species, the rnl gene is composed of four coding regions and the cox1 gene is composed of five coding regions (Hwang et al., 2013). According to the presence of the intron, the 53 individuals displayed five structural patterns of the rnl gene and three of the cox1 gene (Figure 3). The combined structural variation of the rnl and cox1 genes exhibited six genotypes in total (one genotype of RZ-58 and five genotypes of the wild collections) among the whole collections (Supplementary Table S2). Only one individual sample, QD23 (Qingdao23), had an intron between each exon of the rnl gene. The others had one or two absent introns (Figure 3A). All 53 individuals and the reference RZ-58 lacked an intron between exon 1 and exon 2 in the cox1 gene (Figure 3B), which was present in the mtDNA of Pyropia haitanensis (NC_017751) and Pyropia tenera (NC_021475). Moreover, each intron contained one full-length intronic ORF. These ORFs encode reverse transcriptases that function in intron mobility and as maturases in RNA splicing. Furthermore, according to the different intronic structural patterns in the population, we annotated the exon–intron boundaries of the rnl and cox1 genes from P. yezoensis mtDNA. The splice junctions of introns 1 and 2 of the rnl gene are a GA/AT and CG/AA motif, respectively (here, we defined the junction start/end as the intronic start/end bases). The splice junctions of introns 2, 3, and 4 of the cox1 gene are a CA/CC, GT/AT, and GT/AC motif, respectively, from the 5’ end to the 3’ end (Supplementary Figure S2).

To further investigate whether there was organellar heteroplasmy within individual P. yezoensis samples, we aligned the NGS data of each individual to their corresponding assembled genomes. The least supported reads of a genotype at each locus was set to four. The aligned results indicated that there was no supporting evidence from these data for heteroplasmy in mtDNA or ptDNA, although four heterozygous loci were found in the ptDNA of all individuals. These two loci were in two direct non-identical repeat regions, which could lead to the misalignment. Validation was then performed using intensive PCR-based sequencing of the repeat regions, which indicated that there were no heterozygous loci. Therefore, we concluded that the P. yezoensis organellar genomes do not display heteroplasmy.

Using the mtDNA and ptDNA of P. yezoensis RZ-58 as references, we identified 849 variants, including 829 SNPs and 20 indels across the whole collection (Table 1). There were 463 SNPs and 10 indels detected in the mtDNA with a variation density of 11.35 per 1 kb. After the location of variation across the genome was determined, 364 variants were found within genic regions (Table 1). The different allele types indicated that C/T and A/G were the major forms of base changes with an overall Ts/Tv (transition/transversion) ratio of 4.08. There were 366 SNPs and 10 indels detected in the ptDNA with a variation density of 1.96 per 1 kb. Variation in the plastid genome consisted of 282 variants in genic regions and 94 in intergenic spacers (Table 1). Similarly, C/T and A/G were again the major allele types with an overall Ts/Tv ratio of 2.56. The Qingdao populations held the highest number of SNPs in their mtDNA and the least number of SNPs in their ptDNA (Table 1).

The nucleotide diversity (pi) of the wild Chinese P. yezoensis was calculated at a mean pi of 0.00361 (mtDNA) and 0.000303 (ptDNA) for the whole collection (Table 2). These results indicated that the mutation rate of the mitochondrial genome is greater than that of the plastid genome. A similar trend has been observed in other species of red algae (Smith et al., 2012; Xu et al., 2018b). Among the subgroups, the wild Qingdao P. yezoensis had the highest diversity (0.00348 for mtDNA and 0.000388 for ptDNA) while Yantai had the lowest (0.00222 for mtDNA and 0.000291 for ptDNA).

The F_ST values across the three populations were estimated using the organellar SNP data (Supplementary Figure S3). The F_ST for all pairs of the populations were different and ranged from 0.016 to 0.272 for the mtDNA and 0.047 to 0.087 for the ptDNA. The F_ST was largest between the wild Qingdao and Yantai populations with average values of 0.272 for mtDNA and 0.087 for ptDNA, which indicates a high degree of genetic isolation between these two groups. The wild Weihai and Yantai populations had the lowest F_ST values (0.0128 for mtDNA and 0.0436 for ptDNA), which could be attributed to a common ancestry or admixture between the populations.

The Tajima’s D value of the organellar genome was also examined to detect balancing selection. The results identified a location fluctuation in 1 kb bins with positive, negative, and zero values across the genome. The Tajima’s D value for the mtDNA indicated that the wild Qingdao P. yezoensis had more positive positions, which indicated a decrease in its population size and/or balancing selection. In contrast, the Yantai groups had more negative positions, which indicated an expansion in population size and/or purifying selection (Supplementary Figure S3). The Tajima’s D value of the ptDNA also shows more positive positions in the Qingdao group and more negative positions in the Yantai group (Supplementary Figure S3). The closer mean of the Tajima’s D value to 0 indicated weak selection acting on the Weihai group (Table 2). However, Tajima’s D test of each population or whole collection was not found to be significant (p > 0.05).

To address the different evolutionary rates of mtDNA and ptDNA, we combined the mitochondrial and plastid SNPs to conduct the cluster analysis. A total of three sub-groups were identified from the phylogenetic tree (Figure 4A). The dendrogram showed that individuals were mixed among the different geographical populations without obvious regional features. The PCA used the whole SNP dataset and was conducted in EIGENSOFT (Figure 4B). The three main groups were inferred from the first two principal components (PCs), which explained more than 51% of the variation. The PCA plot produced similar results, confirming reliability of the clustering results. We also found that the whole collection clustered according to the different haplotypes of intron structures. Cluster 1 was composed of the R011C0100 haplotype (R = rnl, C = cox1, 0 = intron absent, 1 = intron present). Cluster 2 was composed of the R110C0101 haplotype and cluster 3 was composed of the R110C0101, R111C0100, and R010C0100 haplotypes. Individuals from each haplotype group were present within the same sampling areas at sympatric localities, indicating a genetic break of wild P. yezoensis in the absence of geographic barriers. A Bayesian method based on the MCMC algorithm implemented in STRUCTURE was also performed for identification of genetically different groups. We identified the major substructure groups using an optimal K-value of 3 that was based on the highest ΔK (Supplementary Figure S4), corresponding to the groups identified in wild P. yezoensis. This analysis also revealed there was no significant regional geographic structure in the genetic variation within the groups (Supplementary Figure S5).

Variation in the size of P. yezoensis mtDNA was largely due to the different number of mitochondrial group II introns in the rnl and cox1 genes. The differing number of introns was previously reported on in red algae mtDNA at the interspecific level (Hwang et al., 2013; Xu et al., 2018a). In this study, we detected intraspecific variation in intron structure across the 53 Chinese P. yezoensis strains. The combined structural variation of the rnl and cox1 genes resulted in six haplotypes. Recently, Hwang et al. (2018) reported on the genetic features of gene structure in the mitochondrial rnl and cox1 genes of P. yezoensis and found 12 haplotypes from 27 strains from Korea and only one haplotype from 13 strains from Japan. In theory, the combination of the rnl and cox1 genes of P. yezoensis could produce 128 haplotypes. The high intraspecific genetic variation in exon and intron structures could be the more appropriate molecular markers for the identification of species or subspecies. The mitochondrial cox1 gene is considered to be the most suitable DNA barcode for animal and algal species due to its remarkable ability to identify species boundaries (Hebert et al., 2003; Saunders, 2005; Lavinia et al., 2006; Mcdevit and Saunders, 2010). However, the cox1 region is difficult to amplify and sequence compared to other commonly used barcodes for red algae (Sherwood et al., 2010; Du et al., 2015). Results indicate that the failed amplification of cox1 could be caused by the presence or absence of introns. Therefore, we recommend the use of primers that are designed for each exonic region to amplify the cox1 gene.

Heteroplasmy is the presence of DNA from more than one type of organelle in an organism and has been reported in many animals and plants but not yet in this group of red algae (Yamato and Newton, 1999; Chinnery et al., 2000; Steel et al., 2000; Coyer et al., 2004; Frey et al., 2005). However, some of the earlier reports of heteroplasmy may be due to organelle–nuclear DNA transfer (Martin, 2003; Cullis et al., 2009). For instance, the application of NGS followed by read mapping analysis provided no evidence of chloroplast heteroplasmy in sugarcane (Hoang et al., 2015). When assessing heteroplasmy, it is necessary to exclude transferred DNA. In this study, we compared the organelle and nuclear genomes of P. yezoensis in advance, and no transferred DNA fragments could be found (data were not shown). We then re-aligned the NGS data to each of the organellar genomes and found no evidence of heteroplasmy in the mtDNA or the ptDNA. Heteroplasmy could arise from spontaneous mutations or biparental inheritance (Chat et al., 2002; Kmiec et al., 2006). However, a previous study suggested that organelles of P. yezoensis are inherited uniparentally from the female parent (Choi et al., 2010). Moreover, organelle mutations are expected to be purified by genetic drift toward homoplasmy if the mutation is harmful or the two genotypes differ in their replication speed (Hill et al., 2014; Stephan et al., 2015). Therefore, it is plausible that maternal inheritance, an absence of organelle–nuclear DNA transfer and genetic drift could contribute to the homoplasmy observed in the P. yezoensis organellar genomes.

Currently, the complete nuclear genome sequence of Pyropia is unavailable, which is a major constraint for population genetics research for this group. Therefore, information on the polymorphic DNA in organellar genomes can be a useful tool for genetic and evolutionary studies due to the non-recombinant, uniparentally inherited, and effectively haploid genomes (Zeng et al., 2012; Preston et al., 2014). Traditional methods of DNA restriction mapping that use hybridization-based RFLP or SSR have several drawbacks, including low resolution, high labor intensity, and a requirement for large amounts of isolated DNA (Byrne et al., 1993; Wu et al., 2008; Yotsukura et al., 2016). Some of the amplified polymorphic bands were from microorganismal DNA rather than the host DNA due to the number of bacteria or microalgae attached to the surface of the algae. The non-target DNAs from other attached organisms are potentially confounding agents for AFLP DNA fingerprinting. When a reference genome is available, the contaminated DNA can be identified and avoided. Some investigators have predicted that whole-genome analysis will replace other molecular methods such as the AFLP technique. This study represents the first identification of organellar genome SNPs for the analysis of genetic diversity of the Chinese wild P. yezoensis. Compared with existing techniques, this method requires a minimal amount of DNA without any non-target DNA and avoids the laborious experimental procedures that are associated with chloroplast DNA (cpDNA) and mtDNA purification. The sequencing reads of the organellar DNA accounted for 50–90% of the total reads. Considering the small size of the organellar genomes, we could pool the sequencing libraries by adding unique indices. As a result, 200 Mb of sequence bases per individual is enough to achieve a read depth of 400 × for the organellar genomes. Sequencing with higher genomic coverage could decrease the number of missing values and thus obtain even more accurate results. In the present study, we identified a considerable number of SNP markers for the commercially important seaweed P. yezoensis, which led to the identification of 829 SNPs and 20 indels across the diverse collection of specimens. The number of polymorphisms is much larger than the SSRs that were used in previous studies (Kong et al., 2009; Zhang et al., 2009; Hu et al., 2017). Our results demonstrate the high throughput, accurate discovery of organellar DNA polymorphisms, and the low costs associated with this sequencing strategy.

A broad range of molecular markers is an important genetic resource for breeding applications, which include germplasm conservation, hybridization breeding, quantitative trait loci identification, association mapping, and genomic domestication. All of the above efforts require prior information on genetic diversity (Rao and Hodgkin, 2002; Zhu et al., 2008; Rauf et al., 2010). In this study, applying SNP markers to the whole collection and the geographic populations identified relatively high levels of genetic diversity. A P. yezoensis population that was sampled from five regions along the coast of Shandong Province also showed high levels of genetic diversity as measured by AFLP (Cao et al., 2018). Yang et al. (2003) analyzed genetic variation in 11 lines of P. yezoensis using AFLP and concluded that samples from China shared these features of high genetic diversity. Twelve SSR loci were also used to analyze the genetic diversity of P. yezoensis using wild populations from Qingdao, Yantai, and Weihai (Wei, 2012). These previous results indicated the wild Yantai P. yezoensis had the highest genetic diversity, which conflicts with our finding that the Qingdao population exhibited the highest diversity. The differences in these reported results could be explained by the use of different molecular markers but could also be attributed to decreasing diversity caused by the expansion of cultivated fields in the Yantai coast area during recent years, which overlaps with the habitat of wild P. yezoensis. The genetic diversity of plants can be affected by both genetic and environmental factors (Nevo, 2001). While wild P. yezoensis maintains a relatively high level of genetic diversity, its habitat on natural rocks has been gradually destroyed by coastal pollution and human activities. In order to implement long-term breeding strategies, we first need to make efforts to protect existing wild resources, establish core germplasm banks, and update our breeding and cultivation techniques.

A well-supported organellar DNA break was observed in different geographic populations. Some land plants also display high levels of admixture within their analyzed populations (Aranzana et al., 2010; Barać et al., 2017), which could be attributed to a complex breeding history. We speculate that the presence of a genetic organellar DNA break in the absence of a geographic barrier is the result of historic biogeographic events. Previous studies showed that divergent mtDNA emerges in the absence of structured geographic barriers (Miller et al., 2006). The organellar DNA break might be caused by the sundering of a historically continuous distribution of P. yezoensis due to climate change and the retreat to the coast of Shandong Province, which was facilitated by the winter–spring sea currents in the East China sea. Another contributing factor could be the attractiveness of wild P. yezoensis as a source of food for humans and animals, which led to the dispersal of conchospores or monospores to distinct areas. This phenomenon is often observed for some fruit crops (Mariette et al., 2010; Cao et al., 2012). The genetic break of organellar DNA could lead to the inference that a cryptic species complex within P. yezoensis was reproductively isolated, which was previously reported for the Porphyra/Pyropia complex based on the phylogenetic analysis (Niwa et al., 2009). However, if the observed organellar DNA break reflects reproductive isolation between the haplotype groups, then the break pattern should also be reflected in the nuclear genome of the contact-zone individuals. Organellar and nuclear DNA can differ in their relative responses to historic and recent evolutionary processes due to their distinct modes of inheritance. Discordance between cytoplasmic and nuclear DNA patterns has been previously observed in the wild populations of some species (Egger et al., 2007; Yang and Kenagy, 2009). Therefore, genetic information from nuclear DNA should be further investigated to clarify the relative likelihoods of recent gene flow or the retention of ancestral polymorphisms.