Identification of Main Oyster Species and Comparison of Their Genetic Diversity in Zhejiang Coast, South of Yangtze River Estuary

Oysters are an important aquaculture species distributed worldwide, including in Zhejiang Province, located on the east coast of China. Because of the high diversity and complicated introduction history of oysters and their seedlings, there has been much disagreement regarding the origin of each species, and the dominant and indigenous species remain unclear. We sampled 16 batches of oysters from seven sites in three aquaculture bays and found two main oyster species, Crassostrea sikamea and Crassostrea angulata. The former occupied the higher intertidal zone and comprised more than 70% of the cultured oysters. Based on the cytochrome oxidase C subunit I (COI) and mitochondrial noncoding region (MNR), C. sikamea showed higher genetic diversity than C. angulata. The analysis of molecular variance among COI sequences of these species from the Xiangshan Bay populations were comparable to those of other populations and showed that most of the molecular variance was within groups, which was consistent with the low pairwise fixation index FST values. The neutrality test revealed that C. sikamea experienced population expansion events, whereas for C. angulata, the significant Fu’s Fs and non-significant Tajima’s D test results may indicate a possible population expansion event, implying that C. sikamea is likely an indigenous species. The method established based on internal transcribed spacer 1 digestion by the HindIII restriction enzyme is useful for identifying C. sikamea and C. angulata in the local region. The specific primers on the MNR sequence show potential for distinguishing C. sikamea from four other important Crassostrea oysters. These results highlight the abundance of C. sikamea on the Zhejiang coast and lay the foundation for protecting and utilizing the local oyster germplasm resources and for the sustainable development of the oyster industry.


INTRODUCTION
Oysters are important marine aquaculture shellfish that are distributed globally, with an annual production of more than 6 million tons (Botta et al., 2020). There are abundant oyster resources in China, and the dominant species in the northern Liaoning and Shandong Provinces is the Pacific oyster Crassostrea gigas ; in Fujian Province, the Portuguese oyster (C. angulata) (Wang et al., 2010); and in Guangdong and Guangxi Provinces, the Hong Kong oyster (C. hongkongensis) (Lam and Morton, 2003) and Suminoe oyster (C. ariakensis) (Wang et al., 2004). Zhejiang Province is located on the east coast of China, south of the Yangtze River estuary. According to a newly unearthed shell mound, oysters have been utilized for at least 8000 years in the local area (Kaihao, 2020). Important aquaculture bays include Xiangshan Bay, Sanmen Bay, and Yueqing Bay along the coast. The local aquaculture history of the oyster is approximately 800 years old, dating back to the Song dynasty (A.D. 1225). Furthermore, the "Xidian oyster, " cultured in the Xidian region, was listed as a tribute to the royal family of the Qing Dynasty and enjoyed a high reputation. As a highly important farmed shellfish in Zhejiang Province, oyster production reached 223,000 tons in 2019 (Bureau of Fisheries, 2020).
Traditionally, Zhejiang's local oysters are obviously different from the above-mentioned dominant species in terms of shell shape, low growth rate, and small body size, and they have been called monk-hat oysters. Previous research found that the most common small oyster in Zhejiang is C. angulata (Wang et al., 2010) which is thought to be the main oyster cultured in Zhejiang and Fujian Provinces (Li et al., 2015). Kumamoto oysters (Wang et al., 2013) and Suminoe oysters (Wang et al., 2004) are also distributed in this area, although species abundance has not been reported. In addition, the Pacific oyster C. gigas was introduced approximately 20 years ago, in the 1980s, from southern Japan and the spats were thought to have settled, whereas this oyster was gradually abandoned for various reasons (Jin, 2000). Some studies have reported that hybridization occurs naturally between C. gigas and C. angulata (Huvet et al., 2004). Moreover, for approximately 30 years, local oyster farmers generally purchased seedlings from different locations in Fujian Province, of which the majority are C. angulata, and the introduced spat could obviously improve oyster industry yields. These conditions have increased the complexity of oyster species along the coast of Zhejiang. The exchange of seedlings may impact the genetic structure of indigenous populations and reduce genetic diversity (Guo, 2009).
Mitochondrial DNA sequences, such as cytochrome oxidase C subunit I (COI), are preferred genetic markers because they are highly conserved protein-coding genes in the mitochondrial genome of animals (Folmer et al., 1994); thus, they have been widely used for species identification and genetic diversity assessment of mollusks (Hsiao et al., 2016;Sekino et al., 2016;In et al., 2017;Nowland et al., 2019;Özcan Gökçek et al., 2020;Tan et al., 2020;Melo et al., 2021) and other invertebrates (Yue et al., 2020). Moreover, a mitochondrial noncoding region (MNR) lying between the Gly and Val tRNA synthetase regions has the potential to distinguish different oyster populations (Aranishi and Okimoto, 2005) and has recently been used to analyze the population genetics of C. gigas and C. angulata in Europe and Asia (Moehler et al., 2011;Grade et al., 2016).
Owing to high similarities in the appearance and sympatric distribution of up to four oyster species in many regions (Huvet et al., 2004;Wang and Guo, 2008a;Camara et al., 2009;Xia et al., 2009), distinguishing oyster species based on shell morphology is unreliable and misleading. Several methods have been established for oyster species determination (Wang and Guo, 2008a;Cordes and Reece, 2009;Xia et al., 2009;Wang et al., 2014); however, simpler methods with improved targeting for identifying specific oyster species in the local region need to be established.
In view of the high diversity and complicated introduction history of oysters in Zhejiang Province, efforts are needed to maintain the characteristics and diversity of local oysters. In this study, we aimed to clarify the dominant oyster species on the Zhejiang coast, determine the genetic diversity of different oyster populations based on molecular markers, confirm the indigenous species, and establish specific methods for identifying local oyster species. These results may lay the foundation for the subsequent protection and utilization of local oyster germplasm resources and the sustainable development of the oyster industry.

Sample Collection
To identify the main oyster species and their distributions, samples from seven sites ( Table 1 and Figure 1) were collected from July 2019 to August 2020 in the three major aquaculture bays on the Zhejiang coast. Samples from the upper region of the intertidal zone (muddy intertidal) and longline-cultured samples from intertidal wild spats were collected from all the sampling sites. For one site in Xiangshan Bay (Chunhu), samples from the subtidal region were collected from fish culture cages and longline-cultured samples (spat from Fujian Province) were also collected (Figure 2). Thirty oysters of different oyster types were randomly sampled from each location for species identification. To measure growth-related traits, one cultured oyster population from each bay was sampled in August 2020. Shell height, length, and width were measured using vernier calipers (0.01 mm) and total weight was measured using an electronic balance (0.01 g).   After dissection, the mantle edge pigmentation was recorded as yellow or black. The adductor muscle was collected and quickly frozen in liquid nitrogen for subsequent DNA extraction and species identification. For genetic diversity analysis of C. angulata, 88 sequences from Xiangshan Bay (abbreviated as XSB hereafter; n = 37) and Zhangzhou (abbreviated as ZZ hereafter; n = 51), and some COI haplotypes of Langqi dao (LQD, n = 39), Fuqing (FQ, n = 29), and Putian (PT, n = 32) from the Fujian coast were obtained from the NCBI nucleotide database (KP216768-KP216806), and the abundances of each haplotype in the populations were extracted from a previous study (Li et al., 2015). For MNR, 61 sequences from XSB (n = 31) and ZZ (n = 30), as well as some sequences (KY217737-KY217796) obtained from Shantou (ST, n = 27) and Keelung (KL, n = 28) of Guangdong and Taiwan Provinces were used. Haplotype distributions in different populations were included (Grade et al., 2016). For C. sikamea, COI sequences from XSB (n = 63), Sanmen Bay (SMB, n = 46), Yueqing Bay (YQB, n = 50), and ZZ (n = 27), and data from Nantong, China (CHH in the original paper, AB736731-AB736845), and Ariake Sea, Japan (ARK in the original paper, AB736621-AB736685)  were selected for comparison. For MNR, only 147 sequences from XSB (n = 42), SMB (n = 46), YQB (n = 34), and ZZ (n = 25) were included; the geographic information of other samples is listed in Table 1 and Figure 1.
Individual C. gigas oysters sampled from oyster farms in Qingdao, C. ariakensis and C. hongkongensis collected in Yangjiang, Guangdong Province, C. sikamea and C. angulata from XSB, SMB, YQB, and ZZ were selected to verify the method used to distinguish C. sikamea from other species.

DNA Extraction and Species Identification
Genomic DNA was isolated from the adductor muscle using a DNA kit (Tiangen, Beijing, China) according to the manufacturer's instructions. The integrity and concentration of DNA were determined using 1.2% agarose gel electrophoresis and spectrophotometry (Allsheng, Hangzhou, China). Oyster species were identified based on the sequences of nuclear internal transcribed spacer 1 (ITS1) and COI markers. The ITS1 fragment was amplified using ITS1-F (5 -GTTTCCGTAGGTGAACCTGC-3 ) and ITS1-R (5 -ACACGAGCCGAGTGATCCAC-3 ) (Wang and Guo, 2008b). The cycling conditions were as follows: initial denaturation at 94 • C for 5 min, followed by 35 cycles of 94 • C for 30 s, 50 • C for 30 s, and 72 • C for 40 s, with a final extension at 72 • C for 10 min. Some COI fragments were amplified with universal primer pairs: LCO 1490 (5 -GGTCAACAAATCATAAAGATATTGG-3 ) and HCO 2198 (5 -TAAACTTCAGGGTGACCAAAAAATCA-3 ) (Folmer et al., 1994). Polymerase chain reactions (PCR) were performed in a final volume of 20 µL containing 7 µL of distilled water and 10 µL of Taq polymerase premix (TaKaRa, Dalian, China). The cycling conditions included initial denaturation at 94 • C for 5 min, followed by 35 cycles of 94 • C for 30 s, 50 • C for 30 s, and 72 • C for 50 s, followed by a final extension at 72 • C for 10 min. The PCR products were sequenced by Hangzhou Tsingke Company (Hangzhou, China). The sequence was analyzed using BLAST in the NCBI server 1 , and oyster species were identified according to the BLAST results of individual COI sequences. Considering that hybrid individuals have difficulty appearing in nature, two individuals were selected for each species and bay to construct a phylogenetic tree using MEGA (Version 6.0) based on the COI and ITS1 sequences of the same individual to eliminate the possible occurrence of hybrid oysters.

Genetic Diversity Analysis of Main Oyster Species
To analyze genetic diversity, COI and MNR were used as the target regions, and the COI sequences were amplified as described in section "DNA Extraction and Species Identification." Primers for MNR-F (5 -TCACAAGTACATTTGTCTTCCA-3 ) and MNR-R (5 -AACGTTGTAAGCGTCATGTAAT-3 ) were used according to (Aranishi and Okimoto, 2005;Grade et al., 2016). The PCR reaction volume and cycling conditions were the same as those described above.
To analyze the genetic diversity parameters, we calculated the number of polymorphic sites, number of haplotypes, haplotype diversity, average number of nucleotide differences, and nucleotide diversity for COI and MNR using DnaSP software (version 5) (Librado and Rozas, 2009).
Population differentiation was evaluated by computing the global fixation index F ST using the Arlequin software (version 3.5) (Excoffier and Heidi, 2010). To explore the proportion of genetic variation among the populations, we conducted an analysis of molecular variance (AMOVA) in Arlequin. The significance of both the F ST and AMOVA tests was assessed using 1000 permutations. A median joining (MJ) network of haplotypes was constructed using PopART (version 1.7 2 ).
To analyze the demographic history of C. angulata and C. sikamea, Tajima's D (Tajima, 1989) and Fu's F (Fu, 1997) were calculated, and a mismatch distribution using DnaSP was assessed to test whether they experienced population expansion.
Methods for Identifying C. sikamea ITS1 sequences from NCBI were analyzed (KX345544-KX345547 and KF370432-KF370439 for C. sikamea and KU726931-KU726939 and AB735529-AB735533 for C. angulata). The HindIII digestion site AGGCCT appeared in the sequence of C. angulata, but not in that of C. sikamea (Supplementary Figure 2). The ITS1 sequence was amplified as described in section "DNA Extraction and Species Identification." The enzyme digestion reaction system contained 3 µL of the PCR product, 1 µL of HindIII (NEB, Ipswich, United Kingdom), 5 µL of the enzyme buffer, and 11 µL of distilled water at a final volume of 20 µL. The reaction was performed in a water bath at 37 • C for 60 min.
To design the MNR-specific primers, the complete mitochondrial genome sequences of C. gigas, C. angulata, C. sikamea, C. ariakensis, C. hongkongensis, C. iredalei, and C. nippona were downloaded and aligned (Supplementary Figure 2). Specific primers were designed for identifying C. sikamea based on the mitochondrial non-coding region (MNR-F :5 -TCACAAGTACATTTGTCTTCCA-3 , MNR-S-R: 5 -AGGCTTTCACTCCACTTACT-3 , MNR-R :5 -AACGTTGTAAGCGTCATGTAAT-3 ). The PCR reaction volume and cycling conditions were the same as those used for MNR amplification, except that adding three kind of primers instead of two. The product was tested using 1.2% agarose gel electrophoresis for both enzyme digestion and the specific PCR.

Species Distribution and Proportion
There were mainly two oyster species by molecule method in the sampled region: C. sikamea and C. angulata ( Table 2). At all seven sites, all oysters sampled in the muddy intertidal area were C. sikamea. For longline-cultured oysters from the intertidal wild spat of XSB, 67-70% of the oysters were C. sikamea, and the remaining were C. angulata. The longline-cultured oysters originating from spats bought from Fujian Province by local farmers were all C. angulata. Among the subtidal samples collected from fish culture cages, the Portuguese and Kumamoto oysters comprised 67% and 33% of the samples, respectively. In SMB, for longline-cultured oysters from intertidal wild spat, 87%-97% of the oysters were C. sikamea, whereas the rest were C. angulata. In YQB, all longline-cultured oysters from intertidal wild spats were C. sikamea. ITS1 sequencing and the phylogenetic tree eliminated the possibility of hybrid individuals in the three bays (Figure 3).

Production Performance of the Two Oyster Species
The production performance of 1-year-old oysters is presented in Table 3. The shell height, length, and width and total weight of Kumamoto oysters in the three bays ranged from 32.1 ± 4.9 to 35.2 ± 3.3 mm, 22.3 ± 3.3 to 25.8 ± 3.1 mm, 15.1 ± 2.0 to 17.5 ± 1.9 mm, and 6.4 ± 2.2 to 7.5 ± 2.5 g, respectively. Moreover, the corresponding data for Portuguese oysters in XSB were 44.2 ± 4.3 mm, 30.9 ± 2.5 mm, 18.8 ± 1.9 mm, and 14.2 ± 3.3 g, respectively. The growth-related characteristics of Portuguese oysters in XSB were significantly larger than those of the Kumamoto oysters in the three bays. The shell height:length ratio ranged from 1.33 ± 0.19 to 1.46 ± 0.28, which was not significant between the Kumamoto and Portuguese oysters, whereas the shell height:width ratio was much lower for the Kumamoto oyster than for the Portuguese oyster in XSB, although this result was not significance (p = 0.06).

Mantle Edge Pigmentation
There were some differences in the mantle edge pigmentation of the Kumamoto and Portuguese oysters. Kumamoto oysters were 77-87% yellow and 13-23% black. In contrast, for Portuguese oysters, the mantle edges were 90% black and 10% yellow ( Table 3).

Genetic Diversity of the Two Species Based on COI/MNR Sequences
Two mitochondrial sequences (COI and MNR) were amplified for genetic diversity analysis, and a COI fragment of 603 bp was obtained from 188 specimens of the Portuguese oyster. In total, 49 segregating sites and 44 haplotypes were identified. Haplotype diversity and nucleotide diversity were 0.836 ± 0.046 and 0.353 ± 0.039% for XSB and 0.869 ± 0.031 and 0.428 ± 0.036% for ZZ populations, respectively. For Kumamoto oysters, 543 bp sequences from 366 individuals were analyzed, and 105 segregating sites and 145 haplotypes were defined. Haplotype diversity ranged from 0.929 ± 0.027 to 0.983 ± 0.017 for XSB, SMB, YQB, and ZZ populations, and the nucleotide diversity ranged from 0.586 ± 0.053% to 0.727 ± 0.063% (Table 4).

Genetic Differentiation and Population Expansion
To conduct the AMOVA among C. angulata mitochondrial COI sequences, two groups were defined, one included the XSB population and the other group included LQD, FQ, PT, and ZZ populations. The results showed that the majority of the molecular variance was within groups (99.74%), with −0.62% variance among groups and 0.89% among populations within groups. The variance components of the three origin types were not statistically significant (p > 0.05) in the AMOVA analysis. Similar results were obtained for Kumamoto oysters when defining the two groups as XSB with the combination of SMB, YQB, NT, and ZZ, as shown in Table 6. Ten haplotypes were shared among the five C. angulata populations based on the COI sequence, and five of the most shared haplotypes (Hap1-2 and 4-6) were found in 71% (134/188) of all individuals. Several major haplotypes located in the center, with a few singletons or low-frequency haplotypes linked by short branches. The MJ haplotype network illustrates evident differences in the population structure between the two species. Furthermore, two haplotypes were shared among all C. sikamea populations, and 25 haplotypes were shared among at least two populations and were present in 6.3% and 53% of all individuals (Supplementary Table 1). The MJ network showed numerous minor haplotypes with no high proportions, and fewer haplotypes were shared (Figure 4).
Population pairwise F ST values were generally low and not significant among the five tested Portuguese oyster populations based on the COI sequence (Supplementary  Table 2).
The neutrality test showed that Fu's F was significantly different from zero (p < 0.05) and Tajima's D was not significant (p > 0.05) for the Portuguese oyster based on the COI sequences; both parameters were observed to be significantly negative for the MNR fragment. For the four populations of Kumamoto oysters, both Tajima's D and Fu's Fs tests showed  significant results for the two sequences (Tables 4, 5). The mismatch distribution was unimodal for C. angulata in XSB and ZZ and for C. sikamea populations in the six locations (Supplementary Figure 1).

Establishment of Species Identification Method
Considering that there are mainly two oyster species on the Zhejiang coast, the ITS1-RFLP method was established. The ITS1 product was digested with the HindIII restriction enzyme to obtain 200 and 310 bp fragments from C. angulata, and the product could not be confused for C. sikamea, leaving a 530 bp fragment, as shown in Figure 5.
To distinguish C. sikamea from additional oyster species, specific primers were designed based on the highly variable MNR. A 660 bp fragment was successfully amplified for C. sikamea, with a larger unspecific bands for C. angulata, C. gigas, C. ariakensis, or C. hongkongensis on the agarose gel ( Figure 5D). For samples from XSB, SMB, YQB, and ZZ, a specific band was observed for C. sikamea, while a larger band was observed for C. angulata (Figure 5C).

Dominant Species and Its Distribution
In this study, oyster species were identified using molecular diagnosis based on the BLAST result of the COI sequence of  individuals, and we found that C. sikamea was widely distributed and cultivated in large numbers on the Zhejiang coast. C. sikamea comprised more than 70% of the oyster samples in 14 batches collected from the three bays in this area. Zhejiang Province has an annual oyster production of approximately 223,000 metric tons (Bureau of Fisheries, 2020); thus, we speculate that at least 50,000 metric tons of this yield is C. sikamea. Previous studies have mentioned its occurrence in this area, but its abundance has not been reported (Wang et al., 2013). It is the dominant species on the Zhejiang coast, whereas in the Ariake Sea, where it was first found, and in the South Korean peninsula, it coexists with the Pacific oyster and was once thought to be endangered at these locations (Hedgecock et al., 1999). The natural distribution of this species has been reported to include South Korea and southern Japan (Ariake Sea and Seto Inland Sea) (Banks et al., 1994;Hamaguchi et al., 2013), southern China (Wang et al., 2013), and Vietnam . The Kumamoto oysters were distributed especially in the warmer estuary and intertidal zones. For instance, in Yueqing, samples were collected in the estuary of the Qingjiang River, and all oysters were C. sikamea. This species shows a habitat preference for estuaries with mesohaline environments (Camara et al., 2009). In addition, its more frequent occurrence in the middle and high regions of the intertidal zone in the estuaries of the Ariake Sea has been reported (Camara et al., 2009). Furthermore, the Kumamoto oyster occupies the intertidal zone, and the Suminoe oyster occupies the subtidal zone (Xu et al., 2009). These reports are consistent with our result that oysters collected in the upper region of the intertidal zone were all C. sikamea, and of those collected in the subtidal zone, 67% were C. angulata. Studies have found that the most common small-cupped oyster from Zhejiang to Hainan in southern China is C. angulata, which is also thought to be the main oyster cultured in Zhejiang and Fujian Provinces (Wang et al., 2010;Li et al., 2015). We found that the Portuguese oysters indeed exist in a narrow region and comprise 30% of the wild spat in XSB, which is a long and narrow bay with an inland distance of 60 km, and the salinity is higher (average salinity around 25) than that of SNB and YQB (You and Jiao, 2011). Portuguese oysters have been reported to be distributed as far as the Okinawa (Ryukyu Archipelago) and Shikoku Islands in southern Japan, although the authors were unsure whether the oysters were indigenous (Sekino and Yamashita, 2012;Sekino et al., 2016). The latitude of Zhejiang Province is similar to that of the above locations, and Portuguese oysters prefer high salinity and high-temperature environments (Shi et al., 2019); therefore, we cannot exclude the natural distribution of Portuguese oysters here.
The Kumamoto oyster is very small compared to the Portuguese oyster, with the former weighing half as much as the latter in this study, and this is the primary reason why local farmers transport oysters from Fujian Province. Kumamoto oysters tend to be deeper cupped than other species (Hedgecock et al., 1999;Camara et al., 2009). Our results seem to support the opinion that the shell height:width ratio for C. sikamea is smaller than that of C. angulata sampled in XSB. The shell shape seems to be heritable, although some studies have found that it has low heritability (Sang et al., 2019) or is affected mainly by the type of sediment (Mizuta and Wikfors, 2018).
For approximately 30 years, local oyster farmers in XSB have continually transported oyster seeds from Fujian Province, which is the core distribution area of the Portuguese oyster (Wang et al., 2010). As identified by us, artificial spawning seeds purchased from Fujian Province were all C. angulata, and they were usually transported in September and successfully spawned in the next summer in the local bay; thus, the origin of the Portuguese oyster is not clear. Moreover, we found that the mantle edge color was light for C. sikamea and dark for C. angulata. Mantle edge pigmentation has moderate heritability (0.215 ± 0.092) among Pacific oysters (Xing et al., 2017) and successive selective breeding for the golden-shell strains of this species has yielded lightermargined oysters (Han and Li, 2019). Local farmers usually distinguish oysters by mantle color, as they believe that local oysters have yellow mantle edges, whereas imported ones have black mantle color. Therefore, it is necessary to perform genetic evaluations to analyze the indigenous oysters of the two species.

Genetic Diversity of the Two Species and Identification of Indigenous Species
For Portuguese oysters, the nucleotide and haplotype diversity based on the COI sequence of oysters collected in XSB was lower than that of the sampled populations from ZZ as well as four sites on the Chinese coast (Hsiao et al., 2016), and four populations on the Fujian coast (Li et al., 2015). The significant results of Fu's Fs test and non-significant Tajima's D test showed that the Portuguese oysters in XSB may have experienced a population expansion event, and the unimodal distribution of nucleotide mismatch analysis also supported this conclusion. It seems that Portuguese oysters on the coast of southern China form a large group, and gene flow exists (Hsiao et al., 2016).
Taking into account that farmers in XSB have introduced large amounts of artificially spawned Portuguese oyster spats from many geographic locations ranging from the south to the north of Fujian Province for many years, the introduced or invasive species tend to show low genetic diversity. Oysters have great fecundity and high polymorphism in DNA sequences (Li et al., 2018); thus, even if the Portuguese oysters spat that were sampled in this study are not indigenous (the offspring of the introduced individuals), the genetic diversity should not be very low. Therefore, it is difficult to determine whether the oysters were the original or offspring of the introduced populations (Grade et al., 2016;Hsiao et al., 2016). Oyster larvae have a long swimming period of up to 3 or 4 weeks, and there are many islands along the coast of Zhejiang; therefore, it is possible that ocean currents drive larvae to settle (Sekino and Yamashita, 2012;Sekino et al., 2016). Molecular markers at the genome-wide level may be needed to further determine their exact source, such as with Pacific oysters in the Northern Hemisphere, where a clear lineage of different populations has been provided (Sutherland et al., 2020) and oysters in Europe (Vendrami et al., 2018).
For Kumamoto oysters, the nucleotide and haplotype diversity was much higher in our samples than in the Japanese and Korean populations, as well as in the Portuguese oysters based on the COI sequence . High diversity was found in three Kumamoto oyster populations in Zhejiang Province when analyzing the population history of this species along the coastal regions of the Northwest Pacific (Hu et al., 2018). Significant Fu's Fs and Tajima's D test results showed that Kumamoto oysters on the Zhejiang coast experienced a population expansion event. There was no significant genetic differentiation among the XSB population and the other four populations along the coast of China (NT, SMB, YQB, and ZZ). However, populations sampled along the Zhejiang coast were quite distinct from the Japanese populations, and the MJ haplotype network indicated that fewer haplotypes were shared among all populations and that private haplotypes accounted for a larger proportion . Kumamoto oyster populations on the Zhejiang coast show high genetic diversity comparable to that of other populations along the coast of China, and large-scale artificial spawning and breeding programs have not been conducted as has been done with the Pacific and Portuguese oysters.
Portuguese oysters are thought to be distributed south of the Yangtze River (Wang et al., 2010) and have not been found to the north , although sporadic reports exist . The diluted water of Changjiang may be a barrier for marine animal distribution to the north (Ni et al., 2017), whereas the distribution of Kumamoto oysters seems to cross this barrier, as indicated by the non-significant F ST , implying its tolerance to low salinity (Camara et al., 2009). Our results indicate that the Kumamoto oysters have high genetic diversity and are likely an indigenous species along the coast of Zhejiang. Portuguese oysters have relatively low genetic diversity, but their sources require further exploration.

Distinguishing Kumamoto Oysters for Better Utilization of Resources
Conducting species identification to eliminate other oyster species before artificial breeding and exploring simple methods of species discrimination are highly important for protecting and utilizing Kumamoto oyster germplasm resources. Although the shell shape of Kumamoto oysters is somewhat different from that of other species, it remains difficult to classify oyster taxonomy according to morphology (Wang et al., 2004). Our results revealed a sympatric distribution between the Kumamoto and Portuguese oysters and between the Hong Kong and Suminoe oysters on the southern coast of China (Wang and Guo, 2008a;Xia et al., 2009).
The results of ITS1-RFLP showed that after digesting the PCR product of the ITS1 sequence with HindIII, the Kumamoto and Portuguese oysters could be clearly distinguished. The sequence alignment results also supported this finding ( Supplementary  Figure 2), as being an effective method for distinguishing between the Kumamoto and Pacific oysters and their artificial hybrids (Xu et al., 2019). This may reflect the subspecies relationship and sequence similarity of the Pacific and Portuguese oysters. The use of PCR-RFLP analysis of several gene fragments for distinguishing seven Crassostrea oysters from the South China Sea showed wider applicability, whereas it was somewhat difficult to utilize two or more restriction enzymes with one to five fragments after digestion (Xia et al., 2009).
The MNR sequence could be used to distinguish the Kumamoto oysters from the five main species of Crassostrea oysters, and the results of four populations (XSB, SMB, YQB, and ZZ) also demonstrated the versatility of this method. Previous results based on the MNR sequence have shown its potential in discriminating different populations of the Pacific and Portuguese oysters (Aranishi and Okimoto, 2005;Grade et al., 2016). It is not surprising that there are great differences between species because of the high polymorphism of the sequences. A multiplex PCR method based on COI was powerful and widely used for differentiating five Crassostrea species (Wang and Guo, 2008a) and we applied this method at the beginning of the species identification process; however, when distinguishing between Kumamoto and Portuguese oysters in XSB, sometimes the latter species exhibited two specific bands, which muddled the determination. Kumamoto oysters were mistaken as Pacific oysters from the Ariake Sea in Japan and introduced to the United States in the 1940s (Hedgecock et al., 1999). At that time, molecular identification methods were not established; therefore, the American population of Kumamoto oysters has been contaminated by Pacific oysters because of purposeful or unintentional hybridization (Hedgecock et al., 1999;Camara et al., 2009). These newly established species identification methods are more practical locally and it has a certain value for the identification of Kumamoto oysters in other regions.
In this study, we found abundant C. sikamea resources distributed along the Zhejiang coast, and these oysters preferred intertidal and estuarial habitats. It is possible that at least 50,000 tons of the annual oyster production is C. sikamea, making this region the world's largest Kumamoto oyster farming area. Our results indicate that Kumamoto oysters have high genetic diversity and should be an indigenous species along the coast of Zhejiang. Portuguese oysters have relatively lower genetic diversity than Kumamoto oysters, but their sources require further exploration. Established species identification methods are of practical significance and are useful for broodstock discrimination. These results have enriched our understanding of the distribution and biological characteristics of Kumamoto oysters, and the genetic diversity assessment is beneficial for better protection and utilization of local oyster germplasm resources and for the development of the oyster industry.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, MW675439-MW675650 for MNR, MW663490-MW663763 for cytochrome c oxidase subunit I (COX1) gene and MW481329-MW481338 for ITS1.

AUTHOR CONTRIBUTIONS
QX and ZL conceived and designed the study. SL and HX collected and sampled the oysters. SL performed the statistical analyses. SL and QX wrote and polished the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was financially supported by the Project "Studies on Technologies for Disease Mitigation in Marine Mollusk Aquaculture" of "3315" Innovative Team of Ningbo City, the fundamental research foundation for Provincial University supported by Zhejiang Wanli University, Earmarked Fund for China Agriculture Research System (CARS-49), and the Science and Technology Planning Project of Ningbo City (2019C10046).