Correction: Identification of sex-linked markers and genes in Portuguese oyster (Magallana angulata)

Xiaoyan Zhou¹Yue Ning²Jingyi Xie¹Yi Han¹Chaonan Tang¹Caiyuan Su¹Qijuan Wan¹Qisheng Wu²Xiang Guo²Jianfei Qi²YIZHOU KE¹Hui Ge^2*Mingyi Cai^1*

• ¹Jimei University, Xiamen, China
• ²Fisheries Research Institute of Fujian, Xiamen, China

1IntroductionSex is a critical trait in aquatic animals, as sexual dimorphism significantly influences growth, size, shape, color, behavior, and ornamentation (Martínez et al., 2014; Horne et al., 2020; Li et al., 2022). These sex-specific differences have driven global efforts to develop sex-controlled breeding programs aimed at generating populations with a high proportion of either males or females, depending on production goals (Budd et al., 2015; Li et al., 2022). Achieving such goals necessitates a thorough understanding of the genetic architecture underlying sex determination, alongside the identification of reliable sex-specific or sex-linked genetic markers. Recent advances in molecular genetics and genome sequencing technologies have enabled the discovery of numerous sex-specific and sex-linked DNA markers in over 100 aquatic species, including fish, amphibians, crustaceans, echinoderms, and mollusks (Li et al., 2022). These markers have been successfully applied in molecular-assisted sex-controlled breeding programs for species such as Nile tilapia (Oreochromis niloticus), rainbow trout (Oncorhynchus mykiss), and yellow catfish (Pelteobagrus fulvidraco) (Sun et al., 2014; Guiguen et al., 2018; Mei and Gui, 2018). Mollusks,the second-largest animal phylum after arthropods, exhibit extraordinary diverse and plastic sex determination mechanisms (Carlton, 1999; Gofas and Zenetos, 2003; Collin, 2013; Breton et al., 2018). Understanding these systems is critical not only for the development of sex-controlled breeding strategies in aquaculture but also for elucidating the broader evolutionary trajectories of sex chromosomes in invertebrates. Early studies on molluscan sex determination utilized karyotyping (Vitturi and Catalano, 1988; Rolan-Alvarez et al., 1996; Vitturi et al., 1998; Baršienė et al., 2000), sex ratio analysis in family lines, triploidy and gynogenesis (Allen et al., 1986; Guo and Allen, 1994) to infer underlying sex determination mechanisms and sex chromosomal configurations. Recent advances in genomic technology have provided deeper insights into molluscan sex determination. For example, studies in abalones have identified sex-linked markers across three species, revealing interspecies turnover of sex chromosome (Luo et al., 2021; Weng et al., 2022; Tang et al., 2025). Research on scallops has localized sex determination regions associated with transposons, indicating chromosomal transitions (Han et al., 2022). Collectively, these findings highlight the evolutionary diversity and genomic complexity of sex determination systems in mollusks, and emphasize the importance of integrative genomics for unraveling conserved and lineage-specific regulatory mechanisms.Oysters, the world’s largest marine aquaculture group, are distributed across all continents except Antarctica (Forrest et al., 2009). They are predominantly dioecious, capable of bidirectional sex reversal, and occasionally exhibit hermaphroditism (Park et al., 2012). Since the 1970s, researchers have studied the genetic regulation of oyster sex determination. In Crassostrea virginica and M. gigas, paternal effects on sex ratios in half-sibling families were observed, leading to the proposal of a “three-locus” and “single-locus” models of sex determination (Haley, 1977; Haley, 1979; Guo et al., 1998; Hedrick and Hedgecock, 2010). More recently, microsatellite and GBS sequencing have identified major sex-linked QTLs in M. gigas, highlighting the role of genetic factors in sex determination (Guo et al., 2012; Han et al., 2021). The Portuguese oyster (M. angulata), a warm-water subspecies of the Pacific oyster, is distributed along the southeastern coasts of China and Vietnam (Wang et al., 2010). It was inadvertently introduced to France in the late 19th century and later to other European countries for commercial aquaculture (Grade et al., 2016; Hsiao et al., 2016). Portuguese oysters exhibit sexual dimorphism in growth and nutritional traits, making them ideal candidates for sex-controlled breeding (Zeng et al., 2013; Qin et al., 2023; Xie et al., 2024). Furthermore, extreme sex ratios are often observed in some families, such as those highly biased toward males or females, suggesting that it is feasible to use marker-assisted selection to develop high-androgenic or high-femaleogenic lines (Vu et al., 2020). However, the genetic basis of sex determination in Portuguese oyster remains unexplored. In this study, we applied genome-wide association studies (GWAS), comparative transcriptomic analysis, and genetic structure analysis to identify sex-associated loci and genes in a cultured population of Portuguese oysters. The results provide a foundation for elucidating the genetic mechanisms underlying sex determination and advancing sex-controlled breeding strategies for this species.2Article typesOriginal Research3Materials and MethodsSample Collection, and Whole-Genome Resequencing Samples for GWAS analysis and transcriptomics were collected from a culture strain in Weitou Bay, Fujian Province, China, in November 2023. A total of 112 one-year-old oysters (57 females and 55 males), were included in the study. After collection, the samples were brought back to the laboratory, where the sex was identified and recorded through microscopic examination of gonadal tissue smears. Muscle and gonadal tissues were then sampled, rapidly frozen in liquid nitrogen, and stored at -80°C for subsequent DNA and RNA extraction. The adductor muscle samples were sent to Annoroad Gene Technology Co., Ltd. (Beijing, China) for total DNA extraction and whole-genome resequencing at an average sequencing depth of 9×. Sequencing data were processed with Fastp v0.22.0 (Chen et al., 2018) to remove low-quality reads, resulting in clean reads. The clean reads were aligned to the Portuguese oyster genome (GenBank: GCA_025612915.2) using BWA v0.7.17 software (Li and Durbin, 2009). SNPs were identified with GATK v4.1.6.0 (Yang et al., 2011) and filtered using Plink v1.90 based on standard thresholds (MAF ≥ 0.05, missing genotype rate ≤ 5%, and Hardy-Weinberg equilibrium P-value ≥ 1e-3). Missing genotypes were imputed with Beagle v5.1(Browning and Browning, 2016), and high-quality SNPs were used for further analysis.3.1Identification and Validation of Sex-Associated SNPs and genes GWAS was performed with Plink v1.90 (Purcell et al., 2007) using a logistic regression model. A Bonferroni correction was applied to set a significance threshold (0.05/total SNPs). SNPs surpassing this threshold were considered significantly associated with sex, and results were visualized as Manhattan plots using the Cmplot package (v4.5.1) in R. To validate the GWAS results, sex-associated SNPs were used for principal component analysis (PCA) and phylogenetic tree analysis. PCA was conducted with Plink v1.90 (Purcell et al., 2007), and a maximum likelihood (ML) tree in NEWICK format was generated using IQ-TREE v2.2.0 (Minh et al., 2020). The tree was visualized on the iTOL platform (https://itol.embl.de/). SNPs that grouped individuals by sex in both PCA and phylogenetic analyses were deemed reliable.Sex-associated SNPs were functionally annotated using the M. angulata genome (GenBank: GCA_025612915.2) and its annotation file (GCF_025612915.1-RS_2022_12) with the snpEff.jar package. Additionally, the 30 kb upstream and downstream regions of sex-associated SNPs were annotated to identify sex-associated genes.3.2Genome-Wide Analysis of Intersex Genetic Differentiation and Sex-Biased SNP Distribution The FST value, a measure of genetic differentiation, was calculated for 20 kb genomic windows with a step size of 10 kb using VCFtools v0.1.16 (Danecek et al., 2011). Windows with FST values in the top 1% were identified as outliers, indicating potential involvement in sex differentiation. The genome-wide distribution of FST values was visualized using the Cmplot package in R (v4.5.1). Sex-biased SNPs were identified based on genotype frequency differences between sexes. Specifically, SNPs where the reference allele homozygous genotype frequency (0|0) exceeded 85% in females and the alternate allele heterozygous or homozygous genotype frequency (0|1 or 1|1) exceeded 85% in males were classified as sex-biased (Luo et al., 2021). The genome-wide distribution of these sex-biased SNPs was visualized in 20 kb windows using the Cmplot package in R.3.3RNA Extraction and Transcriptome AnalysisGonadal tissues from four female and four male samples were randomly selected for RNA extraction, library construction, and sequencing, conducted by Annoroad Gene Technology Co., Ltd. Transcriptome sequencing was performed on the Illumina NovaSeq X Plus platform at a depth of 9×. After quality control with Fastp v0.22.0 (Chen et al., 2018) , the clean reads were aligned to the M. angulata genome (GenBank: GCA_025612915.2) using HISAT2 v2.1.0 (Pertea et al., 2016). Alignment files were processed for sorting, indexing and format conversion using SAMtools v1.10 (Li et al., 2009). Gene expression levels were quantified as TPM (transcripts per million) using the run-featurecounts.R script (Liao et al., 2014) (https://github.com/mskcc/RNAseqDB/blob/master/run-featurecounts.R). Differentially expressed genes (DEGs) between sexes were identified with DESeq2 v1.40.2 (Love et al., 2014), with the thresholds |log2 FC| ≥ 1 and P < 0.05. Identified DEGs were further subjected to GO and KEGG enrichment analysis using the ClusterProfiler package (v4.8.2) (Yu et al., 2012).3.4Identification and Variation Analysis of Sex-Determining Candidate GenesAn integrative approach was applied to identify genes involved in sex determination by combining genetic variation, gene expression data, and clustering analysis. Sex-associated genes identified from GWAS were compared with significant DEGs from transcriptome analysis, and overlapping genes were visualized using the VennDiagram package (Chen and Boutros, 2011). The expression levels of candidate genes, measured as TPM, were displayed as heatmaps generated with the pheatmap package (Zhao et al., 2021). For further analysis the genetic divergence pattern of the candidate genes, SNPs within the candidate gene regions were examined extracted and subjected to principal component analysis (PCA) using Plink v1.90. Genetic distance matrices were generated with VCF2Dis v1.47, and neighbor-joining (NJ) trees were constructed and visualized on the iTOL platform (https://itol.embl.de/). Genes that exhibited distinct sex-associated clustering in both PCA and NJ trees were identified as key sex-determining candidate genes.For the representative key gene, SNP variations were further analyzed. The upstream and downstream SNP distributions of these key genes were visualized using ggplot v2.2.1, while local GWAS, heterozygosity, and nucleotide diversity analyses were conducted using Plink v1.90 and VCFtools v0.1.16. Genotype patterns were presented as heatmaps generated with RectChr v1.37 (https://github.com/hewm2008/RectChr/), providing a detailed visualization of genetic variation around the candidate genes.3.5Construction and Validation of Sex-Associated Haplotypes Single-marker association analysis was performed for all SNPs within the critical candidate gene region using GCTA v1.94.1 software (Yang et al., 2011), identifying SNPs significantly associated with sex (P < 0.05). Linkage disequilibrium (LD) blocks were then detected using Haploview v4.2 (Barrett, 2009), and haplotype associations with sex were analyzed using SHEsisPlus (http://shesisplus.bio-x.cn/).To validate the reliability of sex-associated haplotypes, an independent validation population consisting of 50 one-year-old angulata individuals (25 males and 25 females) were collected from Weitou Bay in November 2022. Whole-genome resequencing was conducted, followed by SNP calling and haplotype association analysis using the same protocols as those applied to the discovery population. This approach ensured methodological consistency across datasets and allowed for a robust assessment of the reproducibility of sex-linked genetic signals.4Results4.1Whole-Genome Sequencing and Variant CallingA total of 112 individuals (57 females and 55 males) were subjected to Illumina paired-end sequencing, generating 748.67 GB of raw data. After quality control, 713.71 GB of clean bases were retained (NCBI PRJNA1185961, Supplementary Table S1), with an average sequencing depth of 9×. From the clean data, 8,183,713 high-quality SNPs were identified through SNP calling and filtering, which were used for subsequent analyses.4.2Identification of Sex-Associated VariantsA GWAS analysis using a logistic regression model identified 1,613 SNPs significantly associated with phenotypic sex (P < 6.31 × 10⁻⁹). Among these, 1,599 SNPs (99.13%) were concentrated within a 2.52-Mb region on Chr.9 (7.99–10.51 Mb), six SNPs (0.37%) were located on Chr.6, and the remaining eight SNPs were distributed across Chr.2, Chr.5, Chr.7, Chr.8, and Chr.10 (Figure 1A, Supplementary Table S2).Principal component analysis (PCA) based on the identifiedsex-associated SNPs revealed that PC1 accounted for 77.56% of the variance, effectively separating the samples into two distinct groups: a predominantly male group (92.59% male; 50 males and 4 females) and a predominantly female group (91.38% female; 53 females and 5 males) (Figure 1B). A similar pattern was observed when PCA was conducted using only the sex-associated SNPs on Chr.9, where PC1 explained 76.41% of the variance, segregating the samples into two groups with male proportions of 92.45% and 10.17%, respectively (Supplementary Figure S1B).A maximum likelihood tree constructed using genome-wide sex-associated SNPs grouped the 112 individuals into two major branches (Figure 1C). One branch consisted predominantly of a female monophyletic group (designated FI), containing 42 females and 3 males (93.33% female). The other branch comprised a male-dominant monophyletic group (designated MI), including 49 males and 2 females (96.08% male), along with a polyphyletic group enriched for females (designated FII), which contained 13 females and 3 males (81.25% female). Subbranching within FI and FII indicated high genetic diversity at sex-associated loci. Additionally, a maximum likelihood tree based on sex-associated SNPs from Chr.9 revealed that females clustered near the base of the tree, while males were positioned in more derived branches.To further investigate sex-associated genomic regions, we employed two complementary approaches: (1) identification of highly differentiated genomic windows between sexes, defined as those with intersex FST values in the top 1%; and (2) analysis of distribution of sex-biased SNPs, defined as loci where the homozygous reference genotype (0|0) frequency exceeded 85% in females, while the heterozygous or homozygous alternative genotype (0|1 or 1|1) frequency exceeded 85% in males. Using the FST approach, 564 highly differentiated windows were identified, with 369 windows (65.4%) located on Chr. 9. The remaining windows were distributed across the other nine chromosomes, ranging from 2 to 52 windows per chromosome (Figure 2A, Supplementary Table S3). For sex-biased SNPs, 82 loci were identified (Figure 2B, Supplementary Table S4), 80 of which (97.56%) were concentrated in the 8.19–11.24 Mb interval on Chr. 9. This interval overlapped with the distribution of sex-associated SNPs identified through GWAS (Figure 2B, Figure 3E-F). Notably, no loci meeting the opposite criteria (0|0 in males and 0|1 or 1|1 in females) were detected. Moreover, a local zoomed-in view of the sex-linked region on Chr. 9 revealed a plateau-like distribution pattern of both sex-associated and highly differentiated SNPs (Figure 3A–D) , suggesting the existence of a relatively broad and stable sex-linked genomic region.4.3Annotation of Sex-associated VariantsSNP annotation revealed that 1,484 of 1,613 sex-associated SNPs (92.00%) were located within genes, while 129 SNPs (8.00%) were in intergenic regions. Among the intragenic SNPs, 678 (45.69%) were in introns, 452 (30.46%) in upstream regions, 177 (11.93%) in exons, 80 (5.39%) in downstream regions, and 97 (6.54%) in untranslated regions (UTRs) (Figure 4). Gene annotation within 30 kb upstream and downstream of these SNPs identified 57 genes, including 54 on Chr.9 and 3 on Chr.5 (Supplementary Table S5). Of these, 20 genes contained missense mutations (Supplementary Table S6). 4.4Comparative Transcriptomic Analysis of Male and Female GonadsTranscriptomic analysis of gonadal tissues (testis and ovary) from four male and four female individuals generated a total of 51.64 Gb of raw sequencing data. After quality filtering, 50.82 Gb of clean data were retained, with Q30 scores exceeding 94.53% for all libraries, indicating high sequencing quality (Supplementary Table S7). Pearson correlation analysis revealed higher similarity within the same sex compared to between sexes, underscoring significant differences in gene expression patterns between testis and ovary (Supplementary Figure S2).Differential expression analysis identified 8,635 differentially expressed genes (DEGs) between testis and ovary (|log2 FC| ≥ 1, P < 0.05), with 3,792 genes upregulated in the ovary and 4,843 genes upregulated in the testis (Figure 5). GO enrichment analysis showed that ovary-upregulated genes were primarily associated with chromosomal separation, chromatin binding, nuclear division, DNA replication, repair, and conformational changes (Supplementary Table S8, Supplementary Figure S3A). In contrast, testis-upregulated genes were enriched in processes related to G protein-coupled receptor signaling, transmembrane signaling, and cytoplasmic composition (Supplementary Table S9, Supplementary Figure S3B). KEGG pathway enrichment analysis revealed distinct functional roles for DEGs in two sexes. Ovary-upregulated genes were enriched in pathways related to DNA replication, repair, transcription, and the cell cycle (Supplementary Table S10, Supplementary Figure S3C). Testis-upregulated genes were associated with transport processes, transcription factors, ribosome biogenesis, and G protein-coupled receptor signaling pathways (Supplementary Table S11, Supplementary Figure S4D). Additionally, key sex differentiation pathways, such as Wnt/β-Catenin, MAPK, and cAMP/PKA signaling, were also found in KEGG enriched pathways (Supplementary Table S10-S11), highlighting the potential roles of DEGs in regulating sex determination and differentiation. 4.5Identification of Candidate GenesAmong the 57 sex-associated genes identified through GWAS analysis and SNP annotation, only 12 were differentially expressed between male and female gonads, including 10 male-upregulated genes and 2 female-upregulated genes (Figure 6, Supplementary Table S12). All 12 genes were located within the sexually-differentiated region on Chr. 9. Notably, ABCF2 and SPATA48 were consistently and highly expressed in the ovary, while SCP, 5-HTR, ADPK1-a, ADPK1-b, and CCKAR were stably and highly expressed in the testis (Figure 6B).To assess whether the 12 identified genes exhibit sex-associated genetic divergence pattern, NJ trees and PCA analyses were independently constructed based on SNPs within each gene. The results revealed that SNPs from 5-HTR, ADPKD1-a, SCP, and SLC28A3 demonstrated strong ability to distinguish individuals by sex in both NJ trees and PCA analyses, whereas SNPs from the other genes did not show similar discriminatory power (Figure 7, Supplementary Figures S4-S5). These findings suggest that these four genes may represent the earliest instances of genetic differentiation between sexes and may play a pivotal role in the sex determination process of Portuguese oysters. Consequently, 5-HTR, ADPKD1-a, SCP, and SLC28A3 are proposed as key candidate genes for future studies on the mechanisms of sex determination in this species (Table 1).4.6Haplotype Construction and Association Analysis of the 5-HTR GeneAmong the four key candidate genes, 5-HTR exhibited the strongest association with sex and contained a missense mutation, making it a representative gene for detailed SNP variation analysis (Table 1). This gene is located on Chr.9 between 8.37 and 8.43 Mb, comprising five exons and four introns (Figure 8A). Within this gene, 1,552 SNPs were identified, of which 239 (15.40%) were significantly associated with sex (Figure 8B). Local GWAS analysis revealed that sex-associated SNPs were distributed throughout the gene and its flanking genomic regions, with only 150 SNPs located in exon. Notably, a nonsynonymous mutation in exon 5 (Chr 9: 8,369,495 A>G) resulted in an amino acid substitution from isoleucine (Ile) to threonine (Thr) (Supplementary Table S6).Interestingly, males exhibited higher genetic diversity (π) and observed heterozygosity (Ho) within the 5-HTR gene, particularly in the second intron region, suggesting the potential existence of a primitive male heterozygous pattern (Figures 8C, E). Genotype heatmaps further demonstrated that males had significantly more heterozygous locithan females across the four candidate genes (Figure 8E, Supplementary Figure S6). To validate the reliability of the resequencing data, Sanger sequencing of five female samples targeting the 5-HTR gene at Chr. 9: 8,369,495 A>G was performed, confirming the accuracy of the genotyping results (Supplementary Figure S7 and S8).Single-marker analysis identified 26 SNPs within 5-HTR that were significantly associated with sex. These SNPs formed six haplotypes, each containing sex-biased alleles, with association significance ranging from 9.41 × 10⁻³² to 6.92 × 10⁻¹³ (Table 2,Supplementary Figure S9 and S10). Validation in an independent population confirmed that all six haplotypes were significantly associated with sex. However, differences in significance levels and allele frequencies between populations suggest that the genetic basis of sex determination exhibits both consistency and variability across populations (Figure 9).5DiscussionOysters are dioecious, with bidirectional sex reversal potential across breeding seasons. In 1977, Haley identified paternal effect in half-sib families of C. virginica and proposed a “three-locus model”, in which oyster sex is determined by the additive effects of three loci (Haley, 1977; Haley, 1979). Subsequently, Guo et al. also observed paternal effect on sex ratios in half-sib families of M. gigas and suggested a “single-locus model”, positing that sex is controlled by one pair of alleles where genotype MF denotes males and FF may denote either males or females (Guo et al., 1998). More recently, two independent studies identified a major QTL linked to sex determination in M. gigas, using microsatellite markers and genotyping-by-sequencing (GBS), respectively (Guo et al., 2012; Han et al., 2021). Consistent with these findings, our study identified a major sex-associated genomic region on Chr.9 in M. angulata, supporting the hypothesis that sex determination in oysters involves a key genetic locus, particularly in certain families or populations.Our study, along with the findings of Han et al., identified a major sex-determining locus on Chr.9 in Portuguese oysters and Pacific oysters, respectively (Han et al., 2021). However, the previously identified sex-linked genes PSMB3 and GPER1 were not located on Chr.9in Portuguese oysters, suggesting that the primary sex-linked loci on Chr. 9 may differ between these two species. This observation aligns with the well-documented variability in sex-determination systems among closely related species and even among populations within the same species, as observed in plants (Liston et al., 2014; VanBuren et al., 2015; Massonnet et al., 2020), bony fish (Lee et al., 2003; Franchini et al., 2018) , amphibians (Jeffries et al., 2018), and shellfish (Han et al., 2022; Tang et al., 2025). Additionally, reduced significance of the 5-HTR gene in our validation population compared to the GWAS population suggests population-specific differences in genetic factors influencing sex determination. Such differences align with the frequently observed variation in sex ratios among families that were established in breeding programs (Guo et al., 1998; Vu et al., 2020),and highlight the complex and dynamic nature of sex determination mechanisms in oysters. Further studies warranted to comprehensively explorethese inter- and intra- species variations more comprehensively. Chromosomal rearrangements are known to play a critical role in the origin and evolution of sex chromosomes. Examples include fragment duplication on LG9 in Oryzias latipes (Kondo et al., 2006), transpositions in salmonids (Bachtrog et al., 2008), and Y chromosome inversions in Drosophila miranda that restricted recombination with the X chromosome and accelerated differentiation (Bachtrog et al., 2008). In our study, the sex-linked region on Chr.9 exhibited distinct boundaries, elevated FST values, and a high density of sex-dimorphic loci over a 2.5 Mb interval, suggesting that one or more chromosomal rearrangements may have contributed to the formation of this sex-linked genomic region. We also observed multiple low-coverage gaps within this region, which is potentially due to the limitations of Illumina sequencing in high GC content, centromeric, or heterochromatic regions (Henikoff et al., 2001; Treangen and Salzberg, 2012; Ross et al., 2013). These findings highlight the need for further investigation into the structural genomic mechanisms underlying sex determination in oysters.Sex-determining genes typically exhibit two critical characteristics: (1) the presence of sex-associated genetic variations, such as SNPs exhibiting significant differences between sexes; and (2) differential expression in male and female gonads, reflecting their functional role in sex differentiation. These dual features have been validated in well-known sex-determining genes, such as dmy in medaka (Matsuda et al., 2002), hsd17b1 in Japanese sea bream (Koyama et al., 2019), and fshr in flathead grey mullet (Ferraresso et al., 2021). By integrating GWAS and transcriptomic data, we identified 12 candidate genes that exhibit both sex-associated SNPs and differential expression. In addition, sex-determining genes are often located in the earliest evolutionary strata of sex chromosomes, representing the initial and strongest differentiation between sexes (Ferguson-Smith, 2006). Consequently, we further refined our candidate gene analysis by examining the association of their genetic structure with sex, narrowing down the 12 candidates to 4 key genes: 5-HTR, ADPKD1-a, SCP, and SLC28A3. Among the four candidate genes identified, 5-HTR, ADPKD1-a, and SCP-a appear to participate in a regulatory framework involving neuroendocrine signaling and calcium dynamics. 5-HTR encodes a serotonin receptor, which regulates intracellular Ca²⁺ levels, activates the cAMP/PKA signaling cascade, and influences downstream transcriptional programs (Nichols and Nichols, 2008). In invertebrates such as pearl oysters and giant clams, serotonin (5-HT) signaling promotes oocyte maturation and gamete release(Wang and He, 2014; Yang et al., 2024), while in mosquitoes, it regulates ovarian development (Ling and Raikhel, 2018). In vertebrates, particularly fish, 5-HT has also been shown to modulate GnRH neuron activity and gonadotropin secretion, indicating a conserved role in the neuroendocrine regulation of reproduction (Prasad et al., 2015). ADPKD1-a belongs to the polycystin family and encodes a component of mechanosensitive calcium channels (e.g., Polycystin-1/2) that mediate Ca²⁺ influx in response to hydrodynamic or osmotic stimuli, thereby influencing cell fate decisions and gene expression (Wang et al., 2019; Lemos & Ehrlich, 2018). SCP-a, which encodes a sarcoplasmic calcium-binding protein, contributes to intracellular calcium buffering and signal precision, maintaining the spatial-temporal integrity of Ca²⁺ signaling (Asunción-Alvarez et al., 2024). Together, these genes may form a calcium-centered regulatory module that transduces environmental cues into developmental responses. Supporting this, studies in alligators and fish have shown that calcium signals activated by environmental changes regulate sex-determining genes such as Dmrt1 and Cyp19a1a via MAPK, CREB, and redox-sensitive pathways (Castelli et al., 2020; Heck & Woodward, 2025). In parallel, SLC28A3-like, a sodium-dependent nucleoside transporter, may link reproductive metabolism to sex-specific cell proliferation, particularly under high biosynthetic demand during gonadal development (Gray et al., 2004). Collectively, these findings suggest that 5-HTR, ADPKD1-a, and SCP-a may constitute a neuroendocrine–calcium signaling module that integrates environmental and genetic signals to orchestrate sex determination, gonadal development, and sexual phenotype maintenance in M. angulata.In the present study, most of the identified sex-associated SNPs (92.0%) were located in non-coding regions, primarily in introns and flanking regulatory areas, suggesting a regulatory role in gene expression rather than direct protein modification. Non-coding SNPs may regulate gene expression through mechanisms such as cis-regulatory elements, splicing regulation, and epigenetic modifications (Pagani et al., 2000; Pastinen and Hudson, 2004; McVicker et al., 2013). In addition, higher heterozygosity observed in males within the 5-HTR region supports a male heterogametic-like segregation pattern, consistent with evolutionary trade-offs in sexual function optimization (Lesaffre et al., 2024).Studies have demonstrated that oyster sex determination is influenced not only by genetic factors but also by environmental conditions. For instance, high temperatures have been shown to accelerate gonadal development in M. gigas, resulting in increased expression of male-biased genes and decreased expression of female-biased genes (Santerre et al., 2013). Similarly, starvation stress delays gonadal development in M. gigas, skews the sex ratio toward males, and alters both the methylation levels and expression of key sex-related genes (Sun et al., 2023). In this study, neither genome-wide sex-associated SNPs nor SNPs within key candidate genes were able to fully differentiate the sexes of the studied individuals. This finding suggests that environmental factors, in addition to geneticfactors, may also contribute to sex determination in the studied population. These results highlight the complexity of oyster sex determination, which likely involves polygenic regulation and intricate gene-environment interactions. 6Conclusion and PerspectivesThis study provides the first population-level genomic insight into sex determination in the Portuguese oyster (Magallana angulata). Through GWAS and transcriptomic integration, we identified a major sex-linked region on Chr.9 and four candidate genes (5-HTR, ADPKD1-a, SCP, and SLC28A3) with sex-associated SNPs and differential expression, suggesting their involvement in sex regulation. Among them, 5-HTR showed strong male-specific heterozygosity and haplotype signals. However, a key limitation is that all transcriptomic samples were derived from mature gonads, which may have missed early-acting sex-determining genes. We also observed reduced validation significance in an independent population, indicating possible population-specific variation. Future studies should focus on early gonadal stages and functional validation using gene editing and expression profiling. Together, our findings lay the groundwork for understanding sex determination in oysters and support the development of sex-controlled breeding strategies in aquaculture.