Evidence for the Functional Relevance of vgll3 and six6 οn Developmental Stages of Commercially Important Fish Species: Gilthead Seabream (Sparus aurata Linnaeus) and European Seabass (Dicentrarchus labrax Linnaeus)

1. Introduction

The SIX homeobox 6 (six6) and vestigial-like family member 3 (vgll3) genes function in an evolutionarily conserved manner in developmental processes and life-history traits across a wide range of taxa, including mammals and teleosts [1,2,3,4,5,6,7,8]. The six6 gene encodes a transcription factor involved in sexual maturation [3,7,8,9]. A locus near the six6 gene has been linked to age at maturity and spawning ecotype in different salmonid species [7,8]. Additionally, six6 has been associated with maturation in mammals, specifically age-at-menarche in humans [2] and puberty timing in beef cattle [3]. The vgll3 gene encodes a transcriptional cofactor involved in adipogenesis and the Hippo signaling pathway [9,10,11] that has been linked to pubertal timing and growth in humans [1]. In Atlantic salmon (Salmo salar Linnaeus), a single nucleotide polymorphism (SNP) in the vgll3 gene correlates with maturation timing [5,6], with different genotypes of this SNP influencing the expression profiles of reproductive axis genes in juvenile males [12]. It was further demonstrated that vgll3 influences the body condition of the species [13,14]. Both six6 and vgll3 are expressed during early development in salmon and are associated with adipogenesis and gametogenesis [9]. This growing body of evidence highlights that six6 and vgll3 play key roles not only in reproductive traits but also in growth and growth-related processes. In aquaculture, growth rate is a critical trait for productivity, and understanding how these genes regulate developmental processes may provide insights into improving growth performance in commercially farmed species [15,16,17]. By leveraging candidate genes identified in one species and assuming their evolutionarily conserved roles, we can investigate their potential for similar effects in other species [8]. It remains an open question whether six6 and vgll3 also influence developmental processes and life-history traits in other economically important teleost fish.

The gilthead seabream (Sparus aurata Linnaeus) and the European seabass (Dicentrarchus labrax Linnaeus) are among the most economically important fish species in the Mediterranean region. A recent study on genetic differentiation between farmed and wild populations of these species investigated how the domestication process impacts the genetic diversity of the six6 and vgll3 genes [18]. Significant differences in genotype and allele frequencies were observed in an exonic SNP on vgll3 and in two intronic SNPs on six6 between farmed and wild populations of both species [18]. Given this evidence and the high level of functional conservation of vgll3 and six6, further investigations, including gene expression profiling across various developmental stages in these species, could yield highly informative insights. Such studies would be particularly valuable for the aquaculture sector and future selective breeding programs targeting maturation and growth performance [8,18,19]. Research in this direction could also enhance our understanding of the genotype- phenotype link and shed light on the genomic architecture of domestication.

To investigate the functional roles of six6 and vgll3 in the developmental processes of European seabass and gilthead seabream, we performed a gene expression analysis using quantitative real-time PCR (qPCR) in early developmental stages. Gene expression studies are essential for understanding the functional impact of genetic variants and for bridging the gap between genotype and phenotype [20,21]. By examining six6 expression across genotypes and developmental stages (larvae and juveniles) in European seabass, and vgll3 and six6 expression patterns across genotypes in gilthead seabream juveniles, we aimed to explore the association of previously identified domestication-related genes with their gene expression patterns [18]. This study sought to provide deeper insights into the roles of these genes in the developmental processes of European seabass and gilthead seabream, as well as their potential evolutionarily conserved functions.

2. Materials and Methods

2.1. Sampling

Fifty larvae at 34 days post-hatch (dph) and 50 juveniles at 71 dph of European seabass, along with 50 juveniles at 69 dph of gilthead seabream, were collected from a fish farm. Both the larval and juvenile stages represent critical time points in the developmental process of fish. All individuals were reared under the same controlled conditions, which included a temperature of 19°C, salinity of 40‰, a pH range of 7.6 to 7.8, and O₂ saturation levels of 6-7 ppm. Samples were preserved in RNAlater stabilization solution (Invitrogen) before being frozen at -20°C.

2.2. Genotyping

DNA extractions from each individual were carried out using the TRIzol reagent (Thermo Fisher Scientific) [22]. The primer sequences, amplicons sizes, and annealing temperatures for each studied genomic region are listed in Table 1 [18]. For PCR, we used a total reaction volume of 25 μL, which included 100 ng of genomic DNA as the template, 0.05 units of Qiagen Taq polymerase, 2 mM dNTPs, 0.25 μL of each primer (100 μΜ), and 2.5 μL of 10x Reaction Buffer (Qiagen, Hilden, Germany). The PCR products were analysed by electrophoresis in 1.5% (w/v) agarose gels. Enzymatic cleanup and Sanger sequencing were outsourced to Genewiz (Leipzig, Germany). SNP genotyping was performed by aligning the resulting sequences with the corresponding reference sequences from GenBank (assemblies S. aurata: GCA_900880675.1; D. labrax: GCA_000689215.1) using the Geneious v.10.2.6 program [23].

2.3. RNA Extraction and cDNA Synthesis

Total RNA was extracted from whole European seabass and gilthead seabream individuals using the TRIzol reagent (Thermo Fisher Scientific, Waltham, USA) [22]. In total, 27 European seabass individuals (nine biological replicates per SNP genotype) and 24 gilthead seabream individuals (eight biological replicates per SNP genotype) were selected for analysis. The quality and quantity of RNA were assessed using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA) and by electrophoresis on 1.5% (w/v) agarose gels. The cDNA synthesis was performed using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan) following the manufacturer’s instructions.

2.4. Quantitative Real-Time PCR

We selected ribosomal protein L13 (rpl13), FAU ubiquitin-like and ribosomal protein S30 fusion (fau), and elongation factor 1-alpha 1 (ef1α) as candidate reference genes for both species. Primers and probes for the reference and target genes were designed using the qPCR Primer & Probe Design Tool (Eurofins, Luxembourg City, Luxembourg) (Table 2). Gene locations and intron/exon information were retrieved from annotations deposited in Ensembl for both species (assemblies S aurata: GCA_900880675.1; D. labrax: GCA_000689215.1). Gene expression profiles were quantified using quantitative real-time PCR with the KAPA^® PROBE FAST Universal reagent (Kapabiosystems, Cape Town, South Africa) in three technical replicates per sample, performed on the ABI StepOne™ platform (Applied Biosystems, CA, USA).

2.5. Analysis of Gene Expression Data

The Ct values were converted to quantities using the comparative Ct method. The most stable reference genes were identified using the NormFinder and geNorm algorithms, which ranked the genes based on their expression stability [24,25]. The normalization factor was calculated as the geometric mean of the two most stable reference genes for each species (rpl13 and fau in European seabass; rpl13 and ef1α in gilthead seabream). The normalized expression levels of the target genes were determined by dividing the raw target gene quantities for each sample by the corresponding normalization factor. The Kruskal-Wallis test was used to compare the normalized expression levels of the target genes among the different SNP genotypes [26]. Additionally, the distribution of the expression data was assessed using the Shapiro-Wilk normality test, followed by the Bartlett’s test for homogeneity of variances [27,28].

3. Results

3.1. six6 and vgll3 Genotypes

We analysed independent datasets from two developmental stages (larvae at 34 dph and juveniles at 71 dph) in European seabass and one developmental stage (juveniles at 69 dph) in gilthead seabream. For both species, we genotyped two intronic SNPs of six6, previously associated with the domestication process, and selected the SNPs with sufficient sample sizes for each genotype (S. aurata 16:14406442 and D. labrax 12:11591053) [18]. No genetic variation in vgll3 was observed in the studied European seabass individuals; thus, the expression profile per genotype was not examined. In contrast, all genotypes for the vgll3 SNP (S. aurata 9:24911884) were identified in the studied juveniles of gilthead seabream [18]. The sample sizes for each studied genotyped SNP are shown in Table 3.

3.2. Expression Patterns of vgll3 and six6 per Genotype in Gilthead Seabream

In gilthead seabream, six6 and vgll3 expression levels were analysed in juveniles (69 dph). For vgll3, the expression levels were highest in the AA genotype, followed by AG, with the GG genotype showing the lowest expression levels (Figure 1a). The differential expression across genotypes was statistically significant (p = 0.040), and pairwise comparisons revealed that the AA genotype was significantly different from the other two genotypes (AA-AG: p = 0.030; AA-GG: p = 0.027). Variance differences were also significant across genotypes (p = 0.002), as well as between AA-AG (p = 0.016) and AA-GG (p = 0.003). A similar trend was observed for six6, where the CC genotype exhibited the highest expression levels, followed by CT, with the TT genotype displaying the lowest (Figure 1b). A significant differential expression was detected between the CT and TT genotypes (p = 0.027), but no significant differences in variability were observed.

3.3. Expression Patterns of six6 per Genotype and Developmental Stage in European Seabass

All observed trends in expression levels among genotypes at both developmental stages were not statistically significant. However, in larvae at 34 dph, the heterozygous genotype (AT) exhibited high variability compared with homozygous genotypes (AT-AA: p = 0.040; AT-TT: p = 0.046; Figure 2a). At the 71 dph juveniles, this pattern had disappeared with all genotypes exhibiting a similar range of expression levels (Figure 2b). Overall, the variation in six6 expression for the heterozygous genotype showed a decreasing trend between larvae and juveniles (p = 0.030).

4. Discussion

In this study, we investigated the gene expression profiles of vgll3 and six6 in the early developmental stages of two commercially important fish species, the gilthead seabream and the European seabass. We aimed to provide the first evidence that prior SNP associations of these genes with domestication [18] are also reflected at the mRNA level. Gene expression profiling was performed using quantitative real-time PCR (qPCR) with probes, due to its high accuracy, low background signal [29], and efficiency [30]. We identified differences in both the means and variances of mRNA expression across genotypes of vgll3 and six6 in gilthead seabream at the developmental stage examined (69 dph juveniles). Additionally, significant variability in six6 expression levels was observed in the heterozygous genotype during the larval stage (34 dph), which was reduced in juveniles (71 dph) in European seabass. These observations support a functional role for vgll3 and six6 in domestication and developmental processes in the studied species. Given the suggested roles of these genes in growth and growth-related processes [13,31,32,33], our findings further underscore their potential relevance to the aquaculture sector. These insights open up new opportunities to enhance desirable traits in cultured populations of the species, ultimately benefiting the aquaculture industry [34,35].

Our findings indicate that gene expression profiles are influenced by both genetic variation (Figure 1) [36,37,38] and developmental stage (Figure 2) [9,39,40]. Specifically, we report a statistically significant difference in expression levels between the two homozygous genotypes (AA and GG) of vgll3 in gilthead seabream (Figure 1a). The AA genotype exhibited significantly higher expression levels and variance compared to the GG genotype, which was associated with domestication in previous DNA-based analysis [18]. A similar pattern was observed in juveniles of Atlantic salmon, Verta et al. (2020) previously found that the late-maturing vgll3 genotype exhibited 14% higher expression and variation within the testis compared to the early-maturing genotype [36]. Beyond their association with maturation, the vgll3 genotypes influenced body condition, with the early-maturing genotype exhibiting higher body condition [13] and seasonally increasing lipid storage in immature Atlantic salmon [37]. The consistent in vgll3 expression patterns among genotypes between juveniles gilthead seabream and Atlantic salmon prompt a key question about the potential role of vgll3 genotypes in physiological processes in gilthead seabream.

The six6 differential expression was statistically significant between the CT and TT genotypes in gilthead seabream, with the TT genotype exhibiting the lowest expression levels compared to the other two genotypes (Figure 1b). Notably, the TT genotype was associated with domestication in a previous DNA-based analysis using farmed and wild populations [18]. Additionally, in Atlantic salmon, six6 expression patterns have been associated with the early developmental process, particularly in sensory, musculoskeletal, and cardiovascular systems [9,33]. However, there is a lack of literature on the effects of six6 genotypes on expression levels and whole transcriptome profiles. The observed differential expression among genotypes in this study suggests a potential role of six6 expression in the juvenile stage in gilthead seabream, warranting further investigation to clarify its specific functions and regulatory mechanisms.

Observations of six6 in European seabass also indicate its potential significance in this species, with the variation in expression levels for heterozygous genotype in larvae being higher than those observed in homozygous genotypes. The expression variance may suggest an association with the molecular basis of complex traits, diseases [41,42], and environmental adaptation mechanisms [43]. Early development is a highly dynamic period in an organism’s life, during which constant changes in gene expression are crucial for regulating processes that govern growth, cell differentiation, organogenesis, and the physiological adaptations of each stage [9,39,40]. A decreasing trend in six6 expression was also noted from the embryonic stages toward fry in Atlantic salmon, highlighting the importance of the gene’s role during the early stages for the hypothalamic-pituitary-gonadal (HPG) axis and eye development [9]. These findings collectively underscore the potential role of six6 in early developmental processes across multiple fish species, warranting further investigation to uncover its specific roles in them.

It remains unclear how the studied SNPs are associated with the observed patterns in DNA and expression levels. The vgll3 SNP in gilthead seabream is a synonymous mutation and the six6 SNPs are intronic in both species [18]. The literature provides evidence that synonymous and intronic mutations can alter the regulation of gene expression and potentially play a role in fish domestication and traits linked to aquaculture [44,45,46]. Synonymous mutations do not alter the protein sequence, but they may affect mRNA stability, translation efficiency, and protein folding [47,48]. Intronic variation can also affect gene expression levels and phenotype by modifying splice site activity, altering cis-regulatory elements like enhancers and silencers, and changing the formation of chromatin loops [49,50,51,52]. Furthermore, it is recognized that alternative splicing contributes to transcriptional diversity in genes associated with developmental processes [53], such as those under study. However, the presence of these SNPs does not necessarily mean they are the causative variants affecting gene expression, as they may be in linkage disequilibrium with the actual causative variant. Further studies are needed to determine if the SNPs have a functional effect or are merely markers for a linked variant.

Our study focused on the expression of targeted genes at specific developmental stages (larvae and juveniles), providing important insights into the roles of vgll3 and six6 in gilthead seabream and European seabass. Notably, research on Atlantic salmon has identified these stages as relevant for studying the expression of these genes [9]. Future research should include more developmental stages and omic-wide approaches, such as genomic and transcriptomics, to confirm our present observations.

5. Conclusions

In conclusion, this study provides new insights into the expression profiles of six6 and vgll3 in developmental stages in gilthead seabream and European seabass, linking these gene expression patterns to domestication-associated SNPs previously identified in DNA-based study. These findings, combined with the evolutionarily conserved functions of these genes and their involvement in developmental processes and life history traits across vertebrates, including teleosts, highlighting the potential roles of these genes in growth, development, and adaptation to domesticated environments. Particularly, the expression differences observed in vgll3 and six6 between gilthead seabream genotypes and the variance of six6 expression between larval and juvenile stages in European seabass, support the notion that these genes influence growth-related processes and physiological traits relevant to aquaculture. Future studies expanding to additional developmental stages and incorporating transcriptome-wide approaches could validate these observations and enhance our understanding of the genetic underpinnings of domestication and growth in cultured fish, ultimately contributing to improved breeding strategies in aquaculture.