A comprehensive survey of the oyster reef ecosystem along Guangdong coastlines was conducted from July 2022 to February 2023, covering 51 sites ( Supplementary Table 1 ; Figure 1 ). Herein, to encompass the oyster’s extensive natural distribution, we adopted a broad definition of oyster reefs for this study. These reefs constitute biogenic structures formed either by dense oysters’ colonies with substantial vertical reliefs —commonly termed oyster reefs or beds—or secondary structures known as oyster aggregations, where adherent oysters settle upon either hard or soft substrates (Beck et al., 2009; Baggett et al., 2014). Site selection was based on historical distribution records in intertidal areas, including mangroves, mudflats, estuaries, and even aquaculture farms (Cheng et al., 2022). Adjustments to the initially selected sites were occasionally made when necessary to account for changes in the actual oyster reef distribution. Two forms of survey, namely qualitative and quantitative survey were conducted at each site during low tide.

At each site, GPS coordinates were recorded to map the distribution of oyster populations within a defined survey polygon. To comprehensively document diversity of oysters in the region, varying morphologies of oysters based on shell features, coloration and diagnostic pattern were collected. Totally 742 samples by qualitative survey were collected from 51 survey sites. In order to describe the original features of oysters’ specimen, specimens were photographed in situ before being removed from their substrate and individually stored in labeled plastic containers. Ultimately, the flesh was separated by sterilized scalpel and then stored at -20°C for subsequent laboratory analysis. All samples of qualitative survey were employed to molecular identification, and haplotypes separation as well as genetic diversity. Shells of all specimens were deposited in the South China Sea Monitoring Center, State Oceanic Administration.

Each oyster reef site was classified into upper, middle, and lower tidal zones based on natural tidal patterns. In the upper and lower tidal zones, three quadrats of almost 0.1 m × 0.1 m were laid out by quadrats placement methodology at random, adhering to the same guideline of intertidal zone research (Lu et al., 2024), while in the middle tidal zone, nine quadrats were collected in same method. The number of quadrats was based on the ground that nine quadrats informed by species accumulation curve was sufficient to collect all species of each site based on previous research ( Supplementary Figure 1 ). All macrobenthic organisms, including oysters and algae within the quadrats, were removed using a brush and screwdriver and fixed in 95% ethanol. A 15 cm ruler was included in each quadrat and photographed to provide a scale reference, and each photograph was analyzed to calculate exact sample area of quadrats using ZEISS Axio Vision image processing software.

Owing to huge amounts of oysters detected by quadrats layout, tentative morphological identification was performed via the quantitative survey. Crassostrea species were identified based on established shell morphology descriptions (Huber, 2010; Wang et al., 2010). Within the genus Saccostrea, species of S. echinata, S. malabonensis, and S. scyphophilla were distinguishable. For example, S. echinata displays dense hyote (tubular) spines on its right valve and S. malabonensis exhibits a curved shell margin coupled with strong ribs forming a distinctive plicated structure, while S. scyphophilla is characterized by an orderly and jagged shell margin and purple spots on right valves. Despite this, owing to vague morphological distinction of other Saccostrea, we were unable to accurately identify them to species-level. Consequently, to ensure consistency, unresolved Saccostrea species were provisionally classified as S. non-mordax. Finally, an integrated approach was employed, wherein morphological results were guided and further validated by molecular results in the qualitative survey.

Oyster abundances were determined based on the number and weight of individuals, which were then transformed into density (ind/m²) and biomass (g/m²) metrics (Baggett et al., 2014). Biomass were measured by whole specimen with shells in wet weight after removing surface attachment and blotting excess moisture. Spatial distribution maps were generated using ArcGIS 10.7 to visualize these data.

Description of distinct shell characteristics for unresolved oysters included external and internal features. External features examined included the area of left valve attachment, the pattern of radial ribs, grooves, marginal plication, and growth lamellae, alongside the occurrence of growth lines and hyote spines ( Supplementary File Figure 4 ). Internal features considered included the presence or absence of chalky deposits and chomata, the shape of adductor muscle scar and ligament (Littlewood, 1994; Sekino and Yamashita, 2016). The shell height was measured by vernier calipers with 0.1 mm accuracy, and the whole weight was measured by electronic scales with 0.01 g accuracy (Zhang et al., 2017). These characteristics facilitated the preliminary identification of oysters during qualitative surveys and were conducive to the discovery of new taxa within the reefs.

For samples collected by qualitative survey, total genomic DNA was extracted from the adductor muscle using the DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer’s instructions. Mitochondrial DNA fragments of 16S rRNA and COI were amplified using the primer pairs ArL(F) (5’-CGCCTGTTTATCAAAAACAT-3’) and BrH(R) (5’-CCGGTCTGAACTCAGAT CACG-3’) (Palumbi et al., 1996), as well as LCO1490 (5’-GGTCAACAAATCATAA AGATATTGG-3’) and HCO2198 (5’-TAAACT TCAGGGTGACCAAAAAATCA-3’) (Folmer et al., 1994). Polymerase chain reactions were performed in a 25 µl reaction mixture containing 13.38 µl ddH₂O, 2.5 µl Ex Taq Buffer (TaKaRa HS), 1.5 µl of each primer (10 mM), 2 µL template DNA, and 0.12 µl Ex Taq HS, using a BIO-Metro thermal cycler. The PCR conditions for 16S rRNA included an initial denaturation at 94°C for 3 min, followed by 40 cycles of 94°C, 45°C, and 72°C for 40, 40, and 45 s, respectively, with a final extension at 72°C for 7 min. For COI, the PCR protocol included an initial denaturation at 94°C for 3 min, followed by 5 cycles of 94°C, 45°C, and 72°C for 40, 40, and 50 s, respectively. This was followed by 35 cycles of 94°C, 51°C, and 72°C for 40, 40, and 50 s, respectively, with a final extension at 72°C for 7 min. PCR products were confirmed using 1% agarose gel electrophoresis and sequenced by Sangon Biotech Company (Guangzhou, China). The obtained sequences were initially analyzed using BioEdit 7.0. Molecular identification was conducted based on BLAST searches against the NCBI database. More importantly, the identification was further confirmed by phylogenetic tree constructed by correct and non-local reference sequence ( Table 1 ). Haplotypes were identified using DNASP v. 5.10.01 (Librado and Rozas, 2009) and used to calculate haplotype diversity. Each identified haplotype was designated following a standardized format, combining “Hap number” with “ML+serial numbers” for qualitative surveys. The compiled sequences of S. sp.1 were deposited in NCBI under accession numbers PQ452111 and PQ459553–PQ459559 for the COI and 16S rRNA genes, respectively. Sequences of S. sp.2 were assigned under accession numbers PQ455391–PQ455392 and PQ459552 for the COI and 16S rRNA genes, respectively.

DNA sequences from this study, combined with reference sequences downloaded from NCBI, were analyzed using Phylosuite v. 1.2.3 and MEGA v. 7.0.26 (Kumar et al., 2016; Zhang et al., 2020) ( Table 1 ). Several plug-in programs were utilized, including MAFFT for sequence alignment, TrimAI for trimming redundant sequence at both ends, Gblocks for removing divergent and ambiguously aligned blocks, and ModelFinder for model selection (Gerard and Jose, 2007; Capella-Gutierrez et al., 2009; Standley, 2013). Following sequence procession, phylogenetic analysis was conducted using Bayesian inference (BI) and ML methods. A Bayesian tree was constructed with a Markov Chain Monte Carlo simulation run for 2 × 10⁶ generations, while ML analysis was performed using IQ-Tree with bootstrap support based on 5,000 replications (Huelsenbeck, 2012). Sequence divergence among haplotype groups was calculated in MEGA v. 7.0.26 using Kimura’s two-parameter (K2P) model (Kimura, 1980).

Molecular analyses and phylogenetic assessments were conducted on representative individuals from oyster samples of distinct shell features. In total, 706 COI sequences and 693 16S rRNA sequences were generated from 742 individuals ( Table 2 ), yielding 240 COI haplotypes and 67 16S rRNA haplotypes ( Figures 2 , 3 ).

Since the topological structures of the IQ-Tree and BI tree were generally congruent at the upper branches, with only minor difference observed at the lower branches at the interspecies level (e.g., the position of S. kegaki and S. non-mordax H) for both mitochondrial DNA, the BI tree was selected for further analysis. The corresponding IQ-Trees are provided in the Supplementary Figures 2 , 3 . Using Striostrea prismatica as the outgroup for BI and ML trees based on COI sequences, 12 oyster species across three genera were identified. The COI fragment used for phylogenetic analysis was 550 bp in length. The BI tree topology revealed that the representative species belonged to three subfamilies. At the species level, the following oysters formed a robust clade: C. sikamea, C. hongkongensis, C. angulata, C. ariakensis, C. dianbaiensis, C. bilineata, S. malabonensis, S. echinata, S. scyphophilla, S. sp. 1, S. sp. 2, and D. cf. crenulifera ( Figures 4A–L ).

Regarding BI and IQ trees based on 16S rRNA, Neopycnodonte cochlearo and Hyotissa hyotis were selected as outgroup species. In total, 12 oyster species across three genera were identified using the 16S rRNA genetic marker. The 16S rRNA fragment used for phylogenetic analysis was 495 bp in length. S. sp. 1 was closely related to S. non-mordax D, forming a sister group with S. non-mordax E, while S. sp. 2 grouped with S. non-mordax H in the 16S rRNA tree. Consequently, the 16S rRNA tree topology was largely congruent with that of the COI tree, except for variations in the interspecies relationships among S. non-mordax lineages within Saccostreinae.

Among the COI sequences, C. sikamea exhibited 112 haplotypes, accounting for 45.71% of the total haplotypes identified in this study. Saccostrea echinata exhibited 47 haplotypes, contributing 19.18% of the total haplotypes. Similarly, S. scyphophilla, with 23 haplotypes, exhibited the highest nucleotide and haplotype diversities with values of 0.0144 and 0.993, respectively. This pattern was consistent in genetic divergence analyses, with S. scyphophilla showing a greater pairwise distance compared to that of other rock oyster species ( Tables 3 , 4 ).

Regarding the 16S rRNA genetic locus, S. echinata and C. sikamea comprised 13 haplotypes each, representing 20.31% of the total haplotypes. S. echinata exhibited the highest nucleotide and haplotype diversities with values of 0.0028 and 0.710, respectively. A similar pattern was observed in genetic divergence analyses, with S. echinata and S. scyphophilla displaying greater pairwise distances compared to that of other rock oysters ( Tables 1 , 4 ).

In this study, we described the shell morphology of the unidentified species, Saccostrea sp. 1 and sp. 2, to tentatively distinguish them from other oyster species. 9 and 23 individuals for COI sequences and 16S rRNA, respectively, were designated to Saccostrea sp. 1. These specimens were collected from rocky shores in Shantou (G3 and G4), Shanwei (G11, G13 and G15), Shenzhen (G21), Chaozhou (G01) and Xuwen (G51). Hap 41 of S. sp.1 for 16S rRNA was identified in all those regions (except Xuwen), making it the most common haplotype of Saccostrea sp. 1. Conversely, two individuals were identified as S. sp. 2, representing 2 and 1 haplotype for COI and 16S rRNA, respectively, via molecular analysis ( Table 3 ). Specimens of S. sp. 2 were sourced from the rocky shore of Jieyang (G08) and Zhanjiang (G45).

The shell morphology of Saccostrea sp. 1 was characterized by its relatively small size compared to that of other rock oysters, with shell heights typically ranging from 20–40 mm ( Figure 4G ). The shell shape was usually elongated and elliptical, tapering towards the posterior end. The left valve was slightly cupped and primarily attached to the substrate. The left valve featured prominent, well-defined ribs that extended to the dorsal shell margin, forming distinct plications. The internal margin of the left valve possessed clearly visible chomata, and chalky deposits were present on the inner surface. The right valve, which featured brownish patches on its surface, was smaller compared to that of the left valve. The growth lines were prominent and concentric, forming a densely brownish growth lamella. Hyote spines were absent on the surface of the right valve. The adductor muscle scar was distinct and horseshoe-shaped, located in the posterior ventral half of the inner right valve. The hinge line was short and obliquely oriented.

The detailed description of S. sp. 2 was limited owing to the scarcity of specimens. The shell shape of S. sp. 2 appeared flat and triangular ( Figure 4I ). There are no hyote spines observed on the outer surface of the right valve in our specimen. The growth lamellae and squamae were densely clustered along the shell margin, with distinct black deposition present. Radial ribs and the grooves between them were relatively faint but noticeable, especially along the shell margin of the left valve. Internally, marginal chomata were clearly present. The internal black and ochre pigment deposition was conspicuous, appearing as large spots. The adductor muscle scar was clearly visible and relatively large.

In total, 33,549 oyster individuals for oyster morphological identification was sampled using a quantitative method to calculate density, biomass, and species dominance. The average density and biomass of C. sikamea were calculated as 1,415.38 individuals/m² and 4,449.42 g/m², making it the most predominated oyster species along the Guangdong shoreline ( Figures 5 , 6 ). This species was found across all survey sites, regardless of habitat type. Following C. sikamea, S. echinata had a mean density and biomass of 281.43 individuals/m² and 1102.29 g/m², contributing 13.11% and 13.42% to the total species richness and biomass, respectively. In this study, C. hongkongensis was the most common and dominant oyster species in the western district of the Pearl River Estuary, with an average density of 93.48 individuals/m² and a biomass of 501.68 g/m². Saccostrea malabonensis, although less abundant with a density and biomass of 40.45 individuals/m² and 370.26 g/m², respectively, was widely distributed across the study area, occurring at 40 of the 53 sites. All ambiguous Saccostrea identified in the quantitative survey (here referred to as S. non-mordax) were present in measurable quantities in our assay. The total number of species recorded at each site ranged from 1 to 6 species by qualitative survey, while G07, G11, G15, and G38 exhibited the highest richness with 6 oyster species found. At least 7 oyster species were found at G07 by quantitative method, being the most abundant area. Our study indicated G12 was the most abundant site with a density of 6779.83 ind/m², while G08 was found to have the lowest abundance with a density of 626.16 ind/m². G32 was the greatest area of biomass with the value reaching around 16000 g/m². In contrast, the biomass of G46 was the lowest with only about value of 2000 g/m² ( Supplementary Table 1 ).

Currently, approximately 88 extant oyster species (Ostreidae) are recognized worldwide (Bayne et al., 2017; Li et al., 2021; Searles et al., 2021). Over 30 oyster species are found along the coast of China (Hu et al., 2019). In our analysis, population composition revealed at least 12 oyster species from three genera (Crassostrea, Saccostrea, and Dendostrea) and three subfamilies based on the phylogenetic systematics of Salvi, Raith, and Guo (Salvi et al., 2014; Raith et al., 2015; Guo et al., 2018) at the oyster reef of Guangdong. Given ongoing debates about the acceptance of new genera within Crassostreinae (e.g., Talonostrea Li and Qi, 1994, and Magallana Salvi & Mariottini, 2016), we continue to use the genus Crassostrea to represent cupped oysters in our study. Furthermore, because C. iredalei is suspected of having cryptic species and remains controversial (Guo et al., 2018), we considered it C. bilineata rather than C. iredalei. Similarly, due to challenges in phylogenetic analysis, Dendostrea was classified as D. cf. crenulifera in this study, as some D. crenulifera sequences uploaded to NCBI may be inconsistent. Further efforts are needed to clarify the relationship within Dendostrea. Our findings also suggest that S. mordax A & B lineage and S. scyphophilla are conspecific, as indicated by the placement of S. scyphophilla within the S. mordax A & B lineage in the phylogenetic tree, with minimal genetic divergence ( Table 4 ). Therefore, it may be more accurate to refer to the entire S. mordax A& B lineage as S. scyphophilla (Raith et al., 2015; McDougall et al., 2024). S. mordax C lineage we haven’t found in our survey was referred to as S. mordoides (Cui et al., 2021).

Genetic validation of rock oysters and phylogenetic divergence underpins conservation, sustainable utilization, and ecological restoration of these species. Saccostrinae remains insufficiently assessed by DNA analysis due to the subtle and unclear morphological distinction among Saccostrea species. Sekino and Yamashita (2016) identify 13 distinct Saccostrea lineages; however, most monophyletic clades have not been assigned valid species names. Currently, only a few of these lineages have been formally named. For example, S. non-mordax B is presumably recognized as S. echinata, while S. non-mordax F and J are assigned to S. malabonensis and S. spathulata, respectively (Cui et al., 2021; Snow et al., 2023).

S. cuccullata was once considered the sole valid Saccostrea species in the Indo-Pacific, and this superspecies was repeatedly used in earlier records of Guangdong province to represent Saccostrea (Harry, 1985; Xia, 2008). However, the geographical range of S. cuccullata is now recognized as being limited to the Atlantic, from the West Indies to the Arabian Peninsula (Huber, 2010), making this name inappropriate for oysters from the Indo-Pacific Saccostrea lineage. Our findings show that the rock oysters we collected did not cluster with S. cuccullata specimens from Madagascar (Volatiana et al., 2015), with a K2P sequence divergence of 5.6%–19.1% in COI and 2.5%–9.4% in 16S rRNA, suggesting that S. cuccullata does not occur along the Guangdong coastline. Another significant issue is the presence of S. glomerata in the northern hemisphere. Historical data have recorded the distribution of S. glomerata, referred to as “tuanju oysters,” along the Guangdong coastline (Phillips and Yim, 1981; Xu and Zhang, 2008; Wang, 2016; Cui et al., 2018). However, Lam and Morton, based on mitochondrial analysis, suggest that the common flattened oyster found in Hong Kong is not S. glomerata. They propose that native S. glomerata might only inhabit the temperate regions of the southern hemisphere (Lam and Morton, 2006). Our findings indicate that S. glomerata does not occur along the Guangdong coastline, as no taxonomic unit clustered with S. glomerata specimens from Southeastern Australia, with a K2P sequence divergence of 9.33%–26.87% in the 16S rRNA gene (Combosch et al., 2017). Considering the flattened shell form, absence of hyote spines, plication of the right valve margin, along with its relatively large size, and attachment to rocks, the oysters previously identified as S. glomerata along the southern China coastline may actually be S. malabonensis. This aligns with our findings, where S. malabonensis was identified as a common species on the rocky shores of Guangdong. Recently, the S. non-mordax lineage F was reclassified as S. malabonensis, supported by robust phylogeny, which matches our findings (Sekino and Yamashita, 2016; Cui et al., 2021). Therefore, we speculate that the occurrence and fertility of S. malabonensis in southern China may have been previously underestimated.

Our findings suggest that S. echinata may represent S. non-mordax A and B, as its tentative haplotypes are nested within lineages A and B in the 16S rRNA gene tree. S. non-mordax A and B exhibit shell morphology with hyote spines and emerge as sympatric species in our study and related research (Lam and Morton, 2006; Snow et al., 2023), with a limited K2P distance between them. This supports the idea that they may share a common lineage.

In our assay, two cryptic species, Saccostrea sp. 1 and sp. 2 were found in oyster reefs along the Guangdong Province. Molecular data indicated that S. sp. 1 and S. sp. 2 were designated to known S. non-mordax lineage D and H, respectively. The average genetic distance between S. sp. 1 and other Saccostrea species ranged from 15.0%–21.7% for the COI gene and 2.5%–23.7% for the 16S rRNA gene. S. non-mordax lineage D showed the closest genetic relationship with S. sp. 1, followed by lineage E, leading us to classify S. sp 1 as lineage D rather than lineage E. This is not the first record of S. sp. 1 in Guangdong, as previous studies have documented its occurrence in Hong Kong, Hainan, and Taiwan (Lam and Morton, 2006), and its distribution has now expanded further into another district of Guangdong. S. sp. 1 appears to be exclusively distributed in adjoining areas of coastal China. The spatial distribution of Saccostrea oysters suggests that the density and biomass of S. sp. 1 may have been underestimated. However, morphological identification of S. sp. 1 remains challenging. The most distinctive shell features of the species, such as its relatively small size and flat shape tapering towards the posterior with a pointed end, are insufficient to distinguish it from other Saccostrea species, highlighting the challenge of identifying species within this genus.

Saccostrea sp. 2 formed a sister clade with known S. non-mordax H and was closely related to S. kegaki. Regarding K2P distance, S. sp. 2 showed the closest relationship to S. non-mordax H, with a genetic distance of 0% for the 16S rRNA gene and 2.00% for the COI gene. The genetic distance supports the intraspecies relationship. In contrast, S. sp. 2 and S. kegaki exhibited greater genetic divergence, with a 4.3% difference in the 16S rRNA gene and 6.3% in the COI gene, reinforcing the distinction between S. kegaki and lineage H. Considering that lineages A, B, and H both share hytote spines, Lam suggests that these lineages align with S. echinata sensu Torigoe (Torigoe and Inaba, 1981). In our phylogenetic analysis, only the IQ-Tree based on the COI gene supported this hypothesis ( Supplementary Figure 2 ). However, another tree construction analysis indicated that lineage H did not cluster with lineages A or B but instead formed an independent lineage with S. kegaki as the basal clade ( Figures 2 , 3 ).

A putative quantitative survey by morphological analysis were conducted to study oyster composition of Guangdong province. Previous studies of China’s southern coastline (Zhou and Jr., 2003; Lu et al., 2024) and the 2007 national coastal survey (Wu et al., 2025) in Guangdong Province reported substantial populations of C. ariakensis and C. angulata (Wang et al., 2010; Xia et al., 2014). In contrast, our study revealed that C. sikamea and S. echinata were the most predominant and endemic species at coastal of Guangdong Province, while historically dominant species like C. ariakensis and C. angulata were comparatively rare. Similarly, C. sikamea was also identified as the dominant oyster species along the Zhejiang coast to the north of Guangdong (Liu et al., 2021). These discrepancies may result from earlier studies relying on morphological traits or limited spatial sampling. Alternatively, climate change, widespread mariculture, and habitat alteration threatening both wild and farmed oyster populations over the past few decades could explain shifts in oyster species composition in this region. On the other hand, C. sikamea owning to its euryhaline and eurythermal characteristic possessed a specific capacity for adaption to coastline of China. A sufficient genetic variability and strength of gene flow also promote its performance of regional acclimation (Liu et al., 2021). Although historically confined to southern China, C. sikamea has recently expanded beyond its traditional range, now aggressively invading artificial shorelines in northern China (Hu and Dong, 2022). Its remarkable regional adaptation and acclimation may play a key role in its expansion. As such, our study provides an updated baseline for detecting future changes of this species and its distribution along the Guangdong coastline. From this view, the current situation underscores the need for proactive measures to safeguard and preserve oyster diversity. Additionally, C. talonata was found in mangrove of Beibu gulf and coastline in Hainan Island by previous surveys (Li and Qi, 1994; Lu et al., 2024), indicating that it was a common oyster in South China Sea. There are no specimens of C. talonata were detected in our survey though. The result indicated C. talonata was unusual species at oyster reefs in Guangdong province, which was identical to Xia’s result (Xia, 2008). The contrast might be due to its sparse distribution and its preference to subtidal area where it adhered to the surface of mesh cages or shells (Wang, 2016).

Notably, wild-type C. hongkongensis were found between Yangjiang and Zhuhai (G24 - G37) mostly situated at Pearl River Delta in our study. However, this area suffered from habitat losses, tourist pressure, and contamination of terrestrial runoff as a result of rapid urban expansion. In consideration of the restricted area of this species and its essential position for traditional aquaculture in Guangdong Province, more efforts are needed to monitor its population dynamic process and genetic resources. Similarly, there was an ancient and naturally occurring oyster reef first documented at site G47 of Xuwen. However, the low diversity and richness of this site might suggest the poor condition recently. Establishing protected zones of this area for strengthening management may be a viable solution. In addition, our study outlined several critical regions that implied high levels of oyster diversity and abundance. At least 7 oyster species were found at site G07, which might suggest an oyster hotspot of species diversity for this area. Other sites of G38, G12 and G32 were another noteworthy area with high species richness, abundance and biomass. Overall, this study revealed some potential places for site selection of oyster habitat management, protection and restoration efforts. More importantly, it provided insights into where or which oyster species was presumably under population outbreak or population bottleneck.