L. maculatus is a non-endangered, non-protected species in China. All animal procedures were approved by the Animal Ethics Committee of the Chinese Academy of Fishery Sciences (Approval No. 2011AA1004020012).

In August 2024, 150 healthy spotted sea bass at the same growth stage were randomly selected from the same aquaculture ponds of Yueshun Aquaculture Co., Ltd. in Zhuhai city, Guangdong Province, China, where water quality conditions were uniform and individual feed intake could not be quantified. Their morphological characteristics are detailed in Supplementary Table S1. The fish were anesthetized with 200 mg/L MS-222, followed by body weight measurement and phenotypic trait recording in accordance with GB/T 18654.3-2008 (Figure 1), with appropriate adjustments made based on differences in body shape. Caudal fin tissues were excised using sterile scissors, immediately placed in 1.5 mL EP tubes containing anhydrous ethanol, and stored at -20 °C. Dorsal muscle tissues and intestinal contents were separately collected into 2 mL cryovials, flash-frozen in liquid nitrogen, and subsequently transferred to -80 °C freezers for preservation. After sample collection, caudal fin tissues were submitted to Shijiazhuang Boruidi Biotechnology Co., Ltd. for DNA extraction and genotyping. Muscle tissues were sent to Guangzhou Yixi Testing Co., Ltd. for fatty acids analysis using gas chromatography. Lipids were extracted and methylated using the acid hydrolysis–extraction method in accordance with GB 5009.168-2016, followed by fatty acid analysis using gas chromatography–mass spectrometry (GC–MS; Agilent 7890–5975). Chromatographic separation was performed on an Agilent J&W DB-FastFAME capillary column (30 m × 0.25 mm, 0.25 μm). The injector temperature was set at 250  °C, with an injection volume of 2.0 μL and a split ratio of 50:1. Helium was used as the carrier gas at a constant flow rate of 2.0 mL/min. The oven temperature program was as follows: initial temperature at 130 °C held for 0.5 min; increased to 180 °C at 10 °C/min and held for 0.5 min; then increased to 225 °C at 5 °C/min and held for 1 min. Mass spectrometric conditions were as follows: ion source temperature, 250 °C; quadrupole temperature, 150 °C; transfer line temperature, 250 °C; electron impact ionization energy, 70 eV; mass scan range, m/z 29–450; solvent delay time, 1.05 min. Fatty acids were identified by comparison with standard mass spectra and retention times, and quantified using an internal standard calibration method.

Intestinal content samples were delivered to Tsingke Biotechnology Co., Ltd. (Guangzhou, China) for high-throughput sequencing of the V3-V4 regions of the 16S rRNA gene. The data processing workflow for 16S rRNA sequencing was as follows: 1) Total DNA was extracted from intestinal contents and conserved-region primers with adapter sequences were used for PCR amplification. Purified, quantified, and normalized amplicons were used to construct sequencing libraries, which were quality-checked and sequenced on the Illumina NovaSeq 6000 platform to obtain raw reads (1.20×10⁷). 2) Raw sequencing data were first filtered using Trimmomatic v0.33 (https://github.com/usadellab/Trimmomatic), and quality-controlled clean reads (1.10×10⁷) were subsequently obtained with Cutadapt v1.9.1. 3) The DADA2 (Louca et al., 2016) pipeline implemented in QIIME2 2020.6 (Bolyen et al., 2019) was used for denoising (with maxEE set to 2, where EE = Σ 10^-Q/10, and all other parameters at default values), merging paired-end reads (minOverlap = 18, maxMismatch = 18 × 0.2), and removing chimeric sequences (consensus model), resulting in high-quality non-chimeric reads (1.05×10⁷) and amplicon sequence variants (ASVs) for downstream analyses. 4) Taxonomic annotation of the ASVs was performed using a naïve Bayes classifier trained on the SILVA reference database. Species abundance tables at different taxonomic levels were then generated using QIIME, and community composition profiles were visualized in R.

Sample genotyping was performed using the “Sea Bass No.1” array via the genotyping-by-target sequencing (GBTS) method. The raw sequencing reads were quality-controlled via Fastp (Chen et al., 2018) software (v0.20.0, parameters: -n 10 -q 20 -u 40) through three filtering steps: 1) removal of adapter sequences; 2) exclusion of paired-end reads containing >10 ambiguous bases (N); and 3) elimination of paired-end reads with low-quality regions (Q ≤ 20) exceeding 40% of the total bases. Cleaned reads were aligned to the spotted sea bass reference genome (GenBank accession: GCA_004028665.1) via SENTIEON BWA software (Jin et al., 2017) in MEM alignment mode. Variant detection was subsequently conducted using the Haplotyper model in SENTIEON DRIVER (parameters: -emit conf 30 –call conf 30 --genotype model multinomial --emit mode GVCF --phasing 0) to generate GVCF files. SNP filtering was implemented through PLINK software (Purcell et al., 2007) with the following criteria: minimum minor allele frequency (MAF) ≥ 0.05 (Chen et al., 2018), missing genotype rate< 20%, heterozygosity rate ≤ 50%, sequencing depth ≥ 5×, and exclusion of non-biallelic loci. The retained loci required successful genotyping in over 80% of analyzed individuals. Finally, population variants were functionally annotated using ANNOVAR software (Wang et al., 2010).

To comprehensively evaluate the genetic variation and population structure of spotted sea bass, this study first performed principal component analysis (PCA) via filtered SNP markers with GCTA software (v1.92.4) (Chen, 2016), visualizing interpopulation genetic differentiation patterns through dimensionality reduction. Bayesian clustering analysis was then conducted via Admixture (v1.3), iteratively testing ancestral population numbers (K = 1–15) and determining the optimal population stratification based on the K value with the lowest cross-validation error (CV error) (Wang et al., 2010). The genetic composition across geographic populations was further quantified using PopHelper (v2.2.7) to generate stacked bar plots reflecting gene flow intensity. Concurrently, genotype data were first converted into a MEGA-format alignment file using a Perl script, with heterozygous sites encoded according to the IUPAC degenerate base standard. MEGA-X was then used to construct a phylogenetic tree based on the Neighbor-joining method and the p-distance model, with 1,000 bootstrap replicates performed to assess the robustness of the inferred topology. The analysis produced a Newick-format tree file containing both branch lengths and node support values. Finally, the resulting tree was visualized in R using the ggtree package. Kinship analysis implemented through GCTA generated genetic relationship matrices, complemented by kinship coefficient frequency distribution plots to infer relatedness among individuals.

Genetic parameters related to growth, gut microbiota, and fatty acid traits were estimated using the GREML approach (Yang et al., 2011) and mixed linear model (MLM) through ASREML (v4.2). The heritability of the above traits was estimated using a univariate animal model, as follows:

Where h² represents heritability, σg2 represents additive genetic variance, and σe2 represents residual variance.

The phenotypic and genetic correlations were estimated using two approaches. First, a bivariate model was fitted with the mmer function from the sommer (v4.3.2) R package. Second, a multi-trait model was implemented using Hiblup software (Linux x86_64 v1.6.0, 2025-09–29 Release), and the standard error of the genetic correlation was calculated using the Delta method.

Where r represents correlation. Cov (x, y) represents the covariance between the two traits, and σx2 and σy2 represents the respective variances of the two traits.

GWAS was conducted using GEMMA software to establish both generalized linear model (GLM) andlinear mixed model (LMM) for analyzing associations between phenotypic traits and genotypes (Ke et al., 2022). The GLM (Equation 1) and LMM (Equation 2) models can be expressed as follows:.

where Y represents the phenotypic vector, X represents the genotype matrix, α represents genotype effect vectors, Q represents the fixed-effect matrix (e.g., population structure, sex (because the sex of the samples was unknown, it could not be included in the model), location, or batch), β represents fixed-effect vectors, K represents the random-effect matrix (primarily the kinship matrix), μ represents random-effect vectors, and e represents the residual vector. The phenotypic variance explained (PVE) by the SNPs was quantified to assess their effect strength using the following formula (Liu, 2023):

where β represents the estimated effect size, MAF represents the minor allele frequency, and N represents the sample size. TASSEL software (Bradbury et al., 2007) was employed to compute two GWAS models: the GLM and LMM. The population structure matrix (Q) derived from the optimal K value identified by admixture analysis and the K generated via GCTA (v1.92.4) were incorporated into the LMM model. The resulting p-values were transformed as -log10 values for visualization in Manhattan and Q-Q plots. A significance threshold of 1×10–⁴ was applied (Liu, 2023), as the sample size (n = 150) and the polygenic trait architecture reduce statistical power, and the Bonferroni correction (P = 0.05/total SNP count) (Ke et al., 2022) would be excessively conservative under these conditions.

Candidate gene mining was performed based on genome-wide significant SNPs identified through GWAS. The core SNP effect regions were defined according to linkage disequilibrium (LD) decay curves, and candidate genes were annotated within 100 kb upstream/downstream intervals of significant SNPs. The genomic coordinates of the candidate genes were preliminarily interpreted using reference genome annotation files. Sequence information was functionally annotated via ANNOVAR software (Wang et al., 2010) to prioritize candidate genes near significant SNPs. To comprehensively characterize gene functions, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted for candidate genes using the OmicShare online platform (https://www.omicshare.com/tools) (Liu, 2023), elucidating their biological roles and pathway associations.

Experimental data are presented as the mean ± standard deviation (Mean ± SD) and graphs were generated using GraphPad Prism software (version 8.0).

In this study, R software (v 3.5.0) was used to detect outliers (based on the 3σ rule, values outside mean ± 3SD) and to conduct statistical analyses on phenotypic traits, gut microbiota, and fatty acid profiles of 150 individuals (Supplementary Table S1). Substantial growth variation was observed, with body weight (BW, 20.25%), snout length (SL, 19.47%), and caudal peduncle length (CPL, 10.30%) showing the highest coefficients of variation (CVs). Gut microbiota exhibited remarkable inter-individual differences, as indicated by extremely high CVs in abundance-based coverage estimator (ACE, 102.80%), Chao1 (103.00%), Feature (OTU/ASV, 103.10%), and the number of Feature (OTU Num, 103.10%). Fatty acid traits exhibited even greater variability, particularly C17:1n-7 (133.30%), C18:1n-9t (201.50%), C18:2n-6t (107.00%), and C22:2n-6 (163.50%). These findings suggest significant breeding potential for both growth performance and fatty acid composition. Normality tests confirmed that most of traits generally followed normal distributions (Figure 2).

Based on high-quality SNP data, we estimated the heritability of growth, gut microbiota, andfatty acid traits (Supplementary Table S2). The results showed that, except for the Shannon index (h²: 0.818 ± 0.164), all measured traits showed moderate-to-low heritability. Phenotypic (r_p) and genetic (r_g) correlation analyses, conducted using the corrplot package in R, revealed strong intra-trait correlations but considerable variants across trait categories (Figure 3; Supplementary Table S2). Overall, most growth traits showed positive correlations (r_p= 0.038–0.994; r_g= 0.047 ± 0.4230 to 0.974 ± 0.044). Within the gut microbiota traits, the number of reads corresponding to each feature (Seqs Num) exhibited negative correlations with other microbiota parameters (r_p= -0.711 to -0.073; r_g= -0.711 ± 0.000 to -0.074 ± 0.003), except for the Simpson index (r_p= 0.007; r_g= 0.003 ± 0.002). For fatty acid traits, both phenotypic and genetic correlations exhibited a diverse pattern (r_p= -0.609 to 0.971; r_g= -0.192 ± 0.019 to 0.971 ± 0.011). Notably, DHA showed relatively high phenotypic and genetic correlation coefficients with EPA (r_p= 0.963; r_g= 0.962 ± 0.000) and C18:1n-9 (r_p= 0.851; r_g= 0.851 ± 0.001). In addition, DHA exhibited relatively strong genetic correlations with ARA (r_g= 0.961 ± 0.001), and C16:0 (r_g= 0.800 ± 0.000). Moreover, the overall pairwise phenotypic and genetic correlations among traits were generally low to moderate in magnitude (r_p= -0.249 to 0.407; r_g= -0.153 ± 0.122 to 0.298 ± 1.033). Overall, these results indicated that co-selecting for growth and fatty acid traits is challenging. Consequently, breeding programs need to adopt specific selection strategies for the simultaneous improvement of multiple traits.

Genotyping was performed on 150 individuals from the same population and developmental stage(Supplementary Table S3), generating 241 Gb of raw data (Supplementary Table S4). After quality control, 229 Gb of clean data were retained, averaging 1.53 Gb per sample. Sequencing evaluation revealed (Figures 4A, B) that read efficiencies ranged from 92.90% to 96.57% (mean: 95.06%) and mapping rates ranged from 95.55% to 97.07% (mean: 96.31%). SNP call rates ranged from 98.99% to 99.90% (mean: 99.66%). Sequencing depth ranged from 79.63× to 174.63× (mean: 118.92×). Quality metrics confirmed high data reliability, with Q20 scores of 99.17%-99.48% (mean: 99.32%) and Q30 scores of 96.39%-97.77% (mean: 97.04%). Bioinformatics processing identified 61,909 mSNPs, and 41,535 high-quality SNPs were retained after stringent filtering for GWAS analysis (Table 1). These SNPs were evenly distributed across 24 chromosomes, with the highest density on chromosome 23 (Figure 5). MAF analysis revealed two major peaks: 36.88% in the 0.4-0.5 range and 31.41% in 0.3-0.4 range (Figure 4D). Genomic annotation indicated that SNPs located in intergenic, upstream regulatory, and downstream regulatory regions accounted for 30.05%, 6.64%, and 5.73% of the total, respectively (Figure 4C).

To assess genetic background variation within the population, multidimensional analyses were performed using 41,535 high-quality SNPs. PCA revealed clear population stratification, with the first three principal components explaining 3.7%, 2.2%, and 1.68% of the variance, respectively. Although most samples clustered in the lower-right quadrant, a small subset showed separation along PC1 (Figure 6A). The phylogenetic tree results supported this structure by separating individuals into several evolutionary branches (Figure 6B). These patterns were further corroborated by CV error trends (Figure 6C) and population admixture analysis (Figure 6D), collectively demonstrating stable genetic heterogeneity within the cohort. The kinship heatmap showed low overall relatedness, with most values ranging from −0.05 to 0.00 (Supplementary Figure S1). Concordance across multiple analytical methods provide robust evidence of population-level genetic divergence, with a consistent tripartite substructure characterizing this breeding population. To minimize potential effects of population stratification, the first three PCs were included as covariates in subsequent analytical models.

LD analysis was performed using PopLDdecay, with systematic estimation of r² values between genetic markers. As shown in Supplementary Figure S2, the LD decay curve exhibited a typical exponential pattern, where r² decreasing rapidly as physical distance increased. LD intensity stabilized at approximately 0.05 when inter-marker distances reached ~100 kb. This rapid LD decay (r²< 0.1 within 100 kb) reflected the historical effective population size, suggesting extensive recombination and relatively low selection pressure. The observed LD patterns offer essential guidance for optimizing marker density in subsequent association mapping studies (Tao et al., 2026).

For the GWAS of growth traits, 31 significant SNPs were identified, with detailed associations shown in Figure 7; Supplementary Figure S3 and Supplementary Table S5. Notably, the locus Superscaffold2_21167501 exhibited pleiotropic effects, associating simultaneously with BW, TL, BL, BH, and caudal peduncle height (CPH). Although these SNPs displayed trait-specific effects, their PVE remained moderate generally (9.89%-14.96%).

For the GWAS of microbiome traits, 35 significant SNPs were identified (Supplementary Table S6; Supplementary Figure S4). Seven loci (Superscaffold1_22134058, 2_16860820, 5_25307233, 6_15942598, 14_2978392, 14_12467889, and 16_10477403) were shared across ACE, Chao1, Feature, and OTU Num. The pleiotropic locus 16_24458793 was associated with ACE, Chao1, Feature, OTU Num, and Species, whereas 19_3679037 showed dual associations with Genus and PD whole tree (Phylogeny-based diversity indices) metrics. These shared markers exhibited moderate PVE values (9.79%-21.49%).

For the GWAS of fatty acid traits, a total of 124 significant SNPs were identified (Supplementary Table S7; Supplementary Figure S5). Several loci were shared across traits: Superscaffold9_11794470 for C14:0 and C18:0; 2_5033478 and 4_14696439 for C15:0 and C18:3n-6; 18_22390868 and 24_19181353 for C15:0 and C16:0; 23_21849556 for C15:0, C16:0, C17:0, C18:1n-9, C18:2n-6, and ALA; 12_4089269 for C18:1n-9, C18:3n-6, EPA, DHA, and DHA/ALA; and Superscaffold17_8900366 for DHA with DHA/ALA. Shared SNPS exhibited relatively low PVE values (9.68%-16.33%).

Collectively, these findings indicate that these traits are controlled by multiple loci and exhibit pleiotropy. The consistently low PVE values further support their polygenic architecture.

Based on LD decay results, r² stabilized at 100 kb, and therefore genomic intervals extending 100 kb upstream and downstream of significant markers were used for candidate gene identification and annotation (Supplementary Figure S2).

A total of 37 candidate genes associated with growth traits were identified (Supplementary Table S8). GO annotation (Supplementary Figure S6) highlighted three core functional modules: signal transduction, metabolic regulation, and developmental regulation. Key molecular functions included CDP-diacylglycerol-glycerol-3-phosphate 3-phosphatidyltransferase activity, oxidoreductase activity (Peptidyl-proline dioxygenase) and antioxidant activity. KEGG enrichment (Figure 8) revealed growth-associated pathways such as Th1 and Th2 cell differentiation/PI3K-Akt signaling, glycosaminoglycan/N-glycan biosynthesis, and Notch signaling, forming a metabolism-development integration network.

For gut microbiota, 40 candidate genes were identified (Supplementary Table S9). GO enrichment (Supplementary Figure S7) indicated primary roles in metabolic regulation, neurodevelopment, and proteostasis. Molecular functions were enriched in kinase activities and nucleotide metabolism. KEGG analysis (Supplementary Figure S8) revealed an integrated immune-metabolic-neurodevelopmental network, with key pathways including mitophagy, folate-mediated one-carbon metabolism, and axon guidance. Enrichment of the bacterial invasion of epithelial cells suggested potential host-microbiota interactions. These findings highlight neuro-endocrine-immune regulation as a major axis shaping gut microbial traits.

For fatty acid traits, 148 candidate genes were identified (Supplementary Table S10). GO analysis (Supplementary Figure S9) revealed predominant involvement in intracellular signaling, metabolic processes, and kinase activity. Molecular functions included cytoskeleton construction, oxidative nucleic acid demethylation, and kinase regulation. KEGG enrichment (Supplementary Figure S10) revealed several core pathways encompassing endoplasmic reticulum protein processing, MAPK/Rap1 signaling, stress response, material and energy metabolism, cell proliferation/apoptosis, and neural signaling. These results reveal a complex interaction network integrating energy homeostasis, stress adaptation, and cell cycle regulation underlying lipid traits.

Integrated multi-omics analysis revealed shared regulatory pathways, while maintaining trait-specific gene combinations. Candidate genes collectively shaped phenotypic variation through coordinated metabolic and signaling cascades. Notably, neurodevelopmental pathways were enriched in microbiota-related genes, expanding the recognized regulatory scope of host genetics over the gut microbiome. Overall, these findings offer multidimensional insights into the molecular mechanisms underlying complex traits in spotted sea bass and establish a theoretical framework for molecular breeding.