Construction of a Growth Model and Screening of Growth-Related Genes for a Hybrid Puffer (Takifugu obscurus ♀ × Takifugu rubripes ♂)

Abstract

The obscure puffer (Takifugu obscurus) is a popular cultured species and accounts for around 50% of the total pufferfish production in China. A hybrid puffer was generated by crossing a female obscure puffer with a male tiger puffer (T. rubripes). Its growth model has not been developed and the genetic basis underlying its growth superiority has not yet been fully investigated. In this study, the growth model and morphological traits of the hybrid puffer were explored. The results indicated that the hybrid puffer exhibited a significant growth advantage compared to the obscure puffer. There were also significant differences in their morphological traits. We conducted genotyping-by-sequencing (GBS) on hybrid and obscure puffer groups, identifying 215,288 high-quality single nucleotide polymorphisms (SNPs) on 22 chromosomes. Subsequently, a total of 13 growth-related selection regions were identified via a combination of selection signatures and a genome-wide association study (GWAS); these regions were mainly located on chromosomes 10 and 22. Ultimately, the screened regions contained 13 growth-related genes, including itgav, ighv3-43, ighm, atp6v1b2, pld1, xmrk, inhba, dsp, dsg2, and dsc2, which regulate growth through a variety of pathways. Taken together, the growth models and candidate genes used in this study will aid our understanding of production characteristics and the genetic basis of growth rates. The hybrid will also be of great significance for the genome-assisted breeding of pufferfish in the future.

1. Introduction

Pufferfishes, especially members of the genus Takifugu, are economically important aquaculture species in East Asia owing to their desirable taste and high nutritional quality. In China, due to their tetrodotoxin (TTX) content, only two Takifugu species have been approved for artificial cultivation and consumed, namely, tiger puffer (Takifugu rubripes) and obscure puffer (T. obscurus), since 2016. Their total output increased to 31,060 tons in 2022 [1]. Normally, the tiger puffer exhibits a very fast growth rate in a farming environment and has the largest body size within the Takifugu genus [2]; however, this species resides obligately in seawater and is only allowed to be farmed in coastal areas. The obscure puffer is a migratory species that has acclimated to completely freshwater regions [3,4], and it is gradually becoming one of the main Takifugu species used for production. Nevertheless, compared with the tiger puffer, the obscure puffer has a smaller body size and a slower growth rate.

Growth is one of the most significant traits in the aquaculture industry, and it is directly related to yield and economic efficiency [5]. At the same time, growth rate is a crucial factor that optimizes feeding rates and breeding density. Therefore, constructing a growth model to estimate or predict fish weight is essential for achieving better control of the output of an aquaculture system [6]. Hybridization, a significant breeding method, has been widely used in aquatic breeding since the 1950s [7]. It can be effective to combine the advantageous traits of both parental species to rapidly improve the traits of the offspring, especially growth [8,9,10]. For example, in terms of growth and stress resistance, a high growth rate and high spot disease resistance were induced by the distant hybridization of Percocypris pingi (♂) and Schizothorax wangchiachii (♀) [11]; further, a high growth rate and hypoxia-tolerance traits of hybrid Pelteobagrus fulvidraco (♀) × Leiocassis longirostris (♂) were observed [12]. Nevertheless, the enhancement of growth increases the difficulty of aquaculture prediction. For example, hybrid Jinhu grouper (Epinephelus fuscoguttatus ♀ × E. tukula ♂) and Hulong grouper (E. fuscoguttatus ♀ × E. lanceolatus ♂) groups showed significantly higher growth than the parental E. fuscoguttatus group in one study, and there were different growth models between the groups, all of which showed allometric relationships between body weight and body length (W = aL^b, b < 3) [13]. Overall, the hybrid group exhibited a significantly improved growth rate, accompanied by alterations in growth models and morphological traits. These observations offer insights into the growth advantages of hybrid groups.

Extensive research has been conducted on hybrid puffers to study factors such as the slow growth rate and diverse osmoregulatory abilities of the pufferfishes. The exploration of hybrid puffers began in the 1960s [14]. Early studies focused on the embryonic development [14], morphological characteristics [15], growth performance, and osmoregulation of the hybrid puffer, with results showing that the hybrid puffer (T. flavidus ♀ × T. obscurus ♂) cannot be hatched in freshwater [16] and that the hybrid puffer (T. flavidus ♀ × T. rubripes ♂) exhibits different morphological characteristics from its parents [15]. In recent years, with the development of molecular biology, research on the hybrid puffer has involved various aspects, such as mitochondrial DNA sequences [17], gut microbiota [18], body composition [19], and screening for salinity regulation genes [20]. However, until now, for the hybrid puffer (T. obscurus ♀ × T. rubripes ♂), there is still a lack of a long-term growth models, growth-related SNPs, and growth-related gene studies.

Recently, numerous studies have demonstrated that variations in growth and morphological traits are associated with genetic diversity among groups [7]. Additionally, genomic analysis methods, such as selection signature analysis [10] and genome-wide association studies (GWASs) [21], have demonstrated that the molecular mechanisms of growth traits can be explored and that the growth-related selection regions, genes, and single nucleotide polymorphisms (SNPs) can be further investigated. For example, regarding the Australasian Snapper (Chrysophrys auratus), selection signature analysis was used to identify growth genes associated with fast-growing groups [22]. Similarly, the growth-related genes of mandarin fish (Siniperca chuatsi), including rnf213, mkk6, and nck2, were screened using a genome-wide association study (GWAS) [23]. Nevertheless, there is still a paucity of analyses of growth performance and growth models for hybrid pufferfish, and their growth-related molecular regulatory mechanisms are still unknown.

In this study, growth models and morphological traits were analyzed between the hybrid and obscure pufferfishes. Subsequently, individuals with different growth traits were genotyped. Concurrently, selection signatures and a GWAS were employed to investigate the growth-related SNPs and selected regions to identify the genes underlying the genetic adaptation responsible for the rapid growth observed in hybrid puffers. Taken together, this study will improve comprehension of the genetic mechanism of pufferfish growth and provide important theoretical support for genome-assisted breeding.

2. Materials and Methods

2.1. Construction of Hybrid Group

A hybridization experiment was conducted on 22 February 2022, at Zhongyang Aquatic Co., Ltd. (Nantong, China), to generate hybrid puffers between randomly selected parental female T. obscurus and male T. rubripes groups (Figure 1). The parental fishes were approximately 3–4 years old, with body weights ranging from 700 to 1500 g. Following a month of conditioning, sexually mature individuals, including randomly selected groups of 10 male obscure puffers, 20 female obscure puffers, and 10 male T. rubripes individuals, were used for breeding under identical conditions (400 m², 18.0 ± 1 °C). Their body weights were 777.63 ± 186.37 g, 859.19 ± 183.33 g, and 1195.50 ± 149.22 g, respectively (Table S1). The selected fish were injected with human chorionic gonadotropin (HCG) and luteinizing hormone-releasing hormone (LHRH-A2) [24,25]. The two-injection method [26] was employed in this study; the first injected dose administered to the females was 200–300 IU/kg of HCG + 50–100 μg/kg of LRHA-A2, while the males were not injected. The second injected dose administered to females was 200–300 IU/kg of HCG + 100–300 μg/kg of LRHA-A2, while the males were injected with half this dose [11]. Following an interval of 24–48 h, the abdomens of the mature parental fish were gently pressed to obtain approximately 500 g of mature eggs and 20 mL of sperm; then, hybrid and obscure puffer groups were created (Figure 1).

Individuals of obscure and hybrid puffers were selected at 0 days post-hatching (dph) after the hatching of fertilized eggs. These puffers were reared at different culture stages (0–270 days post-hatching, dph) for growth comparison. During the initial 30-day period (0–30 dph), two 20 m²-size tanks were used and each pond contained 5000 hybrid puffers and 5000 obscure puffers, with sterilized groundwater (25 ± 1 °C, 0‰) that was changed every 10 days. The fish in both groups were fed rotifers (Brachionus calyciflorus) and brine shrimp (Artemia salina) five times daily at specific intervals (6:00, 9:00, 12:00, 15:00, and 17:00) [27], maintaining a rotifer density of 5000–6000/L and a brine shrimp density of 20–50/L in the water column with each feeding. After 30 days, considering the loss of fish and the breeding density, 3000 healthy individuals each of hybrid and obscure puffers were randomly selected and transferred to same 444 m² tanks, with the same water and temperature. The fish were fed commercial feed (Zhongyang brand, pufferfish variety, containing 40% crude protein) three times daily at specific intervals (8:00, 14:00, and 18:00), the feeding situation was observed for 60 min after each feeding, and the remaining food was removed in time if no fish were eating. Freshwater (25 ± 1 °C and 7–8 mg/L dissolved oxygen) was supplied, and replaced every 10 days.

2.2. Growth and Morphological Traits Comparison

During the experiment, electronic scales and vernier calipers, with precisions of 0.01 g and 0.01 mm, respectively, were utilized to measure the growth traits and morphological traits.

The body weight (BW), body length (BL), and total length (TL) of 100 randomly selected individuals from each group were measured at 30, 60, 90, 120, 150, 180, 210, and 270 dph. The weight gain ratio (WGR), absolute growth rate (AGR), specific growth rate (SGR), and condition factor (CF) were calculated according to the following equations [28,29,30]:

WGR (%) = (W_i − W₀)/W₀ × 100%

(1)

AGR (g d⁻¹) = (W_i − W₀)/(t_i − t₀)

(2)

SGR (% d⁻¹) = (ln_Wi − ln_W0)/(t_i − t₀) × 100%

(3)

CF (%) = (W_i/L_i³) × 100%

(4)

where t is the experiment duration (d), W_i is the weight of t_i, W₀ is the weight of t₀, and L_i is the length of W_i.

The morphological traits of 100 randomly selected samples from each group were measured at 270 dph. The sex of each fish was determined based on a sex-linked SNP in the amhr2 gene [31]. Fifteen morphological traits were measured [30]: BW, TL, BL, caudal peduncle height (CH), head length (HL), snout length (SL), head-behind length (HBL), eye length (EL), eye spacing (ES), nostril spacing (NS), outlet hole spacing (OS), snout cleft (SC), chest length (CL), abdominal length (AL), and caudal girth (CG) (Figure 2). To reduce the effects of excessive individual weight differences, the data obtained were corrected using BL in analyses of morphological traits. The Spearman correlation coefficient was utilized to conduct a correlation analysis of morphological traits among two populations using SPSS (version 26.0). The correlation threshold was established as follows: |r| < 0.4 indicated a low correlation; 0.4 ≤ |r| ≤ 0.7 indicated a moderate correlation; and |r| > 0.7 indicated a high correlation. Based on the results of a normal distribution analysis, a principal component analysis (PCA) of the morphological traits was performed using SPSS (version 26.0). Simultaneously, 75 samples were randomly selected from each group. Additionally, a cluster analysis of two groups was performed using OriginPro (version 2019b) [32] based on the unweighted pair-group method with arithmetic means (UPGMA) and morphological traits.

The functional relationship between BW and BL was calculated according to the following equation [27]:

y = a x^b

(5)

where x is body length; y is body weight, a is the intercept on the y-axis, and b is the allometric growth index.

The von Bertalanffy growth models were calculated according to the following equations [33,34]:

W_t = W_∞ (1 − e^−k(t−t₀⁾)

(6)

L_t = L_∞ (1 − e^−k(t−t₀⁾)

(7)

where W_t and L_t are body weight and body length at age t, respectively; L_∞ is the asymptotic maximum body length; W_∞ is the asymptotic maximum body weight; k is the Brody coefficient (or “growth constant”), describing how fast W_∞ and L_∞ are approached; and t₀ is the theoretical age at size 0 in the von Bertalanffy growth model.

The inflection point age of growth was calculated according to the following equation:

T_t = ln_b/k + t₀

(8)

where T_t is the time taken to reach the growth inflection point; b is the allometric growth index; k is the Brody coefficient (or “growth constant”), describing how fast W_∞ and L_∞ are approached; and t₀ is the theoretical age at size 0 in the von Bertalanffy growth model.

Finally, the Fisher and Bayes discriminant equations for the two groups were calculated using stepwise discriminant analysis based on Wilks’ lambda method [32] and different morphological traits.

2.3. Sample Collection and Sequencing

At the age of 270 days, 100 individuals were taken from both the obscure and hybrid puffer groups. These individuals were then categorized as either fast-growing or slow-growing based on their body weights. Then, 20 fast-growing individuals and 20 slow-growing individuals were selected from each group. A total of 80 pufferfish samples were collected (Table S2). Their growth traits (BW, BL, TL, AL, and CL) were measured, and the fins of each fish were clipped and preserved in 95% ethanol at −20 °C. DNA was extracted using the cetyltrimethylammonium bromide (CTAB) method, while DNA quality was measured using a Qubit fluorometer (Thermo Fisher Scientific, Carlsbad, CA, USA) and a NanoDrop spectrophotometer (Thermo Fisher Scientific).

Based on the size and GC content of the T. rubripes genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_901000725.2/, accessed on 15 July 2023), appropriate endonucleases for genomic fragmentation were predicted. The qualified genomic DNA was digested using restriction enzymes, followed by end repair, the A-tail was added, and an Illumina sequencing connector was added by the NEBNext Ultra DNA Library Prep Kit (NEB, Ipswich, MA, USA). DNA fragments (300–400 bp) were amplified and enriched using PCR. Finally, the PCR products were purified using the AMPure XP system (Beckman Coulter, Brea, CA, USA). Sequencing libraries were examined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and quantified using qPCR. Sequencing was performed by employing a NovaSeq 6000 sequencer using the PE 150 sequencing strategy.

2.4. Quality Control and Genotyping

The reads were subjected to quality control using FASTP software (version 0.18.0) [23,35], Low-quality reads were removed, and low quality was defined as follow: reads aligned to the barcode adaptor; reads with connectors; reads including ≥ 10% unidentified nucleotides (N); and reads for which ≥ 50% of the bases had a Phred quality score ≤ 20. The filtered reads were mapped against the reference genome using the MEM algorithm in BWA (version 0.7.12) [35,36], and the comparison parameter was–k 32 –M. After comparison, the potential PCR duplicates were tagged and removed using PICARD (version 1.129) (Picard: http://sourceforge.net/projects/picard/, accessed on 10 August 2023). SNPs were identified and filtered using GATK (version 3.4-46) [37] Variant Filtration with proper standards (–Window 4, –filter “QD < 4.0 || FS > 60.0 || MQ < 40.0”, –G_filter “GQ < 20”), and functional annotation of the detected variants was performed using ANNOVAR (version 2) [38]. In addition, to enhance the accuracy of this study, we applied stringent filtering criteria to the SNPs using PLINK2 (version 2.0) [39]. Stringent filtering conditions were set as follows: for selection signatures, the missing rate was < 50%, and for GWAS, the filtered SNPs were all indels and non-biallelic SNPs, SNPs with a minor allele frequency (MAF) < 0.05, SNPs with a call rate < 0.9, individuals with a missing data rate > 0.5, and heterozygosity rate > 0.8. Finally, these analyses of SNPs were used in the subsequent GWAS and selection signatures analyses.

2.5. Analysis of Population Structure and Kinship

After SNP quality control, linkage disequilibrium (LD) pruning was performed according to the r² value. LD pruning was conducted with a window size of 50, a step of 5, and an r² threshold of 0.2 [40,41]. According to the filtered SNPs, principal component analysis (PCA) and kinship analysis were performed on the experimental population using GCTA software (v1.93.2) [42]. Simultaneously, an NJ tree (model: p-distance; bootstrap replications 1000) was constructed using MEGA-X software (version 10) [43], and the population structure was constructed using Admixture software (version 1.3) [44]. Since the experimental population comprises individuals with different growth rates and hybrid individuals, it was hypothesized that the samples belong to multiple subgroups. Therefore, the number of subgroups (K) for the samples was assumed to range from 1 to 9. Subsequently, the optimal number of clusters was determined based on the cross-validation error rate (CV error).

2.6. Screening for Selection Signatures and a Genome-Wide Association Study

The F_st (a measure of population differentiation, varies from 0 to 1) and π ratio were used to identify the candidate selected region. PopGenome software (in R package) [45,46] was used to perform sliding-window analysis based on the filtered SNPs, using a window size of 100 kb and a step size of 10 kb according to physical length. This step facilitated a comparative evaluation of intra- and inter-population diversity, encompassing analyses of F_st [47] and the π ratio [48]. The selected regions were determined by comparing the top 1% of results from the F_st and π ratio, ranked by the degree of genetic differentiation.

A GWAS was used to investigate the genetic architecture of the growth-related traits. According to a modified Bonferroni correction method [23,49], the significance and suggestive association thresholds were set at 0.05/N and 1/N, respectively, where N is the number of SNPs used in the association analysis. SNP density maps, Manhattan plots, and quantile-quantile (Q-Q) plots were generated using CMplot in R [23,50].

2.7. Identification and Functional Annotation of Candidate Genes

SNPs were mapped against the T. rubripes genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_901000725.2/, accessed on 22 August 2023). In the selection signatures, the growth-related candidate genes were identified by scanning candidate selection regions; in the GWAS analysis, the genes containing significant SNPs were scanned, and the genes of ± 300 kb of significant and suggestive SNPs were scanned. Subsequently, to enhance the accuracy of those analysis, we performed an intersection analysis of those genes.

All candidate genes were functionally annotated based on the reference genome of T. rubripes posted by the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using OmicShare (www.omicshare.com/tools) and Metascape (https://metascape.org/gp/index.html#/reportfinal/tm8hy07mo, accessed on 30 August 2023), respectively [40]. The calculated p values were subjected to FDR correction, with a threshold of FDR ≤ 0.05. GO terms that met this criterion were considered significantly enriched. Similarly, the calculated p values underwent FDR correction using the same threshold, and pathways meeting this condition were defined as significantly enriched.

2.8. Statistical Analysis

The data were analyzed using SPSS (version 26.0) [32]. Morphological traits were corrected for body weight. The data were checked for the normality of distribution using the Kolmogorov–Smirnov test. Bartlett’s test was used to determine the homogeneity of variance among different groups. Significant differences were analyzed via a t-test, in which p < 0.01 and p < 0.05 were considered highly significant and significant differences, respectively. Data are expressed as means ± standard deviation.

3. Results

3.1. The Disparities in Growth Performance

At 270 dph, the BW, BL, and TL of the hybrid were significantly higher than those of the obscure puffer (Table S3, Figure S1, p < 0.01); the BW and TL of the hybrid reached 266.13 ± 63.90 g and 22.00 ± 1.6 cm, respectively. At 120 dph, due to the health condition of the experimenters, the hybrid puffer did not meet the specified measurement target and therefore was excluded from this analysis. The subsequent analysis remains unaffected. The hybrid achieved a maximum AGR of 1.09 g/d, which is 1.7 times higher than that of the obscure puffer, and the maximum SGR reached 2.11%/d, indicating that the hybrid exhibited a significant growth advantage. The AGR and SGR results indicate that the growth rate decreased with increasing age in both groups, while the hybrid still exhibited a significantly higher growth rate than the obscure puffer (

Table 1. Growth traits of hybrid and obscure puffer groups.

Stage		Obscure Puffer				Hybrid Puffer
Family	DPH	WGR	AGR	SGR	CF	WGR	AGR	SGR	CF
Larval	60	10.48	0.17	4.07%	4.97	17.93	0.3	4.90%	3.56
	90	50.56	0.56	4.38%	4.24	67.22	0.88	4.69%	4.2
Juvenile	120	75.1	0.63	3.61%	3.87	-	-	-	-
	150	92.77	0.62	3.03%	3.85	122.31	0.82	3.08%	4.09
	180	123.99	0.69	2.68%	3.63	184.23	1.02	2.90%	4.32
	210	137.85	0.66	2.35%	3.54	203.16	0.97	2.53%	3.74
Young	270	173.31	0.64	1.91%	3.62	294.29	1.09	2.11%	3.97

DPH, day post-fertilization; WGR, weight gain ratio; AGR, absolute growth rate; SGR, specific growth rate; CF, condition factor.

, p < 0.01). In addition, the CF of the hybrid was lower than that of the obscure puffer prior to 60 dph (p > 0.05), a finding we hypothesize is related to the growth model differences between the hybrid and obscure puffers. In conclusion, the hybrid exhibited a significant growth advantage over the obscure puffer.

3.2. The Variations in Morphological Characteristics

To reduce the effects of excessive individual weight differences, the data were corrected using BL in the analyses of morphological traits. These standardized traits were comparatively analyzed (Table S4), and the results showed that 13 traits were significantly different (BW/BL, TL/BL, CH/BL, SL/BL, HBL/BL, EL/BL, ES/BL, OS/BL, SC/BL, CL/BL, AL/BL, CL/BL, and CF/BL; p < 0.01), while two traits were not significantly different (HL/BL and NS/BL; p > 0.05), between the hybrid and obscure puffers. Spearman correlation analysis of the morphological traits and sex (Figure S2, Table S5) showed that in the obscure puffer, the CH/BL (|r| = 0.50), CL/BL (|r| = 0.46), AL/BL (|r| = 0.53), CG/BL (|r| = 0.47), and CF/BL (|r| = 0.61) were significantly positively correlated with BW/BL (p < 0.01); the EL/BL (|r| = 0.30) was significantly negatively correlated with BW/BL (p < 0.01); the BW/BL (|r| = 0.3) was significantly positively correlated with sex (p < 0.01); the HL/BL (|r| = 0.34) was significantly negatively correlated with sex (p < 0.01); and the SL/BL (|r| = 0.25) and HBL/BL (|r| = 0.28) were significantly negatively correlated with sex (p < 0.05). In the hybrid puffer, CL/BL (|r| = 0.43), AL/BL (|r| = 0.51), CG/BL (|r| = 0.49) and CF/BL (|r| = 0.71) showed a significantly positively correlated with BW/BL (p < 0.01); the CH/BL (|r| = 0.26) was significantly positively correlated with BW/BL (p < 0.05); and the EL/BL (|r| = 0.45) was significantly negatively correlated with BW/BL (p < 0.01); only ES/BL (|r| = 0.26) showed a significant positive correlation with sex (p < 0.05); there were no growth-related traits that had a significant correlation with sex (p > 0.05). In conclusion, the above results indicate that there is a significant difference in morphological traits between the hybrid and obscure puffers, and that the growth of the obscure puffer was significantly associated with sex. These morphological trait differences may allow distinction between these two groups, which are difficult to distinguish by appearance, through morphological features. Additionally, the differences in morphological traits also interact with the growth patterns, indicating the existence of different growth models between the obscure and hybrid puffers.

Based on the morphological traits, a principal component analysis (PCA) was conducted on a mixed group of 150 individuals with a normal distribution (75 hybrid and 75 obscure puffers). The results showed that the morphological traits were grouped into four categories (Figure 3). The cumulative contribution rate was 80.97%. The contribution rate of the first principal component was 39.53%, including AL/BL, CL/BL, CG/BL, BW/BL, CF/BL, CH/BL, and EL/BL, which reflected a pufferfish body shape. The second principal component was 15.48%, including HL/BL and HBL/BL, which reflected the anterior torso ratio. The third principal component accounted for 15.34%, including OS/BL, ES/BL, and SC/BL, which reflected the morphological structure of the head. According to the analysis of PCA1, PCA2, and PCA3, the obscure and hybrid puffers from PCA1 can be better distinguished in terms of body shape. From the aspect of anterior torso ratio in PCA2, the two groups are difficult to distinguish due to a large amount of overlap. Although there are certain differences in the morphological structure of the head from the perspective of PCA3, there were still several individuals that were difficult to clearly classify. In addition, even from the perspective of body shape in PCA1, the two groups are distributed on both sides and can be distinguished well, but there are still a few individuals that overlap. Similarly, there are overlapping individuals in the cluster tree results (Figure S3). We analyzed these overlapping individuals and the results show that the overlapping individuals in the PCA and cluster tree analysis are largely the same individuals. Then, analysis of the morphological traits of overlapping individuals showed that these individuals are not only affected by body weight but also other morphological traits, such as body length and chest length. Among overlapping individuals, the hybrid puffer usually exhibits weight gain and rapid elongation compared to other individuals in the same group. This result further suggests that there are different growth models between the obscure and hybrid puffer groups. In addition, the different growth models were genetically controlled; therefore, it can be speculated that the morphological traits of hybrid puffers among overlapping individuals may be influenced by genetics. The relationship between genetic and morphological traits is analyzed in the following section, “3.5 Analysis of population structure”. In general, there are multiple morphological differences between obscure and hybrid puffers, and the two groups can be well distinguished using morphological traits related to body shape.

3.3. The Growth Model and Discriminant Analysis

The relationship between BL and BW was calculated. The results showed that b = 2.759 and R² = 0.99 for the obscure puffers and b = 3.072 and R² = 0.99 for the hybrid puffers (Figure 4), indicating that the obscure puffer (b < 3) showed negative allometric growth, while the hybrid (b ≈ 3) showed uniform growth. Growth is determined by the regulation of genes. Therefore, we inferred that the differences in growth models are related to genetics, and these different growth models lead to increased morphological trait differences in body shape between the two groups during the growth and development process.

The von Bertalanffy growth models were also calculated (Figure 5). The result showed that the growth inflection point age of obscure puffer was 96 dph, and the corresponding BW and BL were 48.04 g and 10.69 cm, respectively. Meanwhile, the theoretical maximum BW and BL at 270 dph were 165.64 g and 16.64 cm (R² = 0.99), respectively. Similarly, the growth inflection point age of hybrid was 167 dph, and the corresponding BW and BL were 145.93 g and 15.37 cm, respectively; the theoretical maximum BW and BL at 270 dph were 489.23 g and 22.79 cm (R² = 0.99), respectively, in the hybrid.

Based on morphological disparities, stepwise discriminant analysis was conducted using Fisher and Bayes discriminant equations. Stepwise discriminant analysis demonstrated that the BW/BL, EL/BL, SC/BL, AL/BL, CG/BL, and CF/BL ratios exhibited the highest discriminatory power. According to Fisher’s principle, the discriminant equation created was as follows:

Y = 0.886 ∗ BW/BL + 437.567 ∗ EL/BL − 39.976 ∗ SC/BL + 15.646 ∗ AL/BL + 35.454 ∗ CG/BL − 5.744 ∗ CF/BL − 36.175

(9)

where Y is the type of group; Y > 0 indicates the hybrid puffer. Y < 0 indicates the obscure puffer.

According to Bayes’ theorem, the discriminant equations were derived as follows:

Obscure puffer = 25.088 ∗ BW/BL + 19385.331 ∗ EL/BL + 2570.616 ∗ SC/BL + 1001.791 ∗ AL/BL + 1264.575 ∗ CG/BL − 244.468 ∗ CF/BL − 933.933

(10)

Hybrid puffer = 32.477 ∗ BW/BL + 23035.524 ∗ EL/BL + 2237.132 ∗ SC/BL + 1132.311 ∗ AL/BL + 1560.337 ∗ CG/BL − 292.387 ∗ CF/BL − 1233.482

(11)

Individual morphological traits were measured and calculated; if obscure puffer > hybrid, the individual was an obscure puffer and vice versa for a hybrid.

Furthermore, an additional 50 individuals who were not included in the analysis (25 obscure puffers and 25 hybrid puffers) were randomly selected to further validate the equations. The identification rate for both was 96%.

3.4. Quality Control of the Sequencing Reads

A total of 556.50 M reads (79.02 Gbp) were obtained, with 297.02 M reads (42.18 Gbp) and 259.47 M reads (36.85 Gbp) for the hybrid and obscure puffer groups, respectively. The sequencing depths were 96 X and 110 X, respectively. The proportions of high-quality bases (Q score > 20) were 95.90% and 95.73%, respectively. The average GC contents were 45.57% and 45.50%, respectively. After processing and screening, 215,788 and 182,925 high-quality SNPs distributed across 22 chromosomes were included in the selection signatures and GWAS, respectively (Figure 6). The SNPs involved in the GWAS comprised 114,665 transitions and 68,260 transversions. In the transitions, C -> T occupied a maximum of 34,235. In the transversions, C -> A occupied a maximum of 9639. However, the distribution of these SNPs in the genome was not uniform. Following annotation, 23.48% (48,714) of the variants were identified in intergenic regions, 52.30% (108,498) were identified in intronic regions, and only 6.15% (12,759) were identified in exonic regions (Table S6). Among them, 57.02% of the SNPs in the exon region were synonymous mutations, while 40.60% were nonsense mutations and missense mutations (Table S7).

3.5. Analysis of Population Structure

Based on SNPs, the results regarding kinship showed that the genetic relationship coefficient ranged from −0.2 to 0.4, and the two groups contained subgroups with close genetic relationships (Figure S2). The results of the genetic PCA showed that PCA1 and PCA2 accounted for 27.55% of the total phenotypic variation (Figure 6). Specifically, PCA1 explained 19.75% of the total genetic variation and PCA2 explained 7.80%. Meanwhile, the results showed that there is a notable distinction between the obscure and hybrid puffer groups. This is the same as the PCA results for both of their morphological traits. Moreover, the genetic PCA results showed that both groups exhibited subgroup structures (Figure 6). Similarly, the NJ tree result also showed that there are sister groups within the two groups (Figure 6). We combined morphological traits and kinship to analyze the subgroup structures, and the results showed that the subgroup structure was related to morphological traits and kinship. Among the subgroups, one subgroup showed similar growth performance amongst its members, while the growth performance between different subgroups was significantly different. In addition, the analysis of population structure showed that the optimal number of subgroups is K = 4, and there were two subgroups within each of the two groups (Figure S4). In the different clustering results, the obscure puffer can be clearly distinguished, whereas there are few hybrid puffers with the same bloodline as the obscure puffer, indicating that the hybrid puffer and the obscure puffer share a common ancestor. Analysis of the hybrid puffer with an overlapping background shows that most of the individuals with an overlapping background grow slowly and have morphological traits similar to the obscure puffer. Combining the PCA results regarding morphological traits, the results showed that these hybrid puffers were the individuals that were overlapping with the obscure puffer. Therefore, this finding indicates that genetics significantly influence the growth performance of hybrid puffers, and the genetic traits of fast-growing individuals differ significantly from those of slow-growing individuals. Researching these genetics provides support for the analysis of growth mechanisms. Interestingly, when K = 5 or K = 6, the cross-validation error values are close to K = 4. We interpret that this is because the experimental group is built from a smaller parental group. The small parental groups increase the probability of producing multiple lineage groups during the fertilization process. However, excessive subgroup numbers result in an uneven distribution of individuals across different subgroups, which is not conducive to GWAS analysis. Consequently, this study was not adopted. Overall, the results imply that morphological traits are interrelated with genetics, and there are significant genetic differences between the obscure and hybrid puffers leading to differences in their growth and morphological traits. Both experimental groups have certain kinship and subgroup structures; therefore, to improve the accuracy of the GWAS, a mixed linear model (MLM) was employed using GEMMA (version 0.98.1) [51].

3.6. Detection of Genome Selection Signatures

Based on the constructed population pairs exhibiting extreme differences in growth-related phenotypic traits, we used the F_st and π ratio methods to identify the genomic signatures of selection in the hybrid and obscure puffer groups. We obtained a set of 1720 genes from a total of 860 candidate selected regions in the top 1% of selection regions (F_st ≥ 0.60, π ratio ≥ 1.05, Figure 7) by combining the F_st and π ratio.

The F_st method identified a total of 445 candidate selected regions encompassing 823 genes (Table S8), which were mainly enriched on chromosomes 4, 10, and 22, while the π ratio method identified 415 candidate selected regions encompassing 944 genes (Table S8), which were mainly enriched on chromosomes 10 and 22. In addition, a total of 30 intersection candidate regions, encompassing 47 genes, were identified (Table S8).

3.7. GO Term and KEGG Pathway Analysis

GO and KEGG analyses were conducted on the candidate genes within the selected regions to provide insights for subsequent investigations. The results obtained from the GO analysis of the intersection of F_st and π ratio analyses (Figure S5) indicated that 57.45% of the genes (n = 27) were categorized into the biological process (BP) category, with 633 enriched terms, of which 118 were significant; 55.32% of the genes (n = 26) were categorized as cellular components (CC), with 93 enriched terms, of which 11 were significant; and 59.57% of the genes (n = 28) were categorized as molecular function (MF), with 95 enriched terms, of which 14 were significant (Table S9). The top 20 enriched GO terms are shown in Figure 8, ranked by p value. In biological processes, the most abundant terms included regulation of actin filament-based movement (GO:1903115), regulation of muscle contraction (GO:0006937), regulation of blood circulation (GO:1903522), and regulation of muscle system processes (GO:0090257), all of which involved an average of eight genes (p < 0.01). In the cellular components, the cell–cell junction (GO:0005911), cell–cell adherens junction (GO:0005913), and proton-transporting V-type ATPase complex (GO:0033176) all contained an average of four genes (p < 0.05). In terms of molecular function, protein binding was involved in cell adhesion (GO:0098631), insulin-like growth factor binding (GO:0005520), and inorganic cation transmembrane transporter activity (GO:0022890), all of which involved an average of two genes (p < 0.05).

A KEGG pathway enrichment analysis of the candidate genes, ranked by p value, was conducted (Figure 9). The top 20 enriched KEGG terms in the F_st analysis are presented; the growth-related pathways specifically enriched included hematopoietic cell lineage (ko04640), mannose type O-glycan biosynthesis (ko00515), and vitamin digestion and absorption (ko04977) (Table S10). The π ratio analysis (Figure S6) included cell adhesion molecules (ko04514) and cellular senescence (ko04218) (Table S10). In addition, the result of the analysis of all candidate genes indicated that the human diseases category exhibited the highest abundance, and the other pathways are mainly enriched in organismal systems and metabolism (Figure S6).

3.8. Genome-Wide Association Analysis of Growth Traits

The original phenotypic data for each sample were measured and recorded. The growth-related traits exhibited predominantly non-normally distributed patterns according to the Kolmogorov–Smirnov test; therefore, these traits were transformed to achieve normality. Subsequently, a GWAS was conducted on five growth-related traits (BW, TL, BL, CL, and AL) using GEMMA software (version 0.98.1), encompassing 182,925 SNPs.

For the obscure puffer, the GWAS showed 32 SNPs significantly associated (−log ₁₀ > 4.7) (Figure S7, Table S11) and 163 suggestively associated (−log ₁₀ > 3.4) with growth traits. These SNPs were located mainly on chromosomes 1, 21, and 22. Among them, 12 SNPs were associated with multiple traits (

Table 2. The SNPs associated with multiple traits in GWAS analysis in pufferfishes.

Family	Trait	SNP	Chr	Position	Beta	p Value	Region	Gene Name
Obscure puffer	BW and AL	LOC1:2211857	1	2,211,857	−37.33243	1.83 × 10−4	Intergenic	NA
		LOC1:3193101	1	3,193,101	29.33228	2.79 × 10−4	Exonic	rnf213
	BW and TL	LOC22:5379885	22	5,379,885	32.05076	1.04 × 10−4	Intergenic	NA
	TL and AL	LOC1:3253768	1	3,253,768	−1.273611	2.08 × 10−4	Intronic	baiap2
		LOC1:3253792	1	3,253,792	−1.273611	2.08 × 10−4	Intronic	baiap2
		LOC1:3327200	1	3,327,200	−1.273611	2.08 × 10−4	Intronic	wapl
		LOC1:3403622	1	3,403,622	−1.273611	2.08 × 10−4	Intronic	grid1
		LOC1:3479403	1	3,479,403	1.435881	1.04 × 10−4	Intronic	ccser2
		LOC1:3479441	1	3,479,441	1.435881	1.04 × 10−4	Intronic	ccser2
	TL and BL	LOC3:3685861	3	3,685,861	1.141095	3.91 × 10−4	Intronic	cfap74
		LOC3:3685871	3	3,685,871	−1.141095	3.91 × 10−4	Intronic	cfap75
		LOC12:5883873	12	5,883,873	−2.1259	2.70 × 10−4	ncRNA intronic	LOC115251859
Hybrid puffer	BW and CL	LOC1:5146711	1	5,146,711	78.12444	4.94 × 10−5	Intronic	vclb
		LOC20:7178403	20	7,178,403	−58.05161	3.70 × 10−5	Intronic	lingo3
	BW and TL	LOC21:8126695	21	8,126,695	118.4649	5.24 × 10−5	Intergenic	NA

BW, body weight; BL, body length; TL, total length; CL, chest length; AL, abdominal length.

): LOC1:2211857 and LOC1:3193101 were associated with BW and AL, with SNP effects (beta) ranging from −37.33 to 29.33; LOC22:5379885 was associated with BW and TL, with SNP effects of 32.05 and 1.20, respectively; LOC1:3253768, LOC1:3253792, LOC1:3327200, LOC1:3403622, LOC1:3479403, and LOC1:3479441 were associated with TL and AL, with SNP effects ranging from −1.37 to 1.47; and LOC3:3685861, LOC3:3685871, and LOC12:5883873 were associated with TL and BL, with SNP effects ranging from −2.80 to 1.38. Although no significant SNPs related to BW and CL were found, 21 and 55 suggestive SNPs related to them also had certain reference values. In addition, LOC1:3193101 (rnf213) and LOC12:5883873 (LOC115251859) were both significant SNPs (−log₁₀ > 4.7) and associated with multiple traits. A QQ plot of the p values is shown in Figure S8.

Meanwhile, for the hybrid group, the GWAS showed that 1 SNP significantly associated (−log ₁₀ > 5.3) (Figure S7, Table S11) and 40 SNPs suggestively associated (−log ₁₀ > 4.0) with the growth traits. These SNPs were located mainly on chromosomes 1, 10, and 21. Among them, three SNPs were associated with multiple traits (Table 2): LOC1:5146711 and LOC20:7178403 were associated with BW and CL, with SNP effects ranging from −58.05 to 78.12, and LOC21:8126695 was associated with BW and TL, with SNP effects of 118.46 and 2.77, respectively. Although only 1 SNP significant is related to CL was found, 40 suggestive SNPs related to BW, TL, BL, and AL also had certain reference value; a QQ plot of the p values is shown in Figure S8.

To further explore the relationship between morphological traits and genetics, we selected SNPs associated with body weight (chr 21, T -> C) and body length (chr 10, G -> T) in the hybrid puffer for further analysis. The results showed that these SNPs were homozygous in the obscure puffer and were not associated with growth traits. However, in the hybrid puffer, most of them were heterozygous and exhibited fast growth, while only a few individuals were homozygous and exhibited slow growth. Combined with morphological trait analysis, the results showed that only homozygous individuals had morphological traits similar to the obscure puffer, and the morphological trait analysis results of some individuals were overlapping with the obscure puffer. This result is consistent with the analysis of population structure, suggesting that the morphological traits of hybrid puffer are influenced by their parents, and that genetic variants play a regulatory role in morphological trait development.

In conclusion, the morphological traits were significantly correlated with genetics, and the morphological traits of the hybrid puffer were influenced by the obscure puffer. In addition, although the number of GWAS groups involved was limited and the false-positive rate was elevated, we comprehensively scanned all genes within a 300 kb range upstream and downstream of the significant SNPs (Table S11). Subsequently, candidate genes associated with growth and metabolism were selected, thereby providing a valuable gene library for subsequent association analysis. Simultaneously, it is important to direct our attention to chromosomes harboring multiple significant SNPs, such as chromosomes 1, 10, and 22.

3.9. Combining Selection Signatures and Association Analysis

We integrated selection signatures and GWAS approaches to identify the candidate genes associated with growth-related traits in pufferfish. In the selection signatures, the genes involved in processes or pathways related to growth and metabolism were subjected to analysis using the F_st and π ratio as described above. In total, 106 candidate genes were related to growth and metabolism (Table S13), including 56 genes that were independently identified using F_st, while the other 50 genes were identified using the π ratio, and 10 genes intersected. For example, wnt2b, ccnd2, map2k4, ndst2, stat5b, atp6v1b2, and cdk1 were identified using F_st; ngfr, smad4, rnf31, and star were identified using the π ratio; and itgav, ighv3-43, ighm, dsp, dsg2, dsc2, and atp6v1b2 were identified via intersection. Most of these genes were associated with GO terms or KEGG pathways. Notably, wnt2b, ngfr, itgav, and atp6v1b2 have been associated with cellular growth and energy metabolism in the published literature. In the GWAS, a total of 254 candidate genes related growth-related traits were identified, including rnf213, ngfr, smad4, wnt3, wnt5b, and dsp, which are related to metabolism, growth, and immunity (Table S12).

Ultimately, through intersection analysis, a total of 13 growth-related regions and 13 significant potential candidate genes (itgav, ighv3-43, ighm, atp6v1b2, pld1, xmrk, inhba, dsp, dsg2, and dsc2,

Table 3. Growth-related genes identified using both selection signatures and GWAS.

Chr	Gene	Symbol	Description	Selected Region	Fst	π Ratio	SNPs Number	Traits	p Value
1	LOC101079631	itgav	integrin alpha-V	28,970,001–29,070,000	0.71	1.82	-	-	-
5	LOC101067887	ighv3-43	immunoglobulin epsilon heavy chain-like	6,090,001–6,220,000	0.86	1.78	-	-	-
5	LOC101075124	ighm	Ig mu chain C region membrane-bound form	6,120,001–6,200,000	0.86	1.78	-	-	-
-	LOC101072524	atp6v1b2	V-type proton ATPase subunit B, brain isoform	390,001–490,000	0.89	4.45	-	-	-
10	LOC101069674	pld1	phospholipase D1-like	11,240,001–11,340,000	0.35	1.07	6	BL	8.48 × 10−5
10	LOC101073748	xmrk	melanoma receptor tyrosine-protein kinase-like	11,880,001–12,000,000	0.30	1.23	1	TL	6.99 × 10−5
10	LOC101073973	inhba	inhibin beta A chain-like	11,870,001–12,000,000	0.29	1.24	1	TL	6.99 × 10−5
22	LOC115247960	dsp	desmoplakin-like, partial	1,290,001–1,490,000	0.58	1.98	2	TL	5.35 × 10−6
22	LOC101063843	dsp	desmoplakin-like, partial	1,380,001–1,490,000	0.64	1.09	2	TL	5.35 × 10−6
22	LOC115247930	dsg2	desmoglein-2-like, partial	1,340,001–1,49,0000	0.66	1.47	2	TL	5.35 × 10−6
22	LOC101064287	dsg2	desmoglein-2-like, partial	1,250,001–1,450,000	0.57	2.03	2	TL	5.35 × 10−6
22	LOC101064516	dsc2	desmocollin-2-like	1,240,001–1,440,000	0.57	2.04	2	TL	5.35 × 10−6
22	LOC115247939	dsc2	desmocollin-2-like	1,330,001–1,490,000	0.68	1.74	2	TL	5.35 × 10−6

BL, body length; TL, total length.

, Figure S9) associated with pufferfish growth were identified.

4. Discussion

4.1. The Growth Rate and Model

Growth is one of the most significant traits for the aquaculture industry, as it is directly related to economic efficiency [5]. At the same time, the growth rate is a crucial factor that optimizes feeding rates and breeding density. Hybrid puffer has a significant growth advantage over obscure puffer, but the growth rate and model were until unknown. Constructing a growth model to estimate or predict fish weight is essential for achieving better control of the output of an aquaculture system [6], and growth is influenced by several factors. Therefore, we controlled for factors that influence growth, such as environmental conditions, feeding practices, and density. The results show that there is a significant difference in growth and morphological traits between the hybrid and obscure puffers, which might be affected by hybrid dominance. Similar heterosis results have been reported for other aquatic species, such as snow trout [13] and groupers [11].

The growth model is a valuable reference for the aquaculture industry [27]. However, most studies on pufferfish have focused on early-stage and short-term growth models. For instance, previous research has revealed that the early developmental model exhibits negative allometric growth for T. rubripes (b < 3) [27], while the juvenile developmental model demonstrates uniform growth for the obscure puffer (b ≈ 3) [52], indicating a lack of investigation into long-term growth models. Notably, this study shows that the obscure puffer has a negative allometric growth model and the hybrid puffer has a uniform growth model. These novel growth models imply regulatory mechanisms biased towards patrilineal inheritance mechanisms; similar results were found regarding growth rates among catfish (Silurus lanzhouensis) and Pelteobagrus fulvidraco (♀) × Leiocassis longirostris (♂) hybrids [12,53].

4.2. The Growth-Related Candidate Genes and SNPs

Growth is a complex quantitative trait regulated by multiple genes [54]. The analysis of selection signatures is a genomic analysis method that was used in this study to identify the specific regions of higher genetic differentiation between hybrid and obscure puffers, and thereby to elucidate the molecular mechanisms underlying hybridization. Similarly, GWAS constitutes a pivotal approach in aquaculture research, especially in SNP trait association research, such as growth, disease resistance, color differentiation, and hypoxia tolerance [23,55,56,57]. However, both methods have limitations; therefore, the objective of this study was to combine the use of selection signatures and GWAS to investigate the growth-related selection regions and genes in pufferfish.

In the selection signatures, the results demonstrated a close association between the candidate genes derived from the identified selection features, several crucial growth- and metabolism-related GO terms, and KEGG pathways (Figure 8 and Figure 9). Remarkably, the regulation of muscle contraction (GO:0006937) inhibits epidermal growth factor receptor (EGFR)–extracellular signal-regulated kinase (ERK)1/2 signaling under EGF stimulation, which can also enhance growth ability [58], and cellular senescence (ko04218) regulates cellular aging and damage processes [59]. All these pathways are implicated in cellular growth, development, and programmed cell death through indirect mechanisms, highlighting the intricate biological complexity underlying growth traits. In the GWAS, several SNPs were associated with multiple traits. Similar results were found for mandarin fish and Pacific abalone [23], indicating that coordinated regulatory mechanisms may be involved. We hypothesized that growth-related traits in pufferfish are influenced by the coordinated regulation of multiple genes across different chromosomes. Similar observations have been reported in other fish species, such as the Asian sea bass, wherein Lca371 on LG2 was found to be involved in the regulation of TL, BW, and BL [60]. In carp, SNP0626 on LG19 was significantly associated with BW and BL [61].

We identified several genes associated with growth-related traits according to the selection signatures and GWAS, including itgav, dsp, dsg2, dsc2, rnf213, smad4, atp6v1b2, map2k4, stat5b, and inhba. Several of these genes were associated with fish growth in previous studies [21,23]. The itgav gene functions as a major player in regulating cell proliferation, apoptosis, migration, maintaining the integrity and permeability of vascular wall, and other biological functions, and it is related to the growth of osteosarcoma tissue [62]. In zebrafish, the knockout of rnf213 was associated with abnormal vascular development of the head [63], and it was demonstrated that it is associated with fat formation and highly expressed in fat samples from obese populations [64]. Meanwhile, these genes regulate cellular proliferation and apoptosis through complex regulatory networks, such as the Wnt, MAPK, and TGF-β signaling pathways and other synergistic regulatory mechanisms, affecting the growth traits of individuals. For example, dsp has been shown to be a downstream regulator of the arrhythmogenic right ventricular cardiomyopathy (ARVC) pathway (ko05412) and to carry out inhibitory regulation of the Wnt pathway [65], in turn affecting its downstream regulation of the cell-cycle protein (ccnd2), thereby affecting cellular value creation [66]. In the Wnt signaling pathway, wnt3 is upstream of the regulatory genes TCF/LEF, which also regulate TCF/LEF and cellular value-added through ccnd2 [67]. Further, in the TGF-β signaling pathway, smad4 can not only directly regulate downstream gonadal development and apoptosis [68], but has been shown to have a regulatory effect on the cell growth cycle through the Wnt signaling pathway [69,70]. Similarly, ngfr and map2k4, as key genes in the MAPK signaling pathway, regulate the cell cycle, growth, and death [71]. Studies have shown that map2k4 regulates cell growth by activating cell cycle-associated proteins and EMT signals downstream of the PI3K–Akt signaling pathway [72]. Taken together, these candidate genes do not act individually to regulate growth but function through multiple pathways for synergistic network regulation. Although the experimental group used in this study was not a traditional group, such as a full-sib family, but rather a comparison between the hybrid puffer and the parental obscure puffer, these were used because the use of different individuals for the analyses is still of interest for how to select the parental individuals in breeding research work, which provides a theoretical basis for the subsequent study of fish growth and development.

In this study, we identified valuable SNPs and growth-related candidate genes that could contribute to the development of pufferfish growth models. Our findings may be applicable to pufferfish farming and may significantly accelerate the breeding process. However, other relevant physiological regulations in pufferfish breeding and genetics still need to be studied.

5. Conclusions

Our research indicates that the growth of hybrid puffers is considerably superior to that of obscure puffers, and there are notable disparities in their morphological traits. Furthermore, the growth and morphological traits of hybrid puffers are concurrently influenced by the genetic effects of the parental obscure puffers. Growth models of hybrid and obscure puffers in the culture cycle were constructed in this study based on their morphological and growth traits, revealing that there are different growth models between the obscure and hybrid puffers. The obscure puffer exhibits a negative allometric growth model, while the hybrid puffer exhibits a uniform growth model. After integrating selection signatures and GWAS analysis, a total of 13 growth-related regions and 13 significant potential candidate genes (itgav, ighv3-43, ighm, atp6v1b2, pld1, xmrk, inhba, dsp, dsg2, and dsc2) associated with pufferfish growth were identified.