Predominantly Independent Genetic Control Between Growth and Visceral White Nodules Disease Resistance Revealed by High-Density Linkage Map and QTL Mapping in Larimichthys crocea

TL;DR

This study finds that growth and disease resistance in large yellow croaker are largely controlled by separate genes, offering hope for improving both traits in breeding.

Contribution

A high-density genetic map and QTL analysis reveal independent genetic control of growth and disease resistance in Larimichthys crocea.

Findings

53 QTLs for growth and 20 for VWND resistance were identified with limited overlap.

Only two genomic intervals showed co-localization between growth and resistance traits.

Candidate genes in co-localized regions are involved in general cellular processes.

Abstract

The large yellow croaker (Larimichthys crocea) is a key mariculture species in China, however, its industry is threatened by visceral white nodules disease (VWND) caused by the bacterium Pseudomonas plecoglossicida. A significant challenge in breeding is the potential genetic trade-off between growth and disease resistance. To investigate their genetic relationship, we constructed a high-density SNP-based genetic linkage map for L. crocea using a F1 full-sib family (n = 150). The map comprised 24 linkage groups with 32,429 bin markers and an average interval of 0.051 cM. Based on this map, we conducted QTL mapping for one yield trait (body weight), eight morphological traits, and three VWND-resistance traits (survival time, AT; spleen and liver pathogen loads). Phenotypic analysis revealed strong integration among growth traits and a moderate positive correlation between growth traits and AT. QTL mapping identified 53 QTLs for growth (PVE = 0.14–5.83%) and 20 for resistance (PVE = 0.78–8.93%). Notably, only two genomic intervals exhibited co-localization between a morphological trait (AL or BL) and AT, each explaining a modest phenotypic variance (0.66–5.99%). The largest-effect QTLs for growth and resistance were mapped to distinct linkage groups, and candidate genes within the co-localized intervals (Unc5d, SCN5A, HUS1) are involved in fundamental cellular processes rather than core growth or immune pathways. These results suggest that yield, morphological, and VWND-resistance traits in L. crocea are largely under independent genetic control within the studied family, indicating that simultaneous improvement of growth and disease resistance is feasible. This study provides a molecular basis for breeding strategies aimed at overcoming the trait trade-off bottleneck in this economically vital species.

Figures

Full text

1. Introduction

The large yellow croaker (Larimichthys crocea) is a core mariculture fish species along the southeast coast of China, with an annual output of 257,000 tons nationwide in 2022, occupying an important position in the aquatic product industry [1]. In intensive aquaculture systems, visceral white nodules disease (VWND), caused by Pseudomonas plecoglossicida, poses a prominent threat, with a mortality rate of up to over 80%, resulting in substantial economic losses to major producing areas [2]. Unfortunately, pharmacological treatments have limited efficacy once infection occurs. Consequently, breeding disease-resistant strains has emerged as a safer and more sustainable alternative. Growth performance (including yield trait and morphological traits) and disease resistance are core breeding traits in L. crocea, and their genetic interaction (antagonism or synergy) is a critical issue for aquatic breeders. This contradiction is pervasive across animal and plant breeding, with extensive studies confirming the common existence of genetic trade-offs or synergies between growth performance and bacterial disease resistance. In aquatic animals, such trait interactions have been widely documented, reflecting inherent evolutionary and physiological constraints in genetic improvement [3]. Similarly, crop breeding research has consistently reported genetic interplay between growth and resistance, often manifesting as trade-offs that impede simultaneous enhancement of both traits [4,5]. The trade-off between these two traits has become an industrial bottleneck: selective breeding of disease-resistant strains alone tends to cause growth retardation, whereas excessive pursuit of rapid growth reduces population disease resistance. Clarifying the genetic correlation mechanism between them is thus the core premise for breaking this bottleneck.

Conventional approaches for genetic correlation analysis in aquatic animals mainly involve morphological phenotypic selection, family selection, population phenotypic correlation analysis, and quantitative genetic evaluation [6,7,8]. While effective for species with short generation cycles and high environmental controllability, these methods are particularly challenging for L. crocea due to its long generation cycle (2~3 years), high breeding costs, and phenotype plasticity in response to environmental factors such as temperature and diet. Consequently, traditional methods are not only protracted and costly but may also yield biased estimates of genetic correlations due to unresolved genotype-by-environment interactions, hindering the accurate localization of causative loci [9,10]. The advent of molecular markers has enabled more direct genetic analyses. Genome-wide association study (GWAS) has become a predominant tool in L. crocea research, largely circumventing the need for controlled familial construction [11,12,13]. However, robust GWAS requires large and diverse populations, and its correlative nature can struggle to distinguish between tight linkage and pleiotropy—a critical distinction for understanding the genetic relationship between complex traits like growth and disease resistance [1,14]. Linkage analysis based on high-density genetic mapping provides a complementary and powerful strategy for dissecting such trait correlations in structured pedigrees. It can accurately determine the genetic patterns (antagonism, synergy, or independence) between traits through locus overlap, linkage, or segregation, and has been widely applied in aquatic breeding [15,16,17].

Significant progress has been made in constructing genetic maps for L. crocea, laying a solid foundation for subsequent genetic research. However, inherent technical limitations of earlier studies have constrained in-depth investigations into trait correlations. Early AFLP maps developed by Ning et al. [18] and microsatellite maps by Ye et al. [19], with an average interval of 5.4 cM, provided valuable insights for preliminary QTL mapping of growth traits, though their relatively low marker density precluded fine-scale analysis of genetic interactions. The first SNP map constructed by Ao et al. [20] represented a notable advancement, reducing the average interval to 0.54 cM, yet it still faced challenges of marker segregation distortion and incomplete genome coverage. Furthermore, existing studies have focused on single traits either solely analyzing growth [11,19] or targeting a single disease [21,22] without investigating the correlation between growth and visceral white spot disease resistance, resulting in a lack of molecular basis for breeding practice.

Targeting the core breeding challenge of genetic interaction between growth performance and VWND-resistance, this study aimed to (1) construct a high-density SNP-based genetic linkage map for L. crocea using a full-sib F1 family and whole-genome resequencing, (2) perform comprehensive QTL mapping for one yield trait, eight morphological traits, and three VWND-resistance traits, and (3) elucidate the genetic relationship between growth and disease resistance through QTL co-localization and trait correlation analysis. It should be noted that the use of a single F1 family for QTL detection is a well-established approach in aquaculture genetic studies, as demonstrated in various species including Luciobarbus brachycephalus [17], L. crocea [22], and Larimichthys polyactis [23]. While multi-family designs can capture broader genetic diversity, single-family designs offer the advantage of controlled genetic background and reduced environmental noise, enabling precise QTL localization. In this study, although five families were initially generated, only one family retained an adequate sample size for robust QTL analysis due to the inherent challenges of long-generation-cycle aquaculture species, including environmental stressors during extended grow-out periods. Therefore, this study represents an initial step toward understanding the genetic architecture of growth and disease resistance in L. crocea, providing a foundation for future validation across broader genetic backgrounds.

2. Results

2.1. Phenotypic Variation and Normality Analysis of Growth and VWND Resistance Traits

Phenotypic determination was performed on 150 F1 individuals of L. crocea, including 9 growth performance traits, yield trait (Wt), morphological traits (AL, BL, CH, CL, TH, TL, QL, WL), and 3 VWND resistance-related traits (AT, PPLL, PPSL). The phenotypic distribution characteristics and normality test results of all traits are shown in Figure 1 and Table S1, respectively. Figure 1 shows the frequency distribution of phenotypic values of all 12 traits. It can be observed that most morphological traits (AL, BL, CH, CL, TH, TL, WL) and one resistance trait (AT) exhibited an approximate normal distribution, with concentrated phenotypic values and uniform distribution trends. In contrast, the yield trait (Wt) and the other two pathogen load traits (PPLL, PPSL) displayed slight skewness, with phenotypic values skewed to the low-value region, indicating that most individuals had relatively low yield potential, liver pathogen load, and spleen pathogen load.

The normality test results (Table S1) further confirmed the distribution characteristics of the phenotypic data. The skewness values of AL, BL, CH, CL, TH, TL, WL, and AT ranged from −0.30 to 0.19, and the kurtosis values ranged from −0.80 to 0.28, all close to 0, indicating that these traits conformed to a normal distribution (p > 0.05). For Wt, PPLL, and PPSL, the skewness values were 0.49, 0.97, and 1.11, respectively, and the kurtosis values were −0.48, 0.28, and 0.16, respectively, showing mild to moderate positive skewness, which might be related to the genetic variation of these traits. Overall, all phenotypic data showed good genetic variation, with the coefficient of variation ranging from 8.92% (TH) to 87.50% (PPSL), providing a solid foundation for subsequent genetic correlation analysis and QTL mapping.

2.2. Phenotypic Correlation Among Yield, Morphological and VWND Resistance Traits

Pearson correlation analysis was performed to assess the relationships among one yield trait, eight morphological traits, and three disease resistance traits in 150 individuals (Figure 2). A strong, significant positive correlation (p < 0.01) was observed among all morphological traits and between the yield trait (Wt) and each morphological trait, with correlation coefficients (r) ranging from 0.73 to 0.99 (Figure 2). Morphological traits AL and BL showed the highest correlation (r = 0.99). The yield trait Wt was highly correlated with BL and AL (r = 0.97), suggesting potential shared genetic regulation of growth and morphology. Among the resistance traits, survival time after challenge (AT) was significantly positively correlated with spleen pathogen load (PPSL) (r = 0.67, p < 0.01). In contrast, liver pathogen load (PPLL) showed only a weak, non-significant correlation with PPSL (r = 0.16) and no correlation with AT (r = 0.03), indicating distinct regulatory pathways for different resistance components.

Notably, both the yield trait (Wt) and all morphological traits showed positive correlations with survival trait AT (r = 0.45–0.52, p < 0.01) and showed weak positive correlations with PPSL (r = 0.21–0.27, p < 0.05). However, no significant correlations were detected between yield/morphological traits and hepatic pathogen load (PPLL) (r = −0.07–0.03, p > 0.05). These results indicate that individuals with superior growth and morphology exhibited longer survival and moderately higher splenic pathogen burden following challenge, whereas hepatic pathogen load was independent of growth and morphological phenotype.

2.3. Linkage Analysis and Genetic Map Construction

A high-resolution genetic linkage map was constructed for L. crocea based on whole-genome resequencing of a full-sib F1 family, following the pseudo-testcross strategy and using filtered bin markers. The resulting consensus map consisted of 24 linkage groups (LGs) (Figure 3) corresponding to the 24 chromosomes of the species (Table 1). A total of 32,429 bin markers were integrated, spanning a total genetic length of 1417.17 cM, with an average inter-marker distance of 0.051 cM. The linkage groups ranged in genetic length from 45.72 cM (LG14) to 82.43 cM (LG18), and the number of mapped markers per LG varied from 377 (LG18) to 1862 (LG4). The average marker interval across LGs ranged from 0.033 cM (LG1, LG4, LG15, LG16) to 0.219 cm (LG18), indicating high marker density and uniform distribution. More than 99.5% of adjacent marker pairs were separated by less than 1 cM, confirming minimal gaps and comprehensive genome coverage.

Quality Assessment and Collinearity of the Genetic Map

To evaluate the quality and accuracy of the constructed genetic map, we further analyzed the collinearity between the genetic map and the reference genome of L. crocea. By comparing the physical positions of markers on the genome with their corresponding genetic distances on the linkage map, collinearity analysis revealed a strong overall concordance between the two datasets (Figure 4). The marker order within most linkage groups was consistent with the physical order in the genome, with no large-scale inversions or rearrangements observed. Minor local discrepancies in marker order may be attributed to genomic assembly gaps or recombination hotspots. These results confirm the high accuracy and reliability of the high-density genetic map constructed in this study, providing a robust framework for subsequent fine-scale QTL mapping.

2.4. QTL Mapping for Growth-Related Traits

QTL mapping for growth-related traits identified a total of 53 significant QTLs (LOD ≥ 3.0) distributed across 16 linkage groups (LGs), pertaining to the yield trait (Wt) and three key morphological traits: AL, BL, and CH (Figure 5). Specifically, Wt was associated with 12 QTLs, each explaining 0.14 to 5.83% of the phenotypic variance (PVE) (Figure 5C). Notably, five of the Wt-related QTLs co-localized with QTLs for AL and/or BL, particularly on LG6, LG10, LG16, and LG23, suggesting shared genetic regulation between weight gain and body elongation. For the length-related traits, AL and BL shared 8 co-localized QTLs, located on LG6 (2 QTLs), LG10 (2 QTLs), LG16 (2 QTLs), LG20, and LG23 (Figure 5A,B). Among these, the QTL on LG6 (20.8–22.1 Mb) showed the strongest effects on both AL (LOD = 5.75, PVE = 2.56%) and BL (LOD = 6.38, PVE = 1.49%), representing a major locus for longitudinal growth (Table 2). CH was regulated by 14 QTLs distributed across eight LGs (Figure 5D). Four of these co-localized with QTLs for AL and BL, all located on LG6 and LG16, indicating genetic regions that influence both body height and length. The most prominent co-localized region was on LG6 (20.4–28.4 Mb), which contained high-effect QTLs for all four traits (Wt, AL, BL, CH), with LOD scores ranging from 3.24 to 6.98 (Table 2).

Furthermore, analysis of the remaining morphological traits, head CL, QL, TH, TL, and WL revealed that QTLs for these traits also co-localized within the same intervals on LG6 and LG16. Notably, the region on LG16 (15.84–16.23 Mb) contained overlapping QTLs for all nine growth-related traits examined, representing a genomic hotspot with pleiotropic effects on overall body morphology. Overall, LG6 and LG16 were identified as genomic hotspots harboring multiple co-localized QTLs for growth-related traits (Table 2 and Table S2). In particular, LG6 contained a cluster of QTLs with major effects on weight, length, and height, suggesting the presence of a key genomic region governing overall growth performance and body conformation in L. crocea. QTL mapping results for the remaining morphological traits are provided in Supplementary Figure S1.

2.5. QTL Mapping for VWND-Resistance Traits

A total of 20 significant QTLs were identified for the three disease resistance traits: survival time after challenge (AT), relative pathogen load in spleen (PPSL), and relative pathogen load in liver (PPLL) (Figure 6). AT was associated with five QTLs distributed across LG5, LG9, LG16, LG22, and LG23 (Figure 6A), with PVE values ranging from 2.69% to 8.93%. The QTL on LG9 (32.11–32.21 Mb) showed the strongest effect (PVE = 8.93%, LOD = 3.52) (Table 3). For pathogen load traits, eight QTLs were detected for PPSL on LG1, LG2, LG6, LG9, LG14, LG16, and LG22, explaining 1.80–5.58% of phenotypic variance (Figure 6B). The QTLs on LG1 (20.48–20.67 Mb) and LG2 (21.62–22.79 Mb) exhibited the highest PVE values (5.58% and 5.52%, respectively). Seven QTLs were associated with PPLL, located on LG4, LG6, LG9, LG21, and LG22, with PVE ranging from 0.78% to 4.87% (Figure 6C). The most significant PPLL QTL was mapped to LG9 (27.18–27.28 Mb; LOD = 4.04, PVE = 4.19%). Notably, several genomic regions exhibited the co-localization of QTLs for multiple resistance traits. LG9 contained QTLs for all three traits within three distinct intervals: 16.85–17.02 Mb (affecting PPSL), 21.30–27.27 Mb (affecting PPLL), and 32.11–32.21 Mb (affecting AT). On LG16, QTLs for AT and PPSL were clustered within three adjacent intervals (20.17–20.25 Mb, 21.26–21.41 Mb, and 23.88–25.21 Mb). On LG22, overlapping QTLs for AT, PPLL, and PPSL were identified within the region 0.34–2.62 Mb (Table 3 and Table S2). In contrast, no overlapping QTLs were observed between PPLL and AT in other genomic regions, indicating that liver pathogen load may be regulated independently of survival duration under challenge conditions.

2.6. Co-Localization of Growth and Disease Resistance

To investigate the genetic relationship between growth performance and VWND resistance, we analyzed the co-localization of QTLs identified for both trait categories. Two genomic regions were found to harbor co-localizing QTLs influencing growth and disease resistance traits (Table 4). On LG16, a QTL interval spanning 21.36–21.41 Mb was shared between AL and AT, explaining 3.37% and 3.57% of PVE, respectively. Within this interval, a single annotated gene, Unc5d, was identified. Unc5d encodes a netrin receptor involved in axon guidance and has been implicated in immune modulation and cellular stress responses in vertebrates. On LG22, a QTL region between 1.48 and 1.68 Mb co-localized for BL and AT. This locus accounted for 0.66% of phenotypic variance in BL but exerted a stronger effect on AT, explaining 5.99% of variance. Two genes, SCN5A and HUS1, were annotated within this region. SCN5A encodes a voltage-gated sodium channel subunit associated with cardiac electrophysiology and stress adaptation, while HUS1 is a component of the DNA damage checkpoint complex involved in cellular senescence and genomic integrity maintenance.

3. Discussion

This study constructed a high-density SNP-based genetic map for L. crocea and performed systematic QTL mapping for two economically vital growth trait categories: yield (body weight, Wt) and morphological traits (e.g., AL, BL, CH) as well as resistance traits against visceral white nodules disease (VWND). The high-resolution map enabled precise dissection of the genetic architecture underlying these traits. Overall, our results revealed strong genetic integration within growth traits, partial independence among resistance components, and largely separate genetic control between growth and disease resistance, with limited QTL co-localization. These findings provide important insights for breeding programs aiming to simultaneously improve yield, morphology, and disease resistance in L. crocea.

3.1. Phenotypic Correlations Reflect Integrated Growth Traits and Distinct Components of Disease Resistance

The traits analyzed in this study—nine growth-related and three disease resistance-related traits—all showed continuous variation, fitting their quantitative genetic nature. Strong positive correlations were observed among all morphological traits and between these traits and the yield trait (Wt), aligning with findings in other fish species [6,24]. This integration likely reflects the coordinated development program governing body size prior to sexual maturity. Among the disease resistance traits, survival time (AT) was positively correlated with splenic pathogen load (PPSL), consistent with our earlier report [21]. AT is a widely adopted and stable phenotypic indicator that integrates the host’s overall physiological and immune response [25,26]. In contrast, the weak correlation between liver (PPLL) and spleen (PPSL) pathogen loads suggests organ-specific infection dynamics, possibly due to differences in local immunity, metabolism, or tissue tropism of P. plecoglossicida [27].

Notably, the yield trait (Wt) and all morphological traits showed moderate positive correlations with AT and weak links to PPSL, but no significant association with PPLL. This is consistent with previous findings in Mekong striped catfish (Pangasianodon hypophthalmus), where growth performance was reported to have a moderate favorable genetic correlation with survival under pathogen challenge while showing weak correlations with pathogen susceptibility-related traits [28]. This indicates that faster-growing individuals may survive longer under challenge and tolerate slightly higher splenic pathogen loads, without a corresponding rise in hepatic infection. Similar organ-specific pathogen load patterns have been observed in various aquatic and terrestrial hosts, where splenic pathogen accumulation can be dissociated from hepatic infection during the progression of disease, reflecting differential organ-specific immune responses to pathogens [29,30]. In farming environments, individuals with higher growth rates may possess better overall vigor or stress tolerance, contributing to extended survival despite infection. This phenomenon has been widely documented in aquaculture species, including Labeo rohita [31] and Oreochromis niloticus [32], where enhanced growth performance is often associated with improved disease resistance and stress resilience, either through intrinsic physiological vigor or better resource allocation for immune defense [33].

3.2. A High-Density SNP Genetic Map Enables Accurate QTL Mapping in L. crocea

The genetic map constructed in the present study had an average marker interval of 0.051 cM, representing one of the densest linkage maps developed for a marine fish [20,22,23]. High heterozygosity and environmental sensitivity often constrain genetic studies in aquaculture species, as these factors can lead to unstable phenotypic measurements and reduced accuracy of genetic linkage analysis [34,35]. The success of this map demonstrates the utility of whole-genome resequencing combined with a bin-marker strategy, an approach applicable to other aquatic organisms with complex genomes. This strategy has been successfully employed in other marine and freshwater fish species, such as Mylopharyngodon piceus [36] and L. polyactis [23], to overcome genome complexity and construct high-quality genetic maps.

High map density directly improved QTL detection, as dense marker coverage minimizes the genetic distance between markers and target QTLs, thereby reducing the probability of marker-QTL recombination and enhancing the efficiency of QTL localization [37]. All 12 traits yielded significant QTLs, with PVE ranging from 0.14% to 8.93%. The wide range of PVE values suggests that these traits are controlled by multiple minor-effect QTLs, which is consistent with the quantitative genetic characteristics of most economic traits in aquaculture species [38]. Similar findings have been reported in previous QTL studies on aquatic economic traits. Ali et al. [39] detected growth-related QTLs with relatively low PVEs in rainbow trout, and Liu et al. [40] reported analogous results (PVEs ranging from 2.2% to 4.1%) for disease resistance traits in Asian seabass, both indicating that these quantitative traits are regulated by multiple minor effect QTLs, which is consistent with our study results. The efficient detection of these minor effect QTLs in our study, even with PVEs as low as 0.14%, further highlights the advantage of our high-density genetic map combined with precise phenotyping in dissecting the genetic basis of complex quantitative traits in aquaculture species [36].

3.3. QTL Co-Localization Indicates Largely Independent Genetic Control of Growth and Disease Resistance

A key finding of this study is the general genetic independence between growth traits (both yield and morphology) and VWND resistance. Among the 53 growth-related and 20 resistance-related QTLs we identified, only two genomic intervals exhibited co-localization between morphological traits (AL or BL) and AT. Crucially, the QTLs with the largest effects on growth (e.g., on LG6 and LG16) and those with the largest effects on resistance (e.g., on LG9 and LG22) mapped to distinct linkage groups, indicating that the major genetic determinants for these trait complexes are separate.

Furthermore, this pattern of limited genetic overlap is not unique to L. crocea. In agreement with our findings, studies in other aquaculture species have reported similar dissociations. For instance, in Litopenaeus vannamei, genetic correlations between growth traits and resistance to white spot syndrome virus are typically low, enabling independent selection for these traits [41,42]. Similarly, in Salmo salar, genome-wide association studies have identified largely distinct genomic regions controlling growth and resistance to sea lice or pancreas disease, with minimal locus overlap [43]. Moreover, we observed that the genetic pattern aligned well with the phenotypic correlation results. In this study, the moderate correlation between growth and AT may be attributed to a limited number of shared quantitative trait loci (QTLs). Conversely, the lack of correlation with PPLL/PPSL suggests that there are no overlapping QTLs between the two traits. Furthermore, even when the QTLs associated with growth and resistance are located on the same linkage group, they are often situated at considerable physical distances from one another, indicating a high probability of recombination. This further supports the notion of relatively independent regulatory mechanisms governing these traits. Similar genetic mechanisms have been extensively documented in plants [44,45]. In wheat and other cereal crops, overlapping QTLs governing growth and disease resistance traits are exceedingly rare [46]. Even when they are found within the same linkage group, the physical distance between the QTL intervals is typically substantial.

3.4. Candidate Genes in Co-Localized Regions Support a Model of Weak Pleiotropy and Predominantly Independent Genetic Control

From a breeding perspective, this largely independent genetic control indicates that simultaneous improvement of yield, morphology, and VWND resistance is feasible. However, it is still necessary to closely monitor a few co-localized regions (such as LG16 and LG22) to manage potential trait interactions during the selection process. In the present study, the annotated candidate genes within these regions, Unc5d on LG16 and SCN5A and HUS1 on LG22, are not canonical regulators of core growth (e.g., components of the somatotropic axis) or canonical immune pathways (e.g., pattern recognition receptors or inflammatory cytokines) [47,48]. Instead, their known functions are associated with fundamental cellular and systemic processes: neuronal guidance and cellular stress response (Unc5d) [49], electrochemical signaling and systemic stress adaptation (SCN5A) [50], and genomic integrity maintenance (HUS1) [51].

The phenotypic variance explained (PVE) by these co-localizing QTLs further indicates their modest and distinct contributions. The locus on LG16 exhibited comparable effects on total length (AL, PVE = 3.37%) and survival time (AT, PVE = 3.57%). In contrast, the locus on LG22 had a minimal impact on body length (BL, PVE = 0.66%) but a more pronounced effect on survival time (AT, PVE = 5.99%). This pattern is consistent with a model of weak or trait-specific pleiotropy, where genetic variants exert minor, potentially indirect, effects on multiple traits, or where linked genes with distinct primary functions are localized within the same genomic interval [52]. Moreover, the absence of major growth or immune effector genes within co-localization intervals suggests that the observed phenotypic correlations are unlikely to be driven by strong, direct pleiotropy of key developmental or immunological pathways. This finding aligns with reports in other aquaculture species where growth and disease resistance QTLs largely map to distinct genomic regions [41,42,43]. Instead, the shared genetic influence may stem from variants in genes governing systemic homeostasis, basal stress responses, or fundamental cellular functions, which secondarily modulate both growth and defense phenotypes [53]. Collectively, these interpretations reinforce the genome-wide pattern observed here: the vast majority of QTLs for yield, morphology, and disease resistance are genetically distinct. Therefore, the functional annotation and modest effect sizes of the candidate genes within the co-localized QTLs corroborate the conclusion that growth and VWND resistance in L. crocea are primarily under independent genetic control. The identified loci likely represent peripheral, integrative nodes rather than central genetic constraints, a genetic architecture that is favorable for breeding programs aiming to achieve simultaneous genetic gain in both trait complexes without encountering severe antagonistic selection responses [54].

3.5. Study Limitations and Future Perspectives

This study represents an important step toward understanding the genetic relationship between growth and disease resistance in L. crocea through high-density linkage mapping and QTL analysis. However, certain considerations should be taken into account when extrapolating these findings. First, as a single-family study, the QTLs detected herein represent those segregating in this specific family and may not capture all genetic variants influencing growth and disease resistance in the broader population. This is consistent with findings in other species where QTL effects can vary across genetic backgrounds [16]. Nevertheless, single-family QTL mapping remains a valuable approach for initial genetic dissection, and the consistency of our findings with phenotypic correlations supports their biological relevance. Second, most QTLs identified in this study explained relatively low proportions of phenotypic variance (PVE < 5%), consistent with the polygenic architecture typical of complex traits in aquaculture species. Nevertheless, the modest sample size may limit statistical power to detect QTLs with very small effects, and some genuine loci may have escaped detection. Third, although permutation tests were employed to control genome-wide false positive rates-a widely accepted approach in QTL mapping studies [55]. We acknowledge that additional false discovery rate (FDR) control was not applied. Permutation tests were chosen because they account for the correlation structure among linked markers, which FDR procedures may not fully capture [56], and this approach is widely adopted in high-density QTL mapping studies in aquaculture species [17,22,23]. Nevertheless, FDR control could provide complementary information, particularly for prioritizing candidate loci for future validation. Fourth, while candidate genes were identified within co-localized intervals based on genomic annotation, their functional roles in growth and disease resistance remain to be experimentally validated. Future studies incorporating transcriptomic, proteomic, or gene editing approaches will be essential to confirm their biological relevance and facilitate the development of marker-assisted selection programs for simultaneous improvement of growth and disease resistance in L. crocea breeding.

4. Materials and Methods

4.1. Experimental Materials and Full-Sib Family Construction

The experimental fish used in this study were derived from the Daiqu strain of L. crocea. Five full-sibling families were first established via artificial insemination at Xiangshan Harbor Aquatic Hatchery Co., Ltd. (Ningbo, China) in April 2023, following the method described by Ye et al. [21]. Briefly, 10 sexually mature females and 5 males (18 months old, hybrids of wild and farmed populations) were injected with luteinizing hormone-releasing hormone A3 (LHRH-A3) (Ningbo second hormone factory, Ningbo, China) at doses of 1 μg/fish and 0.5 μg/fish (dissolved in 1 mL sterile saline), respectively, for spawning induction. After artificial insemination at a gamete ratio of 1:1, five full-sibling lineages were successfully obtained. A family was selected for the present study. A total of 150 healthy individuals from this family at 8 months of age (average body weight: 22.1 ± 4.2 g) were selected for the bacterial challenge experiment.

4.2. Phenotypic Trait Determination and Data Quality Control

4.2.1. Determination of Yield Trait and Morphological Traits

Nine growth-related traits were determined for each individual, including 1 yield trait and 8 morphological traits: yield trait was body weight (Wt, directly reflecting aquaculture yield); morphological traits included total length (AL), body length (BL), head length (CL), body height (CH), tail height (TH), trunk length (TL), caudal peduncle length (QL), and caudal fin length (WL). These parameters were measured using an automated phenotyping acquisition system as described by Wang et al. [57]. Each trait was independently measured three times by the same operator, and the average value was taken as the final phenotypic value to reduce measurement errors.

4.2.2. Determination of VWND Resistance-Related Traits

The bacterium (P. plecoglossicida) used in this study was isolated from L. crocea affected by a natural outbreak in Xiangshan. Based on 16S rRNA gene (16S) homology and biochemical tests, it was classified as the XSDHY-P strain [58]. The bacterial challenge experiment was performed following the method described by Ye et al. [21] with identical operational procedures. Briefly, bacteria were cultured overnight in TSA liquid medium at 28 °C with 200 rpm shaking for 24 h. The infection concentration was set to 1 × 10^5^ CFU/mL, based on the 96 h LC_50_ previously verified by our team (unpublished data), to ensure high infection success rate and uniformity. A total of 150 experimental fish were intraperitoneally injected with 0.5 mL of the bacterial suspension at the fin base, and the entire injection process was completed within 30 min with water temperature stably maintained. A control group (n = 30) receiving 0.5 mL sterile phosphate-buffered saline (PBS) was included in parallel, and no mortality was observed in the control group throughout the experiment.

Three VWND resistance-related traits were determined: survival time (AT), the relative abundance of P. plecoglossicida in spleen (PPSL), and the relative abundance of P. plecoglossicida in liver (PPLL). Immediately post-injection, all experimental fish were placed under continuous observation. Fish exhibiting a moribund state (defined as complete unresponsiveness to external stimuli) were promptly removed, and the time of removal was recorded as individual survival time (AT). Subsequently, spleen and liver tissues were aseptically dissected on ice, immediately flash-frozen in liquid nitrogen for subsequent pathogen load quantification (PPLL and PPSL), and a fin clip was collected and preserved in absolute ethanol. Observation continued until survival time data had been recorded for all 150 individuals, at which point the experiment was terminated. The identification methods for PPSL and PPLL followed our previously established protocol [21]. Briefly, total DNA was extracted from the entire spleen and liver tissues using the standard phenol-chloroform method. After quality verification by 1% agarose gel electrophoresis, DNA concentration was normalized to 50 ng/μL using a Nanodrop One spectrophotometer (Thermo Fisher Scientific, China). Quantitative PCR (qPCR) was performed with species-specific primers for P. plecoglossicida (pps-F: GGTACTCACCTGGTTGGGCTT; pps-R: GCTTGTCCTTGGTCTGGGAG) and the large yellow croaker single-copy gene gdf8 (gdf8-F: TGGGAGATGACAACAGG; gdf8-R: GCACCGACCAATACACT) as an internal reference. The thermal cycling conditions comprised an initial denaturation at 95 °C for 1 min, followed by 40 cycles of 95 °C for 5 s, 60 °C for 10 s, and 72 °C for 15 s. Relative pathogen load was calculated using the 2^−^ΔΔCt method to quantify P. plecoglossicida proliferation per unit of host tissue. Environmental conditions during the experiment were strictly controlled (water temperature 17~18 °C, salinity fluctuation ± 1‰, dissolved oxygen ≥ 6 mg/L) to ensure consistent challenge pressure.

4.3. SNP Genotyping and Bin Marker Filter

Genomic DNA was extracted from tissue (fin, spleen, liver) using a commercial DNA extraction kit (Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. DNA integrity was verified by 1% agarose gel electrophoresis, purity was detected by Nanodrop 2000 (NanoDrop, Waltham, MA, USA), and concentration was measured by Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA). DNA samples meeting the quality requirements were used for subsequent SNP genotyping.

SNP genotyping was performed on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) using whole-genome resequencing technology. The library construction process was as follows: genomic DNA was fragmented by ultrasonic disruption, the fragmented products were end-repaired and ligated to barcoded adapters, the ligated products were amplified by PCR, and fragments of 300~350 bp were selected to construct the library. After quality inspection, the library was subjected to paired-end sequencing (PE150). Raw sequencing data were filtered to remove adapter-contaminated sequences, low-quality sequences (Q < 20), and sequences with N content > 5%. The filtered clean reads were aligned to the reference genome of the L. crocea Daiqu strain (GenBank accession number: PRJNA1168280) using BWA software (version 0.7.17), and SNP detection was performed using GATK software (version 4.3.0.0) [59]. High-quality SNPs were filtered under the following criteria: missing rate < 20%, minor allele frequency (MAF) > 0.05, and Hardy–Weinberg equilibrium p > 0.001. Prior to genetic map construction, SNP markers were subjected to bin mapping instead of direct use. Briefly, SNPs were first imputed and corrected: a 15-SNP sliding window with a 1-SNP step size was used for chromosome scanning; windows with no less than 13 SNPs typed as “aa” or “bb” were classified as aa or bb, respectively, and other types were imputed and corrected to “ab”. Bins were then defined as continuous non-recombinant SNP clusters with identical genotypes in offspring. Bins shorter than 10 kb and those with severe segregation distortion (chi-square test, p < 0.0001) were filtered out. The genetic material origin of bin markers in offspring was visualized by graphical genotype analysis, and the resulting bins were used as polymorphic markers for genetic map construction.

4.4. High-Density Genetic Map Construction and QTL Mapping

A high-density genetic map was constructed using Lepmap3 software (version 0.2) based on the double pseudo-testcross strategy and the filtered bin markers. Maternal-specific markers (lm × ll) and paternal-specific markers (nn × np) were grouped with LOD thresholds of 3.0 and a maximum recombination rate of 0.4. Marker order was optimized using the Regression Mapping method, and genetic distances between markers were calculated by the Kosambi mapping method [60]. Shared markers (hk × hk) were used as anchors to integrate parent-specific linkage maps into a consensus map. The quality of the linkage map was assessed using metrics including the number of linkage groups, maximum marker spacing, recombination relationship evaluation, and genomic collinearity.

Based on the high-density genetic map and phenotypic data of the F1 generation, QTL mapping was performed for 9 growth-related traits and 3 disease resistance-related traits using the Interval Mapping (IM) method in MapQTL 6.0 software (Kyazma B.V., Wageningen, The Netherlands). The parameter settings were as follows: step size of 0.1 cM and scanning range covering all linkage groups. The LOD threshold for significant QTLs was determined by 1000 permutation tests (significance level α = 0.05), a widely accepted approach for controlling genome-wide type I error in QTL mapping studies [55]. Permutation tests accounted for the multiple testing inherent in genome-wide scans by generating empirical significance thresholds based on the actual data structure without relying on assumptions about the number of independent tests. QTLs with LOD values ≥ 3 were considered significant.

4.5. Data Statistics and Visualization

All data processing and statistical analyses were performed using R software (version 3.5.0) and SPSS 26.0 software (IBM Corp., Armonk, NY, USA) on phenotypic data from 150 experimental samples. Pearson or Spearman correlation coefficients were used for trait correlation analysis, with significance verified by the t-test. Result visualization was implemented using R packages (version 3.5.0) and GraphPad Prism 9.0 (GraphPad Software, Inc., Boston, MA, USA): boxplots were used to display QTL mapping results, phenotypic differences, and marker effects; scatter plots and heatmaps were used to visualize trait correlations.

5. Conclusions

This study successfully constructed a high-density SNP genetic map for the large yellow croaker and employed it to systematically dissect the genetic architectures underlying yield, morphology, and resistance to VWND. The ultra-fine map resolution enabled precise QTL detection, revealing a landscape where the genetic control of growth and disease resistance is largely distinct. The strong phenotypic and genetic integration within yield and morphological traits contrasts with the limited number of co-localized QTLs shared between growth and disease resistance. The candidate genes found within these rare co-localization intervals further support a model of weak pleiotropy, primarily involving systemic stress response pathways rather than direct regulatory hubs for growth or immunity. Our primary conclusion is that the genetic correlation between growth (both yield and morphology) and VWND resistance in L. crocea appears to be weak within this F1 family, consistent with largely independent polygenic architectures. This finding is of significant practical importance as it suggests that antagonistic selection responses are not an inevitable constraint. Breeders can therefore prioritize selection for high yield, desirable morphology, and enhanced disease resistance concurrently, accelerating the development of balanced, high-performance strains. The identified QTL hotspots and candidate genes provide valuable targets for marker-assisted selection and future functional studies.