Identification of SNPs and Candidate Genes Associated with Growth Using GWAS and Transcriptome Analysis in Portuguese Oyster (Magallana angulata)

Abstract

Portuguese oyster (Magallana angulata) is one of the most important shellfish species worldwide. Although significant improvements in growth have been achieved through artificial selection breeding, the genetic basis underlying these traits remains unclear. Thus, this study aimed to (i) estimate variation and heritability for growth-related traits and (ii) identify SNPs and candidate genes associated with growth traits in Portuguese oyster. Five growth-related traits, including shell height (SH), shell length (SL), shell width (SW), whole weight (WW), and soft tissue weight (STW), were measured and analyzed in 114 one-year-old individuals from a cultivated population in Fujian Province, China. Through whole-genome sequencing and genotyping, we obtained 8,183,713 high-quality SNPs. Based on the genomic relationship matrix, heritability for the five traits was estimated, ranging from 0.071 to 0.695. Through genome-wide association analysis (GWAS), a total of nine SNPs were identified as significantly or suggestively associated with one of the growth-related traits, each explaining phenotypic variation ranging from 14.13% to 18.56%. Differentially expressed genes (DEGs) between individuals with extreme phenotypes were identified using comparative transcriptome analysis, ranging from 868 to 2274 for each trait. By combining GWAS and comparative transcriptome analysis, a total of seven candidate genes were identified, with biological functions related to growth inhibition, stress response, cell cycle regulation, and immune defense. The associations between the candidate genes and the growth-related traits were validated by using single-marker association analysis in other populations. Based on SNPs in these candidate genes, 16 haplotypes associated with growth-related traits were obtained. This study contributes to a deeper understanding of the genetic mechanisms of growth traits, and provides a theoretical basis and genetic markers for the breeding of fast-growing strains of the Portuguese oyster.

1. Introduction

Oysters are one of the most important marine bivalves in global aquaculture, cultivated on all continents except Antarctica [1]. Due to the innovations in culture techniques, the global production of oysters has been growing rapidly since 1950. Particularly since 1990, production has increased nearly fivefold due to the advances in larval culture techniques within hatcheries [2]. However, with the rapid development of oyster aquaculture, issues such as slowed growth, decreased resistance, and quality degradation have arisen. Therefore, it is currently crucial to cultivate fast-growing, robust, high-quality oysters for a sustainable aquaculture industry.

In the past decades, researchers have undertaken a series of efforts to improve economic traits such as growth, yield, resistance, and shell morphology [3,4,5,6]. Among them, rapid growth is often considered as the primary goal in genetic improvement programs, since the trait might dramatically influence farmer’s profits [7]. For improving growth traits, selective breeding has been shown to be an effective method. Genetic gains per generation in oysters can reach 7.2% to 25.5% through family selection, and 8.8% to 15.2% through mass selection [8]. However, the genetic mechanisms underlying the growth traits in oysters remain poorly understood.

In recent years, with the accumulation of genomic resources for aquaculture species and the development of whole-genome single-nucleotide polymorphism (SNP) genotyping techniques, genome-wide association analysis (GWAS) has been utilized to explore the genetic architecture and causative loci of quantitative traits including growth traits [9], disease resistance [10], timing of sexual maturity [11], nutritional content [12,13], and tolerance to environmental stresses [14]. At the same time, comparative transcriptome analysis can identify genes and genetic markers associated with economic traits at the gene expression level through differential gene expression [15,16]. Integrating these two methods can enhance the accuracy of candidate gene identification and has been applied in exploring the genetic foundation of growth traits in various aquatic species such as in Salmo salar [17], Larimichthys crocea [18], Haliotis discus hannai [19], and Crassostrea ariakensis [20].

Portuguese oyster (Magallana angulata), also known as Fujian oyster, is a species of the genus Magallana with close genetic similarities to Pacific oyster (M. gigas) [21]. It is an important aquaculture species along the coasts of Asia and southern Europe [22]. Especially in China, Portuguese oyster is the oyster with the largest production [23]. Selective breeding programs to improve growth-related traits in Portuguese oyster have been implemented with notable success. In China, a fast-growth strain of Portuguese oyster named Jinli No. 1 was developed after multiple generations of mass selection, exhibiting a 14% increase in soft tissue weight and a 6% increase in shell width compared to the control strain [24]. In a breeding program in Vietnam, the total weight of Portuguese oyster increased by 17.4% after three generations of selection based on the estimated breeding values (EBVs), with positive correlated changes in soft tissue weight and shell morphology [25]. However, knowledge about the genetic basis of growth traits in Portuguese oyster is still very limited, which has hampered further implementation of molecular and genomic selective breeding. Therefore, in this study, we integrated GWAS and transcriptome analyses to identify SNPs and candidate genes influencing growth-related traits in Portuguese oyster to provide a theoretical foundation and genetic markers for the breeding of fast-growing strains of Portuguese oyster.

2. Materials and Methods

2.1. Sample Collection and Growth Phenotyping

In November 2023, we collected one-year-old experimental samples randomly from a cultivated population of M. angulata in Weitou Bay, Fujian Province, China. The ancestors of this population originated from natural seedling populations in Fujian and Guangdong provinces and have been domesticated for multiple generations along the coast of Fujian Province. The oysters were brought back to the laboratory and then washed to remove the attachments. All individuals were measured for five growth-related traits according to the method described previously [26], including shell height (SH), shell length (SL), shell width (SW), whole weight (WW), and soft tissue weight (STW). The sexes of samples were determined by microscopic examination of gonadal tissue smears. The adductor muscle tissues were quickly frozen in liquid nitrogen for 0.5 h, and then stored at −80 °C for DNA and RNA extraction.

2.2. DNA Extraction, Sequencing, and Variant Calling

The adductor muscle tissues of 114 M. angulata were sent to Annoroad Gene Technology (Beijing, China) Co., Ltd. for total DNA extraction and whole-genome resequencing. Whole-genome resequencing was performed using the Illumina NovaSeq X Plus sequencing platform, with the target sequencing depth of 9×. Preliminary quality control of the raw sequencing data was conducted using Fastp v0.22.0 software [27] to obtain clean reads for subsequent genotyping. Clean reads were aligned to the reference genome of M. angulata (GenBank: GCA_025612915.2) [28] using BWA v0.7.17 software [29]. The resulting SAM files were sorted, indexed, and then converted to BAM files using SAMtools v1.10 software [30]. SNP calling on non-duplicated reads was performed using the HaplotypeCaller module of GATK v4.1.6.0 software and the concordance variants were filtered with the parameter “QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0” [31]. SNPs with minor allele frequency (MAF)  <0.05, missing genotype rate >0.05, and Hardy–Weinberg equilibrium (HWE) p  <  1 × 10³ were filtered out using Plink v1.90 software [32]. Finally, genotype imputation for missing sites of high-quality SNPs was performed using Beagle v5.1 software [33]. The distribution and density of high-quality SNPs across the genome were visualized using the CMplot package [34].

2.3. Population Genetic Structure and Linkage Disequilibrium Analysis

Principal component analysis (PCA) based on high-quality SNP data was conducted using PLINK v1.90 software [32], and population structure of the samples was inferred using ADMIXTURE v1.3.0 software [35]. The range of K values tested was set from 1 to 4, with the optimal K value determined based on the lowest cross-validation error. Using PopLDdecay v3.42 software (https://github.com/hewm2008/PopLDdecay, accessed on 20 November 2024), pairwise linkage disequilibrium (LD) and R-squared correlation coefficients (r²) between alleles at different loci were calculated, and an LD decay curve was plotted.

2.4. Descriptive Statistics and Heritability Estimation of Growth Traits

Descriptive statistics for the phenotypes for SH, SL, SW, WW, and STW were performed using the readxl package in Rstudio v4.5.1, and the Kolmogorov–Smirnov normality test was conducted for each trait. Principal component analysis was conducted based on the phenotypic data of all measured samples using the FactoMineR, factoextra, and psych packages in R studio [36]. The get_eigenvalue function is used to calculate the eigenvalues and contribution rates of each principal component, while the get_pca_var function determines the contribution of each phenotype to the principal components. Pairwise correlations among traits were estimated using the Pearson correlation coefficient [37]. The heritability of each studied trait was estimated using the GREML algorithm in GCTA v1.94.1 software [38], with the SNP-based genetic relationship matrix (GRM) used as the random effect. The formula for calculating narrow-sense heritability is as follows:

$$h_{SNP}^2={\displaystyle \frac{{\sigma }_G^2}{{\sigma }_G^2+{\sigma }_e^2}}$$

where ${\sigma }_G^2$ is additive genetic variance and ${\sigma }_e^2$ is residual variance.

2.5. Genome-Wide Association Analysis (GWAS)

GWAS for the five growth-related traits was performed using the linear mixed model (LMM) in GEMMA v0.98.5 software [39]. Principal components (PCs) and sex were used as fixed effects, while the kinship matrix (K) was used as the random effect to effectively control for false-positive associations. Three models of fixed effects [(1) only sex; (2) sex and 1 PC; and (3) sex and 2 PCs] were evaluated by calculating the genomic inflation factor (λ) using GEMMA v0.98.5 software [40]. The model with the λ value closest to 1 was selected as the optimal model for further analysis [41]. The LMM equation is as follows:

$$\mathrm{Y}=\mathrm{W}\alpha +\mathrm{X}\beta +\mathrm{Z}\mathrm{u}+\epsilon$$

where Y is the vector of phenotypic data for all individuals; W is the matrix of covariates (fixed effects); $\alpha$ is the vector of coefficients for the fixed effects; X is the SNP genotype matrix; $\beta$ is the vector of marker effect coefficients; Z is the random effect matrix derived from the kinship; u is the vector of random effects; and $\epsilon$ is the error vector [39]. To mitigate the overly conservative nature of Bonferroni correction, we set 1/N and 0.05/N as the suggestive and significant thresholds for the genome-wide association study, where N represents the total number of SNP markers used in the association analysis [18].

The suggestive and significant SNPs were annotated using SnpEff v5.2 [42]. Manhattan plots and quantile–quantile (QQ) plots were visualized using the Cmplot package in R (v 4.5.1). Finally, genes containing suggestive or significant SNPs were functionally annotated based on protein sequences or domains using the NCBI NR database (https://www.ncbi.nlm.nih.gov/, accessed on 20 November 2024) and the InterPro protein domain database (https://www.ebi.ac.uk/interpro/search/sequence/, accessed on 20 November 2024). The formula for calculating the phenotypic variance explained (PVE) by SNPs is as follows [43,44]:

$${\displaystyle \frac{2{\widehat{\beta }}^{2}MAF(1-MAF)}{2{\widehat{\beta }}^{2}MAF\left(1-MAF\right)+{\left(\mathrm{s}\mathrm{e}\left(\widehat{\beta }\right)\right)}^{2}2NMAF(1-MAF)}}$$

where $\widehat{\beta }$ is the effect value of the SNP in the GWAS; MAF is the minor allele frequency of the SNP; $\mathrm{s}\mathrm{e}\left(\widehat{\beta }\right)$ is the standard error of the SNP’s effect value; and N is the number of samples analyzed in the GWAS.

2.6. RNA Extraction and Transcriptome Analysis

According to the phenotypic data of the collected samples, the top 10 and bottom 10 for each trait were selected, whose adductor muscles were taken and sent to Beijing Annoroad Company for RNA extraction, library construction, and sequencing. RNA-seq was performed on the Illumina NovaSeq X Plus sequencing platform with the target sequencing depth of 9×. Preliminary quality control of the raw sequencing data was performed using Fastp v0.22.0 software [27]. Clean reads were aligned to the reference genome (GenBank: GCA_025612915.2) [28] using HISAT2 v2.1.0 software [45]. The resulting SAM files were sorted, indexed, and then converted to BAM files using SAMtools v1.10 software [30]. Read count matrices were calculated using the run-featurecounts.R script, and gene expression was quantified using Transcripts Per Million (TPM) values.

Differentially expressed genes (DEGs) between the maximum and minimum phenotypes were identified using the DESeq2 package in R (v1.40.2) [46]. Significant DEGs were defined as genes with |log2 fold change| > 1 and p-value < 0.05. Subsequently, significant DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using the ClusterProfiler package (v4.8.2) [47,48,49].

2.7. Candidate Gene Identification

To further identify key candidate genes associated with growth traits, we integrated GWAS and transcriptome analysis to screen for candidate genes affecting growth traits. Genes within a 100 kb window upstream and downstream of suggestive or significant SNPs were compared with DEGs of extreme phenotypes for growth traits. Overlapping genes were identified as candidate genes. Finally, candidate genes were functionally annotated based on protein sequences or domains using the NCBI NR database (https://www.ncbi.nlm.nih.gov/, accessed on 20 November 2024) and the InterPro protein domain database (https://www.ebi.ac.uk/interpro/search/sequence/, accessed on 20 November 2024).

2.8. Variation of Candidate Genes

SNPs in each key candidate gene were extracted using bcftools v1.9 software [30]. A NJ (Neighbor-Joining) phylogenetic tree was constructed using VCF2Dis v1.47 software (https://github.com/hewm2008/VCF2Dis/, accessed on 20 November 2024) and visualized using the iTOL online platform (https://itol.embl.de/, accessed on 20 November 2024). Groups were formed based on the phylogenetic tree, and the significance of differences in growth-related traits between groups was tested using the T-test. Genotype heatmaps of each gene region were constructed using RectChr v1.37 software (https://github.com/hewm2008/RectChr/, accessed on 20 November 2024) to further analyze the association between genotype and phenotype.

Single-marker association analysis for all SNPs in the candidate gene regions was performed using GCTA v1.94.1 software [38]. For SNPs significantly associated with phenotypes (p < 0.05) in the single-marker association analysis, LD blocks were identified using Haploview v4.2 [50]. Finally, haplotype association analysis of the LD blocks was performed using SHEsisPlus (http://shesisplus.bio-x.cn/, accessed on 20 November 2024).

2.9. Validation of Candidate Genes

We validated the effects of candidate genes on growth traits by using an independent cultured population of Portuguese oyster, which was collected in Weitou Bay in November 2022. We collected the growth phenotypic data and SNP genotype data using the previously described methods. We performed the single-marker association analysis for all SNPs within the candidate gene regions using the GCTA v1.94.1 [38] software to assess the association between genotypes of the candidate genes and growth traits in the validation populations.

3. Result

3.1. Descriptive Statistics of Phenotypes

Table S1 summarizes the phenotypic data of growth traits for a total of 114 samples from a cultivated population of Portuguese oyster, including shell height (SH), shell length (SL), shell width (SW), whole weight (WW), and soft tissue weight (STW), with mean values (Mean ± SD) of 77.20 ± 11.27 mm, 43.03 ± 4.17 mm, 31.82 ± 5.89 mm, 52.96 ± 13.87 g, and 9.55 ± 2.94 g, respectively (

Table 1. Descriptive statistics and heritability estimates for each trait.

	SH	SL	SW	WW	STW
Mean	77.20	43.03	31.82	52.96	9.55
Median	77.42	43.20	31.19	52.29	9.35
Min	56.63	27.91	20.31	28.30	3.76
Max	117.25	52.48	50.64	95.34	20.24
SD	11.27	4.17	5.89	13.87	2.94
CV	14.60%	9.69%	18.50%	26.19%	30.77%
Heritability (h2)	0.690	0.118	0.072	0.695	0.071

, Figure 1A). Kolmogorov–Smirnov normality analysis of the phenotypic data for these growth traits showed that all traits followed a normal distribution (p > 0.05, Figure S1). The coefficients of variation for these five traits were 14.60%, 9.69%, 18.50%, 26.19%, and 30.77%, respectively. Additionally, we observed a significant sexual dimorphism in soft tissue weight, with males showing better growth performance than females (an average increase of 25.13%). In contrast, there were no significant differences between sexes in the other growth traits (Figure 2). Principal component analysis (PCA) based on the phenotypes of these five growth traits showed that the cumulative contribution rate of the first three principal components was 88.1%, with the first principal component contributing the most to the total phenotypic variance at 51.1% (Figure S2). Phenotypic correlation analysis showed a positive correlation between all traits except for SH and SW, with the highest correlation observed between WW and STW (Figure 1B). Subsequently, the heritability was estimated with the SNP-based genetic relationship matrix (GRM) used as the random effect for each trait and sex as covariate for STW, which was sexually dimorphic. Among the traits, WW and SH showed high heritability (h² > 0.6), while SL, SW, and STW had low heritability (h² < 0.2) (Table 1).

3.2. Resequencing, Genotyping, Population Structure, and Linkage Disequilibrium Analysis

Whole-genome resequencing was performed on 114 samples, yielding approximately 677.38 GB of clean bases after quality control, with an average sequencing depth of 9.36 ± 0.8× and a Q30 bases rate of 95.25%. After stringent quality control, a total of 8,183,713 high-quality SNPs were identified for subsequent analysis. The SNP density plot shows that the distribution of SNPs across the chromosomes is uniform (Figure S3A). The average number of SNP loci per chromosome is 722,260, with a density exceeding 13 SNP loci/Kb.

Principal component analysis (PCA) shows that none of the samples exhibit significant population stratification within the group (Figure S3B). A similar pattern was observed in the population genetic structure analysis. All samples were classified into the same category, according to the minimum coefficient of variation (CV) error value (K = 1) (Figure S3D,E). When conducting linkage disequilibrium (LD) analysis on the experimental population, it was found that the average r² value rapidly decayed to half of its maximum value (~0.22) within 42 bp, indicating a very fast decay rate (Figure S3C).

3.3. Genome-Wide Association Analysis (GWAS)

To select the most suitable model for genome-wide association analysis, three fixed-factor linear mixed models (LMMs) were firstly evaluated for each growth-related trait. According to the genomic inflation factor (λ), the LMM that included sex, the first two principal components (2PC), and the kinship matrix (K) was selected for SH and SW; the LMM that included sex, the first principal component (1PC), and K was selected for SL and WW; and the LMM that included sex and K was selected for STW (Figure S4).

Setting the suggestive and significant thresholds for GWAS at 1.23 × 10⁷ (1/N) and 6.16 × 10⁹ (0.05/N), respectively, we identified a total of nine SNPs suggestively or significantly associated with growth traits (Lead SNP), including five for SH, three for SL, and one for STW. The phenotypic variance explained (PVE) of these SNPs ranged from 14.13% to 18.56% (Figure 3,

Table 2. List of significant and suggestive loci associated with the growth-related traits in Portuguese oyster.

Traits	Chr	Position	p-Value	PVE (%)	Location	Gene
SH	1	33,653,680	7.56 × 108	14.27	intergenic	-
	4	11,938,717	3.09 × 108	15.1	intergenic	-
	7	52,563,485	5.05 × 108	14.64	intronic	gapr1
	10	13,473,040	7.78 × 108	14.24	intronic	ppa1
	10	13,563,505	3.90 × 108	14.88	intronic	nwd1
SL	4	11,516,450	8.56 × 108	14.13	exon *	mab21l
	4	11,516,476	7.20 × 1010 *	18.56	exon	mab21l
	5	26,098,444	4.15 × 108	14.81	intronic	dspp
STW	9	12,305,944	1.77 × 108	15.58	intronic	lac25b

In the p-value column, * indicates that the value exceeds the significance threshold, while its absence indicates that it surpasses the suggestive threshold. In the Location column, * indicates that the SNP is a missense variant.

). Of these nine SNPs, five are located in the intron regions of the Golgi-associated plant pathogenesis-related protein 1 (gapr1), pyrophosphatase (ppa1), nacht domain and wd repeat-containing protein 1 (nwd1), dentin sialophosphoprotein (dspp), and laccase25 (lac25b) genes. Notably, two of them are located in the exon region of the mab21l/cyclic gmp-amp synthase-like receptor (mab21l), where Chr4-11516450 is a missense mutation, changing from C to A, resulting in the corresponding amino acid changing from threonine (Thr) to lysine (Lys); the remaining two SNPs are located in intergenic regions (Table 2).

3.4. Transcriptome Analysis Between Extreme Phenotypes of Growth Traits

We collected muscle tissue from individuals with extreme phenotypes of SH, SL, SW, WW, and STW for transcriptome analysis. Table S2 summarizes the phenotypic data of these individuals with extreme growth traits. The average SH, SL, SW, WW, and STW for individuals with the largest phenotypes were 99.71 ± 8.49 mm, 49.54 ± 1.28 mm, 41.43 ± 3.35 mm, 79.58 ± 7.54 g, and 14.56 ± 2.50 g, respectively, while those with the smallest phenotypes were 65.00 ± 2.90 mm, 38.24 ± 1.42 mm, 23.72 ± 2.20 mm, 42.37 ± 3.13 g, and 7.07 ± 0.80 g, respectively. In the SH trait group, we identified 1324 differentially expressed genes (DEGs), of which 986 were upregulated and 338 were downregulated. In the SL trait group, 868 DEGs were identified, including 411 that were upregulated and 457 that were downregulated. In the SW trait group, 1081 DEGs were identified, including 799 that were upregulated and 282 that were downregulated. In the WW trait group, 2274 DEGs were identified, including 1899 that were upregulated and 375 that were downregulated. In the STW trait group, 1970 DEGs were identified, with 1508 being upregulated and 462 being downregulated (Figure 4, Tables S3–S7). The DEGs between extreme phenotypes among the five traits include both trait-specific and overlapping genes (Figure 5; Tables S8 and S9). For downregulated genes, there were 277, 359, 222, 241, and 330 genes identified specifically in SH, SL, SW, WW, and STW traits, respectively, and 1 gene (LOC128185975, unannotated) shared with all five traits. Meanwhile, for upregulated genes, there were 499, 270, 518, 626, and 338 identified specifically in SH, SL, SW, WW, and STW traits, respectively, and 5 genes (potassium channel subfamily K member 2, protein rtoA, sulfotransferase 1C4, CUB and EGF-like domain-containing protein 1, and LOC128177207) shared with all five traits.

GO enrichment of DEGs showed that 224, 757, and 686 GO terms were significantly enriched for SH, WW, and STW, respectively, while no GO term was significantly enriched for SL and SW (Figure S5, Tables S10–S12). By comparing the GO enriched terms among SH, WW, and STW traits, 153 overlapping enrichment terms were revealed, involving functions such as ion transmembrane transport, signal transduction, perception and behavior, hormone regulation, neural signal transduction, reproductive regulation, and so on (Table S17). KEGG enrichment of DEGs for each trait showed that 41, 2, 32, and 49 pathways were significantly enriched for SH, SW, WW, and STW, respectively, while no term was significantly enriched for SL (Figure S5, Tables S13–S16). By comparing the enriched KEGG pathways among SH, WW, and STW traits, 21 overlapping enrichment pathways were revealed, involving functions such as signal transduction, hormone regulation, and neural signaling (Table S18).

3.5. Candidate Gene Identification

Trait-associated mutations may directly alter gene function by changing protein structure, thereby affecting the phenotype of individuals. Therefore, we identified mab21l as a candidate gene, for it contains a missense mutation in exons significantly associated with the SL trait (Table 2 and

Table 3. List of key candidate genes affecting growth traits.

Traits	Lead SNP	Gene ID	Gene Symbol	Description
SH	10-13473040	LOC128167546	sstr2	somatostatin receptor type 2-like
	10-13473040	LOC128165816	v-SNARE	v-SNARE-like, coiled-coil homology domain
SL	4-11516450	LOC128181972	mab21l	MAB21L/Cyclic GMP-AMP synthase-like receptor
STW	9-12305944	LOC128162734	lac25a	laccase-25
	9-12305944	LOC128162736	crfr2	corticotropin-releasing factor receptor 2-like
	9-12305944	LOC128163477	trim36	E3 ubiquitin-protein ligase TRIM36-like
	9-12305944	LOC128163983	hgsnat	Heparan-alpha-glucosaminide N-acetyltransferase-like

).

Trait-associated mutations may also affect gene function by altering gene expression levels, thereby influencing individual phenotypes. Therefore, we identified the intersection of genes that were located within 100 kb upstream and downstream of lead SNPs and DEGs between extreme phenotype groups for each trait as candidate genes. Within 100 kb upstream and downstream of lead SNPs, 48, 15, and 16 genes were identified for SH, SL, and STW traits, respectively (Table S19). The intersection of these genes and the DEGs between the extreme phenotype groups resulted in two (sstr2 and v-SNARE), zero, and four (lac25a, crfr2, trim36, and hgsnat) genes (Table 3, Figure S6A–C).

3.6. Variation in Candidate Genes in the Population

We further examined the genetic structure of candidate genes and their potential impact on growth traits. For sstr2, the phylogenetic tree constructed based on SNPs within the gene revealed three main branches within the population (Figure 6A), with distinct genotypic characteristics for each branch (Figure 6B). Individuals in branches G1 and G3 showed significant differences in SH (Figure 6C and Figure S7), suggesting a potential impact of sstr2 variation on SH traits. For crfr2, the phylogenetic tree constructed based on SNPs within the gene showed four main branches within the population (Figure 6D). The genotypic differences among these branches were visualized using a genotype heatmap, clearly displaying the distinct branches (Figure 6E). There were significant differences in the mean STW values among the branches (Figure 6F). Besides STW, the G2 branch of crfr2 also exhibited relatively higher performance in SH, SL, and WW (Figure S8), suggesting pleiotropy in the regulation of growth traits by crfr2. The remaining candidate genes also showed significant differences among genetic branches in at least one growth-related trait (Figures S9–S13).

Finally, we identified trait-associated SNPs (p < 0.05) within the region of candidate genes through single-marker analysis and used them to construct haplotypes (

Table 4. Haplotype analysis of candidate genes.

Trait	Gene	Hap ID	Haplotype	Count	Beta	SE	p
SH	sstr2	Hap.1	CT	172	−4.15	1.58	9.00 × 103
			TT	11	4.57	3.56	2.02 × 101
			TG	42	4.26	1.94	3.00 × 102
SL	mab21l	Hap.2	GCAT	206	4.47	0.90	2.47 × 106
			TAGG	17	−5.51	0.97	1.07 × 107
		Hap.3	AGCAG	209	4.61	0.89	1.01 × 106
			TAAGA	14	−5.80	1.06	2.88 × 107
STW	lac25a	Hap.4	TGG	205	−1.86	0.57	1.00 × 103
			CAT	20	2.14	0.62	7.75 × 104
		Hap.5	AT	180	−1.98	0.46	3.11 × 105
			GT	7	0.06	1.15	9.59 × 101
			GC	39	2.41	0.49	3.40 × 106
	crfr2	Hap.6	GAC	180	−1.95	0.41	5.02 × 106
			GAT	26	1.39	0.52	8.00 × 103
			ACT	21	1.84	0.62	3.00 × 103
		Hap.7	GTA	200	−1.72	0.54	2.00 × 103
			ACG	26	1.83	0.60	3.00 × 103
		Hap.8	TA	208	−2.07	0.70	3.00 × 103
			AC	20	2.07	0.70	3.00 × 103
		Hap.9	GC	145	−1.57	0.36	3.29 × 105
			GT	41	−0.22	0.49	6.48 × 101
			AT	41	2.34	0.41	8.52 × 108
		Hap.10	CA	206	−1.80	0.58	2.00 × 103
			AG	22	1.80	0.58	2.00 × 103
		Hap.11	GTT	182	−1.48	0.49	3.00 × 103
			TGG	46	1.48	0.49	3.00 × 103
	trim36	Hap.12	CT	207	−2.08	0.69	2.00 × 103
			TA	21	2.08	0.69	2.00 × 103
		Hap.13	AA	206	−1.60	0.65	1.40 × 102
			CT	22	1.60	0.65	1.40 × 102
	hgsnat	Hap.14	CT	210	−2.32	0.64	4.23 × 104
			TC	18	2.32	0.64	4.23 × 104
		Hap.15	CG	205	−1.90	0.57	1.00 × 103
			TT	23	1.90	0.57	1.00 × 103
		Hap.16	TA	192	−1.11	0.52	3.50 × 102
			GG	34	1.07	0.55	5.50 × 102

Beta represents the strength and direction of the association between each haplotype and the phenotype.

). For SH, one haplotype was constructed in sstr2, while no haplotypes were formed by trait-associated SNPs in v-SNARE. For SL, two haplotypes were constructed in mab21l. For STW, 14 haplotypes were constructed across four genes (Figure S14). Through haplotype association analysis, we identified one haplotype significantly associated with SH, two haplotypes significantly associated with SL, and thirteen haplotypes significantly associated with STW (Table 4, Figure S15).

3.7. Validation of Candidate Genes

The results of the single-marker association analysis in the validation populations identified a certain number of SNPs (p < 0.05) associated with growth traits in all seven candidate gene regions (

Table 5. The most significantly associated SNPs with growth traits in the candidate gene regions, based on the validation population.

Traits	Gene	Chr	Pos	REF	ALT	p
SH	sstr2	10	13,403,077	C	T	7.97 × 103
	v-SNARE	10	13,378,409	A	G	1.61 × 102
SL	mab21l	4	11,516,890	T	A	1.64 × 102
STW	lac25a	9	12,223,193	C	T	7.98 × 103
	crfr2	9	12,197,864	A	T	5.15 × 103
	trim36	9	12,328,176	A	T	1.75 × 103
	hgsnat	9	12,363,518	A	T	1.66 × 102

and Table S20), confirming the association of the candidate genes and the growth traits.

4. Discussion

Growth is one of the most interesting economic traits in aquaculture. To support the implementation of genetic improvement programs in oyster, genetic analysis of growth-related traits has been conducted, but mainly in Pacific oyster [26,51,52,53]. Portuguese oyster is also a very important oyster species, but the reports on the genetic analysis of its growth-related traits are still limited [25,54]. Therefore, in this study, we investigate the phenotypic variation and genetic basis of five growth-related traits in Portuguese oyster by combining phenotypic, genomic, and transcriptomic data.

Based on phenotypic data, we revealed high variation (CV values ranging from 9.69% to 30.77%) in the growth-related traits of the investigated population. Especially for STW, the phenotype CV value is as high as 30.77%. Heritability estimation based on SNP-based relationship matrix revealed that WW and SH were high-heritability traits, and SL, SW, and STW were low-heritability traits. High phenotypic variation and significant heritability of the growth trait in Portuguese oyster have also been reported previously. For example, Wu et al. (2019) reported that the CV values ranged from 16.46% to 61.07%, while the heritability ranged from 0.116 to 0.595 for growth-related traits in 342-day-old Portuguese oyster [54]. Vu et al. (2020) reported that CV values ranged from 19.10% to 50.60%, while the heritability ranged from 0.10 to 0.24 in 270-day-old Portuguese oyster [25]. Similar results have been reported in other oysters, such as in Pacific oyster and European oysters (Ostrea edulis) [55,56]. The high phenotypic variation and significant heritability of growth-related traits in oysters indicate a strong potential for improvement through selective breeding.

In addition, we revealed significant sexual dimorphism in the STW trait of Portuguese oyster at one year of age, with the males being approximately 25.13% higher than females. In oyster, only one case of sexually dimorphic growth has been reported. Baghurst et al. (2002) observed that among two-year-old Pacific oysters, females grew faster than males, with an advantage of up to 18% in shell length and up to 15% in total weight [57]. The two studies differ in species, geographic distribution (southern vs. northern hemisphere), and season, making direct comparison challenging; however, both suggest sexual dimorphism in oyster growth. The pioneering findings suggest that more detailed investigations on sexual dimorphism in oyster are warranted to assess the potential value of implementing sex-controlled breeding in Portuguese oyster.

In recent years, GWAS has been widely used for genetic analysis of economically important traits including growth in aquatic animals, such as in Salmo salar [58], Coilia nasus [59], Litopeneaus vannamei [60], Magallana gigas [26], and Patinopecten yessoensis [61]. In the present study, we conducted GWAS on the growth traits of Portuguese oyster for the first time, and a total of nine SNPs associated with growth traits were identified, with PVE values ranging from 14.13% to 18.56%. In studies of growth-related traits in aquatic animals, significant association loci often show low significance and a low PVE. For instance, Zhu et al. (2023) identified 3 significant and 104 suggestively associated loci, with the PVE ranging from 1.56% to 3.68% for body length and weight by GWAS in Golden pompano (Trachinotus ovatus) [62]. Yang et al. (2020) identified 5 significant and 18 suggestive QTL in five growth-related traits in Brown-marbled grouper (Epinephelus fuscoguttatus), with a PVE ranging from 10.24% to 15.46% [63]. Liu et al. (2022) identified 11 significant associated loci in 10 growth-related traits in Pacific oyster (Magallana gigas), with a PVE ranging from 8.09% to 10.81% [64]. The low significance and PVE of the associated loci for the growth-related traits indicates that the traits may be affected by multiple minor-effect loci through a complex genetic network [59].

Among the nine associated SNPs identified in this study, only two were located in exon, and the remaining seven were located in intronic or intergenic regions, suggesting that regulatory elements may play an important role in the phenotypic variation of growth-related traits [65]. Therefore, we identified candidate genes associated with growth-related traits by integrating the results of GWAS and comparative transcriptome analysis between groups with extreme phenotypic values, resulting in six candidate genes including sstr2, v-SNARE, lac25a, crfr2, trim36, and hgsnat. In addition, mab21l was also identified as a candidate gene due to two associated SNPs being identified in its exon, including a missense mutation (Chr4-11516450), even though it was not a DEG between groups with extreme phenotypes. Among the candidate genes, sstr2 encodes the somatostatin receptor (SSTR), which has been confirmed as the primary inhibitor of growth hormone (GH) release in many vertebrates [66,67,68]. However, the mechanisms by which somatostatin (SS) regulates growth hormone (GH) secretion differ among species. For instance, in pigs and baboons, low concentrations of SS can significantly stimulate the release of GH [69,70,71]. In turbot, following exposure to 20 μg/L of tralopyril, a significant increase in total body length was observed in the experimental group, with growth hormone levels 30.11% higher than the control group, along with elevated levels of sstr2 expression [72]. In this study, the expression of sstr2 was significantly upregulated in the maximum phenotype groups of SH, WW, and STW traits, suggesting that sstr2 may positively affect the growth of Portuguese oyster to some extent.

Trim36 is a member of the TRIM ring protein family. It is a microtubule-associated E3 ubiquitin ligase that not only affects cell cycle progression, but is also essential for early embryonic development and plays a key role in regulating cell growth rate and stability. For example, in human ESCC cells, overexpression of trim36 inhibited cell proliferation, while silencing of trim36 led to the opposite effect. In Xenopus laevis, trim36 plays an essential role in the formation of the dorsal–ventral axis [73,74]. In Cynoglossus semilaevis, trim36 is primarily expressed in the gonads and is thought to play a role in sex determination and differentiation [75]. In this study, the expression of trim36 was significantly upregulated in the group with a maximum phenotype for STW traits, suggesting that the gene may also have a growth-promoting effect in Portuguese oyster.

Corticotropin-releasing factor (CRF) is a key regulator of adaptive responses to internal and external stresses, which can activate the hypothalamic–pituitary–interrenal (HPI) axis and subsequent stress response by interacting with CRF receptors (CRFR1 and CRFR2) [76,77]. Under prolonged stress and high cortisol levels, fish redirect their energy allocation from supporting growth to prioritizing survival and restoring normal physiological states. This shift reduces the energy available for growth, thereby slowing the growth rate [78]. In addition, it is known that the CRF system plays a key role in regulating feeding behavior in vertebrates [79,80]. For example, in Dabry’s sturgeon, both short-term and long-term fasting induced a decrease in crfr2 mRNA expression, leading to a 14% reduction in body weight; however, upon resumption of feeding, crfr2 mRNA expression significantly increased, accompanied by an 8% increase in body weight [81]. Similar phenomena have been observed in studies involving mice and frogs [82,83]. In this study, the genetic structure of crfr2 significantly affected the phenotypic values of STW, SH, SL, and WW traits. Moreover, the expression of the crfr2 gene was significantly upregulated in the maximum phenotype group for WW and STW traits. These results indicates that crfr2 is a gene substantially affecting growth-related traits in Portuguese oyster, possibly by adjusting energy distribution or by regulating food intake.

The other four candidate genes, lac25a, hgsnat, v-SNARE, and mab21l, have been reported to be involved in immune defense in shellfish [84,85,86]. For instance, in Patinopecten yessoensis, the expression level of laccase25 was significantly upregulated after infection with Polydora, impacting the immune response of the organism [87]. In Mytilus unguiculatus, Hgsnat was identified as a hub of upregulated genes after being exposed to the toxin-producing algae Alexandrium catenella [86]. In Eriocheir sinensis, SNARE was found to resist Ediocheuris infection by enhancing host–cell phagocytosis [85]. In Litopenaeus vannamei, Mab21cp participated in the antiviral response by regulating the STING pathway [84]. The association of growth-related traits with immune-related genes suggests a potential genetic correlation between growth traits and disease resistance traits in Portuguese oyster. In fact, genetic correlations between growth and disease resistance traits have been observed in fish, both negative [88,89,90] and positive [91,92]. Vehviläinen et al. (2012) proposed that negative genetic correlations may result from genetic compensation between traits, as they compete for limited resources within the organism; meanwhile, positive genetic correlations may arise from general vigor, where individuals that grow faster also perform better in fitness traits [91]. Understanding whether there is a positive or negative correlation between growth traits and fitness traits in Portuguese oyster is crucial for their application in selective breeding programs.

It is noteworthy that there was an overlap among different growth-related traits in terms of DEGs between extreme phenotypic groups, as well as in the GO terms and KEGG pathways that they enriched. The overlapping pathways are involved in cation transmembrane transport, G protein-coupled signaling, hormone secretion and action, neuropeptide signaling, and more (Tables S17 and S18), which may be due to the pleiotropic effects of genes and pathway functions. For example, G protein-coupled receptors (GPCRs) transmit numerous extracellular signals by coupling with G proteins and inhibitory proteins to trigger intracellular signaling, participating in various physiological processes including perception, behavior and emotion regulation, immune system activity regulation, autonomic nervous system transmission, and homeostasis regulation [93,94]. This overlap in gene expression networks may explain correlations between traits, indicating that selection for one trait can impact others [95]. This pleiotropic genotypic differentiation provides breeders with an opportunity to improve multiple traits simultaneously.

5. Conclusions

There were significant additive genetic variations in the SH, SL, SW, WW, and STW traits of the investigated population of Portuguese oyster, with heritability ranging from 0.071 to 0.695. By using GWAS and transcriptome analysis, a total of 9 SNPs, 7 candidate genes, and 16 haplotypes associated with growth-related traits were identified. These candidate genes are involved in growth regulation, cell cycle regulation, stress and feeding regulation, and immune defense. In addition, the overlaps of gene regulatory networks between different grow-related traits were initially revealed. These findings will help deepen the understanding of the genetic mechanisms underlying growth traits and provide a theoretical basis and genetic markers to improve growth-related traits in Portuguese oyster.