Genotype‐Dependent Transcriptome Divergence Associated With Variation at vgll3 in Juvenile Gilthead Seabream (Sparus aurata)
Aristotelis Moulistanos, Elisavet Kaitetzidou, Styliani Minoudi, Konstantinos Gkagkavouzis, Efthimia Antonopoulou, Alexandros Triantafyllidis, Spiros Papakostas
TL;DR
This study shows how a genetic variant in the vgll3 gene affects gene expression in juvenile gilthead seabream, influencing growth and maturation traits important for aquaculture.
Contribution
The study reveals genotype-dependent transcriptome divergence linked to the vgll3 gene in gilthead seabream, highlighting its role in growth and maturation.
Findings
The GG genotype of vgll3 is associated with distinct transcriptome profiles and lower vgll3 expression.
Genes like amh, cacng1b, and casq2 show differential expression linked to vgll3 genotype differences.
The findings suggest a conserved role of vgll3 in growth and reproductive maturation pathways in fish.
Abstract
Early developmental processes significantly influence growth and maturation patterns, aquaculture traits that are critical for physiological adaptation and productivity. The vestigial‐like family member 3 gene (vgll3) plays a key role in growth and maturation across diverse taxa, including mammals and teleost fishes. A single‐nucleotide polymorphism in vgll3 (SNP vgll3 ) shows evidence of selection under aquaculture conditions in gilthead seabream ( Sparus aurata ), as demonstrated by previous genome scan and targeted transcriptomic (qPCR) analyses. This study investigated how different SNP vgll3 genotypes (AA, AG and GG) affect gene expression in juvenile gilthead seabream. Genotype‐dependent regulatory signatures were identified, as the transcriptome profiles (over 240 quantified transcripts) of the farming‐associated GG genotype, which also showed significantly lower vgll3…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3| Gene name ( | Ortholog name ( |
| log2 fold change |
|---|---|---|---|
| LOC115588759 |
| 2.381e‐18 | −1.65 |
| LOC115586181 |
| 6.592e‐12 | −3.49 |
| LOC115583046 |
| 3.482e‐11 | −1.83 |
| LOC115587310 |
| 7.221e‐11 | −1.54 |
| LOC115568588 |
| 5.000e‐10 | −1.40 |
|
|
| 8.973e‐9 | −1.49 |
|
|
| 5.703e‐6 | −1.55 |
| LOC115585027 |
| 7.610e‐6 | −1.53 |
|
|
| 9.118e‐6 | −1.67 |
| Term name |
| Term size | Number of identified genes |
|---|---|---|---|
|
| 5.476e‐6 | 58 | 8 |
|
| 1.073e‐3 | 102 | 8 |
|
| 1.664e‐3 | 89 | 7 |
|
| 1.696e‐3 | 130 | 8 |
|
| 1.696e‐3 | 2517 | 40 |
|
| 3.572e‐6 | 84 | 9 |
|
| 2.538e‐5 | 20 | 5 |
|
| 2.137e−3 | 172 | 10 |
| Term name |
| Term size | Number of identified genes |
|---|---|---|---|
|
| 4.814e‐4 | 833 | 5 ( |
|
| 1.120e‐6 | 32 | 3 ( |
|
| 2.295e‐7 | 73 | 4 ( |
|
| 2.189e‐4 | 122 | 3 ( |
|
| 2.040e−6 | 172 | 4 ( |
|
| 4.920e‐5 | 113 | 3 ( |
- —Hellenic Foundation for Research & Innovation (H.F.R.I.)10.13039/501100013209
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsReproductive biology and impacts on aquatic species · Developmental Biology and Gene Regulation · Aquaculture Nutrition and Growth
Introduction
1
Understanding the molecular mechanisms underlying phenotypic diversity across species and environments remains a central challenge in biology. Although protein‐coding mutations have been traditionally recognised for altering gene function, accumulating evidence emphasises that regulatory variation can significantly influence physiological processes, shaping organismal fitness and phenotypic traits (Ahi et al. 2022; Ahi, Verta, et al. 2025a; Albert and Kruglyak 2015; Jiang et al. 2020; Moulistanos et al. 2024; Verta et al. 2020). Identifying such regulatory variants is challenging due to their subtle effects and the intricate interactions among non‐coding elements within complex regulatory networks (Paul et al. 2014; Rojano et al. 2019). In this context, transcription factors are critical as integrators of regulatory signals, orchestrating gene expression patterns that determine phenotypic outcomes (Weidemüller et al. 2021).
Research in commercially important fish species has contributed substantially to our understanding of gene regulation and phenotypic diversity. Comparative transcriptomic studies in Atlantic salmon ( Salmo salar ) and Nile tilapia ( Oreochromis niloticus ) have uncovered conserved molecular responses to environmental and physiological challenges, advancing core biological knowledge well beyond immediate aquaculture applications (Ellison et al. 2018, 2020; Beemelmanns et al. 2021; Krasnov et al. 2021). Transcriptome profiling in farmed salmon has revealed regulatory networks underlying metabolism, immunity, and development (Jiang et al. 2020), while genomic analyses in intensively farmed teleosts have identified extensive genetic variation and regulatory elements associated with key life‐history traits (Sinclair‐Waters et al. 2020; Moulistanos et al. 2024). Collectively, these systems have emerged as powerful comparative models, broadening the phylogenetic and mechanistic scope of molecular ecology, and providing fundamental insights comparable to those derived from classical model organisms.
In gilthead seabream ( Sparus aurata Linnaeus, 1758), elucidating the molecular factors controlling development and growth is crucial within the framework of life‐history theory, in which traits are defined by the timing of key developmental transitions such as somatic growth, maturation, and sex differentiation (Chavanne et al. 2016; Teletchea 2021). Genetic and physiological variation influencing these transitions generates life‐history trade‐offs, such as the allocation of resources between growth and reproduction, variation in age and size at maturity, and alternative reproductive pathways. Growth‐related traits in gilthead seabream emerge from the interplay between environmental, genetic, and physiological factors, each contributing to the complex regulation of muscle development—a major determinant of body size and mass (Birnie‐Gauvin et al. 2021; Mohammadabadi et al. 2021; Moya et al. 2019; Rodgers and Gomez Isaza 2024). In this context, aquaculture represents a powerful and highly directional selective regime, in which growth‐oriented management and controlled rearing conditions can systematically alter developmental timing and resource allocation, thereby reshaping life‐history trade‐offs practice (Chavanne et al. 2016; Janssen et al. 2017; Milla et al. 2021; Teletchea 2021). Such conditions provide a unique opportunity to investigate how regulatory genetic variation influences the molecular mechanisms underlying these trade‐offs (Roff 2007).
A particularly promising candidate gene implicated in these processes is the vestigial‐like family member 3 (vgll3). vgll3 encodes a transcriptional cofactor that regulates various life‐history traits, including growth and maturation, across diverse taxa, underscoring its evolutionarily conserved role. In mammals, vgll3 has been associated with pubertal height growth and age at menarche in humans (Cousminer et al. 2013; Perry et al. 2014), as well as body weight regulation and gonadal adipose content in mice (Halperin et al. 2013). In Atlantic salmon ( Salmo salar Linnaeus, 1758), a single nucleotide polymorphism (SNP) in vgll3 correlates with body condition (Debes et al. 2021), alongside its well‐documented role in maturation timing (Barson et al. 2015). In non‐salmonid species such as zebrafish ( Danio rerio , Hamilton, 1822), vgll3 has been implicated in growth during early development (Pennonen 2017). Recent studies in gilthead seabream and European seabass ( Dicentrarchus labrax Linnaeus, 1758) suggest that vgll3 may represent a putative farming‐associated locus, reinforcing its conserved function in determining growth regulation among aquaculture species (Moulistanos et al. 2023).
Importantly, vgll3 is a well‐characterised transcriptional cofactor within the Hippo signalling pathway, a deeply conserved regulatory cascade that integrates signals governing cell proliferation, organ growth, energy balance and environmental sensing. VGLL3 interacts with TEA domain (TEAD) transcription factors to regulate a broad spectrum of downstream genes (Hori et al. 2020; Mesrouze et al. 2020; Mesrouze and Chène 2024), thereby influencing developmental timing and tissue growth. Across vertebrates, the Hippo pathway has emerged as a central regulator of life‐history traits, including growth rate, maturation timing, and reproductive investment traits that are central to pace‐of‐life variation (Ahi, Panda, and Primmer 2025). Recent work positions vgll3 as a critical mediator linking Hippo signalling to maturation and growth decisions, particularly when environmental conditions shape developmental trajectories (Barson et al. 2015; Ahi, Verta, et al. 2025a, 2025, 2025b; Halperin et al. 2013; Debes et al. 2021). Because the pathway responds to cues such as energy availability, nutritional status and population density, variation at vgll3 can have pronounced effects under aquaculture conditions, where growth‐focused selection and rearing environments interact to influence life‐history trade‐offs. Situating vgll3 within the Hippo signalling network provides a mechanistic basis for interpreting genotype‐specific transcriptional differences and highlights how variation at this locus can propagate through broader regulatory networks rather than act in isolation.
Building upon these findings, a synonymous SNP within vgll3 (SNP_ vgll3 : AA, AG and GG) in gilthead seabream shows marked differences in genotype and allele frequencies between wild and farmed populations (Moulistanos et al. 2023). In Greek populations, the G allele is more frequent in farmed fish, and analyses across the Mediterranean indicate that this pattern is consistent across multiple wild and farmed populations, spanning genetically differentiated populations rather than being restricted to a single locality (Moulistanos et al. 2023; Peñaloza et al. 2021; Moulistanos 2025). Although the degree of independence among farmed populations is not fully known, and some may trace back to shared hatcheries, these patterns support the G allele as a candidate for genotype‐specific divergence under aquaculture conditions. This positions vgll3 as a compelling target for investigating molecular signatures associated with putative domestication‐related traits, without assuming strong selection in all farmed populations. A qPCR analysis revealed significant differences in vgll3 expression among SNP vgll3 _ genotypes in juvenile gilthead seabream, with the GG genotype, the most frequent genotype in farmed populations, showing lower expression (Moulistanos et al. 2025). Interestingly, in Atlantic salmon, lower vgll3 expression is associated with a higher probability of male maturation and increased fish length during the first year of development (Verta et al. 2020). Although these cross‐species parallels are intriguing, the gilthead seabream data, derived from whole‐body juvenile samples at a single time point, provide only limited support for functional conservation of vgll3 between species.
In this study, we conducted a comprehensive genotype‐specific transcriptome analysis of juvenile gilthead seabream at 69 days post‐hatch (dph) to investigate the molecular consequences of the farming‐associated SNP in vgll3 (SNP_ vgll3 _). We analysed four individuals per genotype (AA, AG and GG) previously characterised through qPCR for distinct vgll3 expression levels (Moulistanos et al. 2025). Our primary goal was to identify and characterise genotype‐specific transcriptional profiles through differential expression analyses, thereby elucidating the biological processes, molecular functions, and regulatory pathways associated with vgll3 variation. We hypothesised that distinct vgll3 genotypes would be accompanied by broader transcriptomic shifts, reflecting altered developmental programs potentially shaped by farming‐related selection pressures. Ultimately, we aimed to move beyond single‐gene analyses to uncover extensive regulatory networks underpinning traits of biological, evolutionary, and economic significance.
Materials and Methods
2
Sampling
2.1
In total, 50 gilthead seabream juveniles, all from a single cohort and raised in the same aquaculture facility, were collected at 69 days post‐hatch (dph). Detailed pedigree information, including family structure and number of broodstock parents, was not available due to restrictions imposed by the collaborating hatchery. This limits our ability to assess potential relatedness effects. Nevertheless, all sampled individuals were reared in a single cohort under identical environmental conditions, minimising environmental variance as a confounding factor in the observed transcriptional differences. Juveniles were kept under controlled conditions, with a temperature of 19°C, salinity of 40‰, pH range of 7.6 to 7.8, and O_2_ saturation levels of 6–7 ppm. After sampling, the juveniles were immediately preserved in RNAlater stabilisation solution (Invitrogen, Waltham, USA) and were stored at −20°C until further processing.
Genotyping
2.2
Single‐individual DNA extraction was performed using the TRIzol reagent (Thermo Fisher Scientific) (Rio et al. 2010). The extracted DNA was then used for PCR amplification of the vgll3 SNP region using the primers 5′AACGTCTATCACCCTCACCC3′ and 5′ACCAAACTGACGTCTTTGCT3′ (Moulistanos et al. 2023). The total reaction volume of 25 μL consisted of 100 ng of genomic DNA as template, 0.05 units of Qiagen Taq polymerase, 2 mM dNTPs, and 1 μM of each primer and 2.5 μL of 10× Reaction Buffer (Qiagen, Hilden, Germany). The PCR conditions were as follows: initial denaturation at 94°C for 2 min, 35 cycles of denaturation at 94°C for 30 s, annealing at 63°C for 40 s and extension at 72°C for 1 min, and a final extension at 72°C for 10 min. The PCR products were evaluated by electrophoresis in 1.5% (w/v) agarose gels and were afterwards outsourced for enzymatic cleanup and Sanger sequencing (Genewiz Company, Leipzig, Germany). Genotyping of SNP_ vgll3 _ ( S. aurata 9:24911884) was performed by aligning the resulting sequences with the reference sequence from GenBank (assembly S. aurata : GCA_900880675.1), using the Geneious v.10.2.6 program (Kearse et al. 2012).
RNA Extraction and Sequencing
2.3
Total RNA was extracted from whole juvenile gilthead seabreams using TRIzol reagent (Thermo Fisher Scientific, Waltham, USA) (Rio et al. 2010). In total, 12 individuals (four biological replicates per SNP_vgll3_ genotype: AA, AG, GG) previously genotyped and characterised via qPCR for vgll3 expression by Moulistanos et al. (2025), were selected for mRNA sequencing. RNA quality and integrity were confirmed by electrophoresis in 1.5% agarose gels and spectrophotometric quantification (Nanodrop 2000, Thermo Fisher Scientific, Waltham, USA). The RNA samples were outsourced (Novogene, Europe, UK) for strand‐specific mRNA library preparation using poly‐A selection, followed by paired‐end sequencing (150‐bp reads) on the Illumina NovaSeq platforms.
mRNA Transcriptome Analysis
2.4
Raw RNA‐Seq data were quality trimmed with Trimmomatic v.0.38 (Bolger et al. 2014) to remove low‐quality reads and adapter sequences. Trimming was performed in paired‐end mode with the following parameters: ILLUMINACLIP:Illumina adaptors file.fa:2:30:10; SLIDINGWINDOW:3:20; MINLEN:50. Filtered reads were mapped to the gilthead seabream reference genome (GCF_900880675.1), using HISAT2 (Kim et al. 2019). Genotype‐specific transcriptional profiles were visualised using Principal Component Analysis (PCA), implemented through the R function ‘prcomp’. Gene expression data were log_2_‐transformed and standardised by Z‐score normalisation across samples and were further visualised using hierarchical clustering with Euclidean distance and complete linkage in a heatmap generated by the ‘pheatmap’ R package. Differential gene expression analysis was conducted with the Bioconductor package ‘metaseqR2’ (Fanidis and Moulos 2021; Moulos and Hatzis 2015), with *p‐*values adjusted using the Benjamini and Hochberg method to control the false discovery rate (FDR). Transcripts were considered differentially expressed if the adjusted p‐value (p adj) was below 0.05. All visualisations were generated in R v.4.4.1 (R Core Team 2021).
Functional Annotation and Enrichment Analysis
2.5
Sequences of differentially expressed genes (DEGs) were retrieved from the gilthead seabream reference genome assembly GCF_900880675.1. These sequences were used to identify zebrafish orthologs through local BLAST searches (https://www.ensembl.org/Multi/Tools/Blast). Zebrafish orthologs were subsequently employed for functional annotations and enrichment analyses due to their superior annotation quality and completeness compared to other teleost fish species (Primmer et al. 2013). For each DEG, the highest scoring zebrafish gene was selected as putative ortholog, using an E‐value threshold of 10^−3^.
Gene ontology (GO) and pathway enrichment analyses were conducted using g:Profiler, employing the zebrafish genome as a reference (https://zfin.org/). The Benjamini–Hochberg method was applied to correct for multiple testing (p adj < 0.01). To pinpoint the biologically meaningful (henceforth referred to as ‘overrepresented’) GO terms, the two‐stage algorithm of g:Profiler was utilised (Kolberg et al. 2023). This algorithm effectively reduces redundancy by clustering similar GO terms into biologically meaningful groups, enhancing interpretability. It also highlights the most statistically robust and ‘overrepresented’ terms, providing clearer insights into the underlying biological processes and functions (Kolberg et al. 2023). Promoter sequences of all DEGs were further analysed for transcription factor binding site (TFBS) enrichment using g:Profiler with the TRANSFAC database.
Interactions among zebrafish orthologs were initially analysed based on co‐expression using STRING version 12.0 with default settings (medium confidence = 0.4) (Szklarczyk et al. 2022), with the aim of capturing coordinated transcriptional responses among DEGs and identifying modules of genes that are transcriptionally regulated in a concerted manner. Because co‐expression reflects shared regulatory control and functional coupling at the expression level, this approach provides a data‐driven view of network organisation directly grounded in the observed RNA‐seq profiles. In parallel, a complementary integrative network was constructed by combining co‐expression evidence with experimentally validated and database‐curated interactions. This integrative network extends the co‐expression‐based framework by incorporating known functional relationships, enabling the identification of additional connected genes and potential hub nodes among the DEGs, and providing a broader functional context for interpreting the observed transcriptional modules.
Results
3
All analysed samples yielded high‐quality transcriptome data, with each sample generating over 48 million reads. After quality filtering and alignment to the reference genome, more than 90% of the reads were uniquely mapped. This high mapping efficiency ensured a robust dataset suitable for downstream differential expression and clustering analyses (Supplementary File 1).
The PCA of the four replicate transcriptomes per SNP_ vgll3 _ genotype revealed two distinct groups at the 95% confidence interval, separating the AA and GG genotypes. The transcriptomic profiles of heterozygous AG individuals overlapped substantially with those of the GG genotype (Figure 1a). Hierarchical clustering analysis further supported the PCA results, demonstrating that samples grouped primarily by genotype (Figure 1b). The heatmap showed distinct expression patterns, with AA and GG genotypes forming separate clusters, while AG samples exhibited an intermediate expression profile more closely aligned with the GG genotype. Together, these analyses indicated significant genotype‐dependent differences in transcriptomic profiles of juvenile gilthead seabream.
(a) Principal Component Analysis (PCA) of gene expression data for different gilthead seabream juvenile SNP vgll3 genotypes. Ellipses indicate the 95% confidence intervals for each genotype. (b) Heatmap displaying Z‐score normalised gene expression with hierarchical clustering of samples. Rows represent genes, and columns represent individual samples coloured by genotype. Purple to orange colour intensity reflects expression levels, with purple representing higher expression and orange representing lower expression. Dendrograms above the heatmap (samples) and to the left (genes) illustrate similarity in expression profiles, clustering individuals and genes based on shared expression patterns.
Consistent with these observations, the highest number of DEGs was identified in the AA versus GG genotype comparison (244 DEGs; Supplementary File 2), considerably exceeding those found in other genotype comparisons (135 DEGs in AA vs. AG, and 41 DEGs in AG vs. GG; Figure 2a). Intriguingly, the majority of DEGs were underexpressed in the farming‐associated GG genotype (218 underexpressed vs. 26 overexpressed; Figure 2b). The adjusted *p‐*values and log_2_ fold changes for the nine most statistically significant DEGs are summarised in Table 1. Notably, 167 DEGs identified in the AA versus GG comparison were unique to this contrast, representing 68.5% of the total DEGs identified across all genotype comparisons (Figure 2a). Among these unique DEGs, 154 were found downregulated in the farming‐associated GG genotype.
(a) Venn diagram illustrating the distribution of differentially expressed genes (DEGs) across the three pairwise comparisons between SNP vgll3 genotypes in gilthead seabream juveniles. Numbers represent DEGs unique to each comparison as well as those shared among comparisons, highlighting both specificity and overlap. (b) Volcano plot depicting DEGs identified in the comparison between AA and GG genotypes. Genes significantly under‐expressed in the GG genotype (p adj < 0.05) are shown in orange, while significantly over‐expressed genes (p adj < 0.05) are shown in purple. Genes without significant differential expression (p adj > 0.05) are indicated in grey.
TABLE 1: The nine most statistically significant differentially expressed genes (DEGs) identified in the AA versus GG gilthead seabream juveniles genotype. All log2 fold change values are negative, indicating lower expression of these genes in the farmed GG genotype compared to AA genotype. p adj is the adjusted p‐value (Benjamini–Hochberg correction for multiple testing), and log2 fold change is the logarithm (base 2) of the fold change in expression between AA and GG genotypes.
The enrichment analysis of the 244 DEGs identified in the AA versus GG comparison revealed 30 significant Gene Ontology (GO) terms, with eight terms considered most ‘overrepresented’ gene sets (Table 2). These terms predominantly indicate roles related to tissue structure and organisation, with notable examples including contractile fibre (GO:0043292; p adj = 3.572e‐6) and collagen trimer (GO:0005581; p adj = 2.538e‐5) emphasising structural and functional aspects of muscle tissues. Notably, pathway enrichment analysis also identified muscle contraction (R‐DRE‐397014; p adj = 2.137e‐3) as the sole significantly enriched pathway (Table 2). Analysis of the promoter sequences for motif enrichment revealed significantly enriched transcription factor binding sites (TFBSs) for several transcription factors, including Pax‐6, MEIS1B/HOXA9, DLX2 and FOX family members (Table S1), suggesting potential upstream regulators of the observed gene expression changes.
TABLE 2: The eight most ‘overrepresented’ Gene Ontology (GO) terms and molecular pathways associated with differential expression between AA and GG SNP vgll3 genotypes. Term size refers to the total number of genes linked to each GO term or pathway in the reference database. Number of identified genes indicates the subset of genes from each term that were identified as differentially expressed genes in the analysis. p adj is the adjusted p‐value (Benjamini–Hochberg correction for multiple testing).
Interactome‐based analysis of DEGs in the AA versus GG comparison further provided insights into functional interactions. Specifically, this analysis identified six distinct gene interaction networks. Among these networks, we focused exclusively on the largest network comprising 12 genes, as networks with fewer genes (2–4 genes each) offered less comprehensive biological interpretation. In this 12‐gene network (Figure 3a), all genes were found to be consistently under‐expressed in the farming‐associated GG genotype compared to the AA genotype (Figure 3b). The most significant and ‘overrepresented’ GO terms included calcium ion binding (GO:0005509; p adj = 4.814e‐4), regulation of muscle system process (GO:0090257; p adj = 1.120e‐6), and sarcomere (GO:0030017; p adj = 2.295e‐7), thus underscoring their roles in muscle physiology (Table 3). Additionally, pathway enrichment highlighted two relevant pathways, namely muscle contraction (R‐DRE‐397014; p adj = 2.040e‐6) and cardiac conduction (R‐DRE‐5576891; p adj = 4.920e‐5), further emphasising the muscle‐related biological impact of the identified DEGs (Table 3). Expanding this network using an enhanced interactome that incorporated experimentally validated database‐curated interactions added seven additional genes, and the resulting network remained similarly enriched for functional GO terms (Figure S1).
(a) Gene interaction network derived from differentially expressed genes (DEGs) in the AA vs. the farming‐associated GG genotype comparison. Nodes represent genes and edges interactions. Edge thickness corresponds to the likelihood of interaction, with thicker edges indicating higher confidence (score range: 0.407–0.732). Node colours distinguish functional categories: green indicates calcium regulation (cacng1b, casq1a, casq2, pvalb2, pvalb4), brown represents actin‐myosin dynamics (tnnc2, tnnt3b, tnnt2e, mybphb, tpma), and grey denotes muscle tissue integrity and signalling (nrap, myoz1a). (b) Heatmap illustrating relative expression levels of the genes included in the interaction network. Rows represent DEGs, while columns represent individual samples grouped by genotype. Colour intensity indicates expression levels, with purple representing higher expression and orange representing lower expression. Dendrograms above the heatmap (samples) and to the left (genes) show clustering based on expression profiles similarity.
TABLE 3: The six most significant and ‘overrepresented’ Gene Ontology (GO) terms and molecular pathways associated with the gene interaction network of differentially expressed genes between AA and GG SNP vgll3 genotypes. Term size refers to the total number of genes linked to each GO term or pathway in the reference database. Number of identified genes indicates the subset of genes from each term that were identified as differentially expressed genes in the analysis. p adj is the adjusted p‐value (Benjamini–Hochberg correction for multiple testing).
In the AG vs. GG genotype comparison, 41 DEGs were identified (Figure 2a). Enrichment analysis of these DEGs revealed significant GO terms predominantly associated with vision and perception processes, including structural constituent of eye lens (GO:0005212; p adj = 1.034e‐13), lens development in the camera‐type eye (GO:0002088; p adj = 3.220e‐10), and visual perception (GO:0007601; p adj = 2.004e‐8). These findings were largely driven by a substantial proportion (12 out of 41) of crystallin genes among the DEGs in this genotype comparison. In the AA versus AG comparison, 135 DEGs were identified (Figure 2a). However, the GO and pathway enrichment analysis did not reveal any statistically significant results.
Three DEGs were common across all genotype comparisons (Figure 2a): cacng1b (calcium channel, voltage‐dependent, gamma subunit 1b), igfn1.1 (immunoglobulin‐like and fibronectin type III domain containing 1, tandem duplicate 1), and lsp1a (lymphocyte‐specific protein 1a). Notably, in the comparison of homozygous genotypes (AA vs. GG), these genes also ranked among the most statistically significant DEGs (Table 1), with cacng1b further involved in the interaction network depicted in Figure 3a. These three genes showed the lowest expression levels in the GG genotype in the AA versus GG genotype comparison (Table 1).
Discussion
4
Our findings support a genotype‐dependent regulatory role for vgll3 in gilthead seabream at the studied developmental stage, highlighting its potential connection to rapid evolutionary changes under aquaculture conditions (Moulistanos et al. 2023). In Atlantic salmon, vgll3 acts as a master regulator, influencing thousands of genes related to sexual maturation and growth traits, such as body length and body condition (Ahi et al. 2022; Ahi, Verta, et al. 2025a; Debes et al. 2021; Kurko et al. 2020). In mice, vgll3 similarly influences key developmental and physiological processes, including muscle development (Figeac et al. 2019; Takakura et al. 2023). Consistent with this evolutionarily conserved role, our analysis revealed significant enrichment of muscle‐related biological processes among the DEGs identified in gilthead seabream (Tables 2 and 3). Collectively, these parallels suggest that vgll3 exerts comparable regulatory functions across vertebrates, including across divergent teleost lineages, although species‐specific life‐history traits (e.g., protandry in gilthead seabream, and anadromy in Atlantic salmon) may modulate how these mechanisms are manifested.
Specifically, we observed markedly divergent transcriptome profiles among SNP_ vgll3 _ genotypes previously associated with aquaculture‐driven selection in gilthead seabream (Moulistanos et al. 2023). Biological replicates of the farming‐associated GG genotype, which exhibits lower vgll3 expression (Moulistanos et al. 2025), consistently displayed transcriptomic profiles that were clearly distinct from those of the homozygous AA genotype (Figure 1). The transcriptomes of the heterozygous AG genotypes clustered closely with those of the GG genotypes, whereas their divergence from the AA genotypes was less pronounced than the contrast observed between the GG and AA genotypes (Figures 1 and 2). These findings align with previous qPCR analyses, which demonstrated similar vgll3 expression levels between AG and GG genotypes, but significantly higher in AA fish (Moulistanos et al. 2025). Collectively, these results support a dominance effect of the G allele on downstream gene regulation, potentially overriding the influence of the A allele on transcriptional profiles (Jiang et al. 2020).
Differential expression and interactome analyses comparing the homozygous genotypes (AA vs. GG) revealed an enrichment of genes associated with muscle physiology, supporting a role for vgll3 in muscle development and growth in gilthead seabream (Tables 2 and 3). Importantly, our study focused on 69 days post‐hatching (dph) juveniles, a critical developmental stage marking the transition from hyperplastic (fibre proliferation) to hypertrophic (fibre enlargement) muscle growth in this species (Rowlerson et al. 1995). Consistent with previous findings identifying vgll3 as a key modulator of skeletal muscle myogenesis (Figeac et al. 2019) and slow‐twitch muscle differentiation (Takakura et al. 2023), our interactome analysis identified multiple DEGs directly linked to muscle function (Figure 3). These included genes involved in calcium regulation (cacng1b, casq1a, casq2, pvalb2, pvalb4) (Chu et al. 2015; El Ghaleb et al. 2022; Freise et al. 2000; Infante et al. 2011), actin‐myosin dynamics (tnnc2, tnnt3b, tnnt2e, mybphb, tpma) (Ohtsuki et al. 2021), and muscle tissue integrity and signalling (nrap, myoz1a) (Baker et al. 2010; Casey et al. 2023; Luo et al. 2018; Zhou et al. 2022) (Figure 3). Taken together, these findings suggest that vgll3 is associated with muscle growth in gilthead seabream by modulating calcium homeostasis, contractile dynamics, and muscle structural organisation during a key developmental stage.
Our data also suggest a potential role for vgll3 in extracellular matrix (ECM) organisation, in line with previous evidence linking vgll3 to ECM‐related processes, such as collagen synthesis and tissue integrity (Haque et al. 2023; Horii et al. 2023; Kular et al. 2014). Given the essential role of ECM in muscle function and growth (Csapo et al. 2020), vgll3‐mediated regulation of ECM components further supports its involvement in calcium homeostasis and contractile machinery, together shaping muscle development. This regulatory axis has important ecological and aquaculture implications: for example, sustained swimming stimulates muscle contraction and increases white muscle mass and fibre density in gilthead seabream, thereby enhancing growth performance (Ibarz et al. 2011; Moya et al. 2019). These mechanisms may underlie fitness‐related variation and help explain the putative signatures of selection observed at vgll3 in farmed populations (Moulistanos et al. 2023).
Beyond muscle growth, vgll3 may also have a role in sexual maturation and sex determination in gilthead seabream, highlighting its broad regulatory significance. Among the nine most statistically significant DEGs between AA and GG genotypes (Table 1), amh, cancg1 and igfn1.1 have documented roles in puberty and sex determination in various teleost species. Specifically, amh represses male puberty and is crucial for Sertoli cell differentiation (Morais et al. 2017; Verta et al. 2020), with reduced expression in the farming‐associated GG genotype paralleling observations in Atlantic salmon, where early‐maturing vgll3 genotypes also showed reduced amh expression in juvenile testes (Verta et al. 2020). Additionally, cacng1b and igfn1.1 have established links to gonadal differentiation in other fish species as cacng1b is down‐regulated by miR‐202 specifically in male gonads of stickleback ( Gasterosteus aculeatus Linnaeus, 1758) (Kaitetzidou et al. 2022), while igfn1.1 was identified as a marker gene for genetic sex identification in Spotted knifejaw ( Oplegnathus punctatus Temminck & Schlegel, 1844), based on male‐specific insertion information (Ma et al. 2022). Although sex‐specific differences in vgll3 function have been studied extensively, such as allele‐dependent transcriptional differences between males and females in Atlantic salmon (Ahi, Verta, et al. 2025, 2025b), our study does not include sex comparisons given that all gilthead seabream juveniles are males due to the species' protandrous nature. Nevertheless, vgll3 exhibits a conserved role in regulating maturation and growth, highlighting its pleiotropic and evolutionarily conserved functions.
These functional observations also raise the question of how local genetic variation can propagate such widespread transcriptional effects. The SNP_ vgll3 _ variant previously linked to aquaculture selection in gilthead seabream is likely associated with altered vgll3 expression (Moulistanos et al. 2023, 2025); however, the pronounced downstream transcriptional divergence among genotypes suggests that this local regulatory effect propagates broadly through vgll3's role as a transcriptional cofactor. Similar scenarios, in which a local regulatory variant influences extensive gene networks via a key regulatory gene, have been described in other systems (Nica and Dermitzakis 2013; Ahi, Verta, et al. 2025a). For example, Small et al. (2011) demonstrated that variants can act as master trans regulators, modulating the expression of hundreds of downstream genes, while Liu et al. (2022) showed that eQTLs in plants can similarly control extensive gene networks through key regulatory loci. In our dataset, within this regulatory network, map4k4 was differentially expressed, highlighting it as an example among the DEGs and suggesting that vgll3 modulates Hippo signalling through upstream kinase regulation (Meng et al. 2015). In addition, transcription factors predicted to bind enriched motifs in DEGs, including SMAD3 and TCF7, are directly or indirectly linked to Hippo signalling and YAP/TAZ‐mediated transcription (Qin et al. 2018; Liu et al. 2024). In this context, vgll3 may function as a regulatory hub through which amplified effects are exerted on growth‐, muscle‐ and maturation‐related pathways.
Limitations and Future Directions
5
Our study provides novel insights into vgll3‐dependent regulatory pathways in gilthead seabream, but several considerations should be noted. First, the use of four biological replicates per genotype may limit statistical power and the generality of the findings, though clear genotype‐dependent transcriptomic patterns were still observed. Second, our analyses focused on 69 dph, a critical stage marking the transition from fibre proliferation to hypertrophy in gilthead seabream (Rowlerson et al. 1995), which has lasting impacts on growth trajectories and physiological adaptations. Evidence from Atlantic salmon has indicated that the vgll3 genotype influences gene expression in a differential fashion across developmental stages, with distinct transcriptional differences observed between mature and immature male salmon based on their vgll3 genotype (Ahi, Verta, et al. 2025a). Third, the absence of phenotypic data (e.g., growth or muscle metrics) for sequenced individuals further limits genotype–phenotype inferences.
In addition to these experimental design‐related constraints, whole‐body transcriptomic analyses are strongly influenced by dominant tissues such as muscle, which constitute a large proportion of juvenile body mass and therefore may disproportionately shape transcriptomic profiles (Shi et al. 2025). This dominance can mask signals from other tissues, including those involved in reproduction. Consequently, while our results suggest vgll3 may contribute to both muscle physiology and reproductive processes, they primarily reflect organism‐level transcriptional effects rather than tissue‐specific mechanisms. To resolve the spatial and temporal regulation of vgll3, future studies employing tissue‐specific or single‐cell transcriptomics across developmental stages will be essential.
Interpretation of these transcriptional patterns is further complicated by genetic and regulatory considerations related to gene expression. In particular, regulatory inference for vgll3 is challenged by the synonymous SNP (SNP_ vgll3 ) analysed in this study. Although synonymous variants do not alter protein sequences, they can affect gene regulation by influencing mRNA stability, translation efficiency, or codon usage (Ramírez‐Bello and Jiménez‐Morales 2017; Robert and Pelletier 2018; Mitra et al. 2016; Oelschlaeger 2024). Such effects could potentially modulate VGLL3 protein availability for interaction with TEAD transcription factors, thereby influencing muscle differentiation and growth pathways (Hori et al. 2020; Mesrouze et al. 2020; Mesrouze and Chène 2024; Joshi et al. 2017; Tsika et al. 2008). Importantly, while SNP vgll3 _ represents a plausible candidate underlying the observed GG‐associated downregulation of genes, the current data do not allow us to establish a direct causal relationship. The SNP may instead act as a linked marker within the vgll3 genomic region, with contributions from nearby regulatory variants or epistatic interactions. Disentangling these possibilities will require targeted functional assays to directly test causality.
At the same time, these genetic considerations intersect with functional annotation, which in this study relied on Danio rerio orthologs that provide relatively well‐curated Gene Ontology information among teleosts. Nevertheless, the enrichment of GO terms, pathways, and interaction networks should not be interpreted as direct mechanistic evidence linking vgll3 to somatic growth, reproductive maturation, or domestication‐related traits, but rather as statistically supported hypotheses that require experimental validation, as cautioned in studies of selection and functional inference (Pavlidis et al. 2012). The patterns revealed in enrichment analyses reflect meaningful biological trends in the gilthead seabream, and the conclusions of this study remain robust despite the inherent limitations of cross‐species functional inference.
These considerations are further contextualised by the broader evolutionary framework of the VGLL gene family. The gilthead seabream gene investigated, vgll3, belongs to the vestigial‐like (VGLL) gene family, which has undergone multiple duplications in vertebrates, giving rise to paralogues (vgll1‐4) with varying levels of annotation across species. While sequence similarity supports assignment of the seabream locus to the vgll3 clade, functional divergence among paralogues cannot be fully excluded. Accordingly, cross‐species functional comparisons of VGLL3 should be interpreted with appropriate caution. We used zebrafish orthologs for functional annotation due to their well‐curated Gene Ontology and pathway resources (Primmer et al. 2013). This approach is most informative for identifying broad, evolutionarily conserved biological processes, rather than precise species‐specific functions. Nevertheless, the genotype‐associated transcriptional signatures observed associations in gilthead seabream are broadly consistent with the reported of VGLL3 in growth, maturation, and developmental regulation in other vertebrates (Barson et al. 2015; Cousminer et al. 2013; Debes et al. 2021; Halperin et al. 2013; Pennonen 2017; Perry et al. 2014), supporting the hypothesis of partial functional conservation. Future work incorporating dedicated phylogenetic analyses or schematic representations of VGLL family evolution across teleosts would help clarify orthology relationships and strengthen evolutionary and functional interpretation.
Future studies incorporating larger sample sizes, multi‐stage, tissue‐specific, and functional analyses, along with phenotypic measurements, will be essential to fully elucidate vgll3's regulatory roles and its contribution to ecologically relevant traits. Specifically, allele‐specific expression analyses, CRISPR/Cas9‐mediated functional validation, and proteomic profiling could clarify the impact of candidate variants such as SNP_ vgll3 _ on protein function and downstream pathways. Additionally, longitudinal studies across developmental stages and environmental conditions would help reveal how *vgll3‐*mediated regulatory networks respond to ecological pressures. While exploratory, these findings provide a foundation for understanding the evolutionary and functional significance of vgll3 in gilthead seabream and establish a framework for future research on its role in muscle growth, reproductive development, and potential applications in aquaculture breeding programs.
In conclusion, our study demonstrated a genotype‐dependent role for vgll3 in regulating transcriptomic divergence during early development in gilthead seabream. The data provide strong evidence for genotype‐dependent differential expression associated with regulatory variation at the vgll3 locus, with downstream transcriptional patterns consistent with vgll3 acting as a central regulator of multiple genes. Distinct transcriptomic profiles observed among SNP_ vgll3 _ genotypes, combined with previous qPCR expression analyses, indicated vgll3 as a significant regulator of developmental and physiological pathways linked to putative domestication‐related traits. Differential expression analyses highlighted enriched muscle‐related pathways and suggested that vgll3 also appears to play an important role in reproduction processes, reflecting its pleiotropic and evolutionarily conserved functions. While our findings at 69 days post‐hatch provided valuable insights, they also underscore the need for tissue‐specific, multi‐stage studies incorporating phenotypic measurements to fully elucidate vgll3's comprehensive regulatory role. Collectively, vgll3 emerges as a promising candidate for further investigation in relation to muscle development, growth, and sexual and reproductive maturation in gilthead seabream, although additional studies are required to clarify its phenotypic impact and potential utility in selective breeding.
Author Contributions
Aristotelis Moulistanos: data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); visualisation (lead); writing – original draft (equal). Elisavet Kaitetzidou: data curation (equal); formal analysis (equal); methodology (equal); validation (equal); writing – review and editing (equal). Styliani Minoudi: data curation (equal); formal analysis (equal); methodology (equal); writing – review and editing (equal). Konstantinos Gkagkavouzis: methodology (equal); project administration (equal); writing – review and editing (equal). Efthimia Antonopoulou: investigation (equal); methodology (equal); supervision (equal); writing – review and editing (equal). Alexandros Triantafyllidis: conceptualisation (equal); funding acquisition (equal); investigation (equal); methodology (equal); project administration (equal); supervision (equal); writing – review and editing (equal). Spiros Papakostas: conceptualisation (equal); data curation (lead); formal analysis (lead); funding acquisition (equal); investigation (equal); methodology (equal); project administration (equal); resources (lead); supervision (equal); validation (equal); visualisation (equal); writing – original draft (equal); writing – review and editing (lead).
Funding
This study was conducted under the project ‘SEaLIFT: SystEms Biology Modelling of Key LIFe History Traits for Sustainable Aquaculture Production in the Mediterranean Region’, funded by the Hellenic Foundation for Research & Innovation (H.F.R.I.) Grant Number 00414.
Ethics Statement
The ‘Principles & Procedures Governing the Operation of the Research Ethics Committee’ of the University in which the research was conducted (Aristotle University of Thessaloniki). This regulation is drafted by the provisions of Article 68, Law 4485/2017, and Articles 21–27, Law 4521/2018 of Greece, an EU member country. Article 25 of the regulation (‘Research Involving the Use of Animals’) states that the implementation of our experimental protocols did not use endangered species or wildlife animals (i.e., aquaculture‐reared fish) and is part of routine practice in aquaculture practice (Article 25, section 5.d); such procedures do not require approval.
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting information
Supplementary File 1. Raw gene expression count matrix obtained from RNA‐seq analysis of gilthead seabream juveniles, used for downstream differential gene expression analyses.
Supplementary File 2. Genomic coordinates and differential gene expression results of gilthead seabream transcripts identified from RNA‐seq analysis, including chromosome position, gene annotation, strand information and DESeq2 statistical outputs (p‐value, FDR, and fold‐change estimates).
Figure S1: Complementary gene interaction network derived from differentially expressed genes (DEGs) in the AA versus farming‐associated GG genotype comparison. The network was constructed using STRING evidence based on experimentally validated interactions, curated databases, and co‐expression. Nodes represent genes and edges represent interactions, with edge thickness indicating interaction confidence. Node colours are consistent with Figure 3a, with pink indicating seven additional genes incorporated in the enhanced network (smarce1, acta1a, tpm2, myl2b, tmod4, mylz3, rasgrp4). The network showed highly significant functional enrichment for muscle‐related processes, including sarcomere (GO:0030017, p adj = 5.3 × 10^−10^), muscle system process (GO:0003012, p adj = 1.1 × 10^−7^), myofibril assembly (GO:0030239, p adj = 2.6 × 10^−6^), and regulation of muscle contraction (GO:0006937, p adj = 2.6 × 10^−6^), as well as calcium ion binding (GO:0005509, p adj = 1.7 × 10^−4^), supporting a central role for contractile and excitation–contraction coupling mechanisms in VGLL3‐associated genotype differences.
Table S1: Transcription factor binding site (TFBS) enrichment analysis of promoter regions of DEGs in the AA vs. GG comparison. Motifs significantly enriched in DEG promoters were identified using g:Profiler with the TRANSFAC database. The table lists the transcription factor name, consensus motif sequence, Term ID, and adjusted p‐value (padj) for each TFBS.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Ahi, E. P. , B. Panda , and C. R. Primmer . 2025. “The Hippo Pathway: A Molecular Bridge Between Environmental Cues and Pace of Life.” BMC Ecology and Evolution 25, no. 1: 35. 10.1186/s 12862-025-02378-8.40275190 PMC 12020181 · doi ↗ · pubmed ↗
- 2Ahi, E. P. , M. Sinclair‐Waters , J. Moustakas‐Verho , S. Jansouz , and C. R. Primmer . 2022. “Strong Regulatory Effects of vgll 3 Genotype on Reproductive Axis Gene Expression in Juvenile Male Atlantic Salmon.” General and Comparative Endocrinology 325: 114055. 10.1016/j.ygcen.2022.114055.35580687 · doi ↗ · pubmed ↗
- 3Ahi, E. P. , J. P. Verta , J. Kurko , A. Ruokolainen , P. V. Debes , and C. R. Primmer . 2025. “Hippo‐vgll 3 Signaling May Contribute to Sex Differences in Atlantic Salmon Maturation Age via Contrasting Adipose Dynamics.” Biology of Sex Differences 16: 23. 10.1186/s 13293-025-00705-8.40176157 PMC 11966934 · doi ↗ · pubmed ↗
- 4Ahi, E. P. , J.‐P. Verta , J. Kurko , et al. 2025 a. “Gene Co‐Expression Patterns in Atlantic Salmon Adipose Tissue Provide a Molecular Link Among Seasonal Changes, Energy Balance and Age at Maturity.” Molecular Ecology 34, no. 15: e 17313. 10.1111/mec.173131-16.38429895 PMC 12288797 · doi ↗ · pubmed ↗
- 5Ahi, E. P. , J.‐P. Verta , J. Kurko , et al. 2025 b. “Brain‐Associated Alterations of Hippo Pathway Transcription in Early Maturing Atlantic Salmon.” BMC Ecology and Evolution 25, no. 1: 53. 10.1186/s 12862-025-02398-4.40419966 PMC 12105140 · doi ↗ · pubmed ↗
- 6Albert, F. W. , and L. Kruglyak . 2015. “The Role of Regulatory Variation in Complex Traits and Disease.” Nature Reviews Genetics 16, no. 4: 197–212. 10.1038/nrg 3891.25707927 · doi ↗ · pubmed ↗
- 7Baker, J. , G. Riley , M. R. Romero , et al. 2010. “Identification of a Z‐Band Associated Protein Complex Involving KY, FLNC and IGFN 1.” Experimental Cell Research 316, no. 11: 1856–1870. 10.1016/j.yexcr.2010.02.027.20206623 · doi ↗ · pubmed ↗
- 8Barson, N. J. , T. Aykanat , K. Hindar , et al. 2015. “Sex‐Dependent Dominance at a Single Locus Maintains Variation in Age at Maturity in Salmon.” Nature 528, no. 7582: 405–408. 10.1038/nature 16062.26536110 · doi ↗ · pubmed ↗