A collection of 13 Saccharina cultivars was evaluated in the present study. Six cultivars came from Changdao, and 7 cultivars came from Rongcheng ( Figure 1 ). These cultivars were generated from different breeding methods. Detailed information on the regional distribution of the 13 cultivars is provided in Table 1 .

Approximately 3- to 5-cm young healthy leaves were collected, frozen in liquid nitrogen, and used for genomic DNA isolation. Each sample was ground into a powder using a mortar, and genomic DNA was isolated and purified via the Plant Genomic DNA Kit (TIANGEN BIOTECH CO., LTD). The purity and concentration of extracted DNA were qualified using a spectrophotometer (NanoDrop-2000); 117 qualified DNA samples were obtained, and 13 unqualified DNA samples were not used. Qualified DNA samples were diluted to 50 ng/μl. The high-throughput sequencing method is individual sequencing.

Genomic DNA extracted from sporophytes were genotyped using the reported SLAF (Specific-Locus Amplified Fragment) method (Sun et al., 2013). Restriction enzyme combinations were tested and selected following the four criteria: (I) low proportion of restriction fragments comprising repeat sequences; (II) even distribution of restriction fragments through chromosomes; (III) simulated fragments match uniquely to the reference genome; and (IV) ample number of SLAF tags. RsaI and HaeIII were selected as the fittest restriction endonucleases to obtain evenly distributed SLAF tags and to avoid repetitive SLAF tags.

To obtain the SLAF tags, genomic DNA of each accession was digested with RsaI and HaeIII. Then, the target fragments (fragment sequence of 314–394 bp) for SLAF library construction were selected for 3’ end reparation, dual-index paired-end adapter ligation, and PCR amplification. The selected SLAFs were sequenced with 126 bp × 2 pairs by Illumina high-throughput sequencing platform (Illumina, Inc., USA). After sequencing, the Q30 and GC content were checked for sequencing quality.

The files of SLAF reads were aligned to the reference genome of S. japonica (https://www.ncbi.nlm.nih.gov/genome/?term=txid88149[orgn]) using the software BWA (Li and Durbin, 2009). Oryza sativa ssp. japonica DNA (http://rapdb.dna.affrc.go.jp/) was used as a control to assess the normal rate of enzyme digestion. The two sequenced reads were compared through SOAP (Arabidopsis Genome Initiative, 2000) software. SNPs were developed based on the polymorphic SLAF tags information via the software of GATK (McKenna et al., 2010) and SAMtools (Li H. et al., 2009). SNPs predicted by both GATK and SNPs were screened with the criteria of minor allele frequency (MAF) > 0.05 and integrity > 0.5.

Based on the SNPs identified in this study, genetic distances were calculated between the cultivars using the neighbor-joining method. Phylogenetic tree was constructed using the MEGA5 software (Tamura et al., 2011). The K-means clustering analysis using admixture software (Alexander et al., 2009) was conducted assuming that the K-value ranged from 2 to 20. A principal component analysis (PCA) was performed using EIGENSOFT (Price et al., 2006).

Prior to population genetic analysis, VCF tools (Version 0.1.14) were used to rigorously obtain the most informative SNPs (Eimanifar et al., 2018). The nucleic acid diversity (Pi), heterozygosity (Ho), and expected heterozygosity (He) of SNP loci in each population were separately calculated using PLINK2 (Purcell et al., 2007).

The efficiency of two terminal comparison is 92.27%, the efficiency of enzyme digestion is 87.67%, and the read length inserts are distributed within the expected range. The efficiency of database building comparison is normal, SLAF library building is normal, and sequencing quality is normal.

In total, 277.27 M reads were obtained from the raw data of the SLAF library. The Q30 ratio was 87.84%, and the GC content was 48.49% ( Table S1 ). A total of 178,047 SLAFs were anchored to the whole Saccharina genome, including 124,745 polymorphic SLAFs (approximately 70.1%), and the average sequencing depth per sample is 12.98× ( Table S2 ). Using SAMtools and GATK for realignment and detection, a total of 2,847,332 SNPs from the original reads were developed ( Table S3 ). After the low-depth SNPs were filtered out, 340,398 SNPs were successfully encoded.

The phylogenetic relationships between the 117 entries from 13 cultivars, based on the set of 340,398 SNPs, are illustrated in Figure 2 . The cluster results showed that DF1 and DF4 are pure and independent cultivars. DF2 were almost all dispersed between the DF6 and DF8 groups, YZ were almost all dispersed between 407 groups, and DF3 were almost all dispersed between the BN and DBC groups, suggesting genetic permeation between these cultivars. Most of DF5 are in one group, with only a few clustering into the DF1 group. DBC and HK dispersed widely.

The first three PCs accounted for 10.07% of the population-wide SNP variation; the first, second, and the third principal components explained 4.34%, 3.14%, and 2.59% of the genetic diversity, respectively ( Figure 3 ), and the hybrid level among the cultivars was high.

The PCA results confirmed that “HK” was widely dispersed from the YZ and XBN, providing evidence of a large diverse genetic background. DF4 is a pure and independent cultivar. DF2, DF6, and DF8 are mixed groups. DF3, BN, and DBC are also mixed. DF1 and DF5 in relationship are near with each other.

The estimated membership fractions of the 13 cultivars for different values of K ranged from 2 to 20, and the maximum likelihood revealed by the population structure showed an optimum value of 2 (K = 2) ( Figures 4 , 5 ), which indicated that the entire population could be categorized into 2 groups ( Table S4 ). Group 1 contained 3 cultivars, DF2, DF6 and DF8. Other cultivars belong to group 2.

Prior to population genetic analysis, we performed secondary filtering of the SNP data following the three criteria: (I) SNP loci that could be classified in <80% of individuals; (II) SNP loci with minor allele frequency (MAF) <0.05; and (III) SNP loci with an allele gene greater than 2.

The Pi distribution of 407, DF2, and DF8 was similar with average Pi values of 4.24×10^-5, 4.27×10^-5, and 4.28×10^-5 respectively ( Figure 6 ), which are lower than the average (the average Pi = 4.38×10^-5). The highest Pi is 5.10×10^-5in DF5.

We found from Table 2 and Figure 7 that most observed heterozygosity (Ho) was lower than expected heterozygosity (He) at the population level and that the inbreeding coefficient Fis was high. DF5 has the highest Ho, indicating that the population heterozygosity was high. DF1 and DF3 have the lowest Ho. XBN and DF5 have the highest F.

In this study, the concordance of the model-based structure analysis reveals two groups in the population with the phylogenetic tree and PCA clustering. However, there is no correlation between molecular classification and main character classification ( Table 3 ). DF2 and DF6 belong to one group according to molecular analysis. DF2, DF3, and DF6 belong to one group according to the characteristics of high light resistance. DF1, DF6, and YZ belong to one group according to the characteristics of temperature resistance. HK, XBN, DF2, and DF3 belong to one group according to the characteristics of high yield. Because of the mixing in the actual farming process, we cannot accurately evaluate the phenotypic characteristics of 407, DF4, DF5, and DF8. XBN is a cultivar mainly cultured in Rongcheng, and the phenotypic characteristics can be obtained through personal investigation.

Generally, the diversity of phenotypic characteristics of Saccharina cultivars is high, but some characteristics are likely hereditable, and the correlated genes have changed. However, some traits are more affected by circumstance factor, and genetic models were more complex too. Thus, the reasonable classification of cultivated strains is worth further study.

The parents of most Chinese Saccharina cultivars are from S. japonica; long-term selection results in genetic variability. Across the 13 Saccharina cultivars examined in this study, we observed that the highest He and Ho value is 0.295, which is lower than natural kelp population in Japan and Russia (Zhang et al., 2017). However, another study showed that the genetic diversity of cultivars group was higher than the natural kelp populations (Li et al., 2017). The AFLP analysis of cultured and wild population of kelp in Dalian coast showed that the level of genetic variation of two kinds of population has no significant difference (Wang, 2009).

We can discuss the reasons from two aspects: one is that the identification of germplasm of Saccharina cultivars still relies on the experience of technical staff, to increase the impact of the intermixed germplasm; the other is that the analysis of genetic diversity detected by different molecular markers may show different results (Zhang and Yang, 2008). For the kelp analysis of genetic diversity, different analysis methods should be applied and compared; meanwhile, the phenotypic characteristics should be an important reference. Besides, for the cultivars, strong genetic coherence is beneficial for application, but is likely harmful for continuous expression of some excellent characteristics, and it is worthy to find a balance between the two sides in a further study. The reason for high Fis is most probably that the direct ancestor of the cultivation cultivars of kelp was the S. japonica of Japan in China (Shan et al., 2010). It is also necessary to continuously detect genetic structure changes in combination with apparent traits.