Integrated Analysis of Skin Development and Pigmentation in Amphiprion frenatus: Insights from Histology, Transcriptomics, and Translatomics

1. Introduction

Teleosts exhibit a remarkable diversity of body color patterns that serve various functions, including camouflage, communication, and mate attraction [1]. The formation of these color patterns involves a complex interplay of genetic, developmental, and environmental factors. Histologically, the skin of teleosts comprises three main layers: the epidermis, the dermis, and the subcutaneous layer. The epidermis is the outermost layer and primarily comprises stratified squamous epithelial cells, which provide a protective barrier against the external environment. The skin of teleost fish is closely linked to pigment cells known as chromatophores. They can change shape and size, allowing them to adjust the intensity and distribution of pigments in response to various stimuli such as light, stress, and social interactions.

Chromatophores are derived from neural crest cells and consist of six subtypes [2]. The main types include (1). Melanophore: These contain melanin, which is responsible for black and brown colors. Melanophores can expand or contract, allowing the fish to darken or lighten its skin. (2) Xanthophores: These chromatophores contain carotenoids and produce yellow colors. They often interact with other chromatophores to create a variety of shades. (3) Erythrophores: Like xanthophores, erythrophores contain carotenoids but primarily produce red hues. They can influence overall coloration by layering with other chromatophores. (4) Cyanophores: These contain blue pigments, although they are less common. They contribute to the blue coloration found in some species. (5) Iridophores: These chromatophores contain guanine crystals that reflect light, creating iridescent effects. The colors they produce can vary depending on the angle of light, enhancing the overall vibrancy of a fish’s coloration. (6) Leucophores: These reflect light and provide a white or light coloration, often enhancing the brightness of other pigments. The diversity of coloration in fish is controlled by many genes related to pigment synthesis, transport, and cell morphology [3].

Numerous studies have identified how genetic factors influence coloration in specific species. Zebrafish research has provided valuable insights into pigment cell development, and more recent multi-omics studies have begun to characterize pigment pathways at both the transcriptomic and translational levels. In cichlid fishes, genetic studies demonstrate how opsin genes contribute to variations in light perception and body coloration, which facilitates environmental adaptation. Teleost body coloration is the result of complex genetic interactions between pigment-producing cells (chromatophores) and regulatory genes. These genetic factors are finely tuned through evolutionary processes to help teleosts adapt to their environments, attract mates, or deter predators. Genes such as MC1R [4], TYR [5], ASIP [6], etc. play fundamental roles in determining specific colors and patterns in different teleost species, and their variability within populations leads to striking diversity in fish coloration.

In clownfish, coloration plays a critical role in their ecology and behavior [7]. The molecular basis of color variation in Amphiprion species is still under investigation [8,9]. Recently, de novo transcriptome sequencing has been used to identify genes involved in pigment synthesis in clownfish [10]. Several recent studies have used transcriptomic and genomic approaches to investigate the molecular basis of coloration in anemonefishes. For example, a chromosome-level genome assembly of Amphiprion clarkii revealed duplications in pigmentation-related genes [8]. At the same time, transcriptome analyses in Amphiprion ocellaris linked RNA methylation changes to variations in skin color [11]. Additionally, transcriptomic datasets in Amphiprion percula have served as valuable resources for exploring the expression of genes related to coloration [10].

In this study, we aimed to investigate the molecular and cellular basis of fish body color by combining fish body color ontology, histology, gene expression analysis, and transcriptomic and translatomic approaches. Using histological techniques, we examined the underlying cellular structures responsible for coloration. To further understand the genetic regulation of these color traits, we analyzed gene expression profiles and identified key genes involved in pigmentation. Finally, we examined the transcriptome and translatome to determine the molecular networks and regulatory pathways that govern color patterns at the genetic and protein levels. Analyzing both transcriptomes and translatomes is crucial for comprehending the comprehensive scope of gene regulation, particularly in systems as dynamic and responsive as teleost skin. Comparing RNA expressions with protein synthesis can reveal the precise molecular mechanisms that control complex traits, such as body coloration, and show how gene regulation is fine-tuned through transcriptional and translational processes. This integrated approach will provide valuable insights into the genetic, cellular, and molecular mechanisms driving fish coloration and its evolutionary significance.

2. Materials and Methods

2.1. Experimental Fish

Mature tomato clownfish (A. frenatus) were paired in each 50 L aquarium with several radiation sea anemones. The temperature and salinity of tank water were 27±1 ℃ and 34±1 psu, respectively. The fish were acclimated to a photoperiod with 9:15 h light: dark. The fish were fed fresh shrimp, squid, and a commercial feed (GROBEST GROUP LIMITED) ad libitum. After spawning, the newly hatched larvae were fed rotifers and, latterly, copepods. The EFRC Institutional Animal Care and Use Committee approved all experiments and procedures performed by the guiding ethical principles.

2.2. Skin Histology for Light Microscopy

A. frenatus were euthanized with MS222 in seawater and dissected to collect skin samples from the dorsal region [12]. The sampled tissues were fixed in Bouin’s solution and embedded in paraffin. The samples were sectioned at 5 µm, stained with hematoxylin and eosin (H&E) for general morphology, and observed under a light microscope.

2.3. Scanning Electron Microscopy

Larva A. frenatus, each about 10 mm in total length, were euthanized with MS222 in seawater and cut into three parts, head, middle trunk, and caudal parts, which were about 3 mm in length each, and were preserved in a fixative solution (0.2 M phosphate buffer solution (PBS) at pH 7.2, 1% glutaraldehyde and 2.5% paraformaldehyde) at 4℃. Targeted portions were preserved in a fixative solution for an additional 24 hr. After washing with 0.1 M PBS, they were post-fixed with a chilled 1% osmium tetroxide (OsO4) solution for 2 hours, followed by a final wash in 0.1 mol/L PBS. The fixed portions were dehydrated in an ethanol series and then placed into a critical-point dryer (k850, Quorum, United Kingdom) for a final 15 min of critical-point drying. The dried samples were mounted onto aluminum stubs and coated with a gold layer using an ion sputter (SC7620, Quorum, United Kingdom). Finally, samples were observed with a Hitachi SU1510 scanning electron microscope (Tokyo, Japan).

2.4. RNA Extraction and cDNA Synthesis

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions [10]. First-strand cDNA was synthesized using Superscript III First-Strand Synthesis Supermix (Invitrogen) with oligo (dT)20 primers. This cDNA (1 μg/40 μL Milli-Q water) was used for RT-PCR and quantitative RT-PCR (qRT-PCR) analyses.

2.5. Sequencing and Analyzing Pipeline

The Illumina HiSeq platform is widely used for RNA sequencing (RNA-Seq) because of its high throughput and accuracy. RNA is extracted using commercial kits, ensuring it is free of DNA, proteins, and other contaminants. Integrity is assessed using the Bioanalyzer RNA integrity number (RIN) or an equivalent score, which measures the quality. A RIN ≥ 7 is preferred for most RNA-Seq experiments. RNA is converted into a library of cDNA fragments for sequencing. The steps include: RNA selection/enrichment, fragmentation, reverse transcription, adapter ligation, and PCR amplification. Electrophoresis-based methods are used to verify the quality and fragment size distribution of the final cDNA library. Library concentration is quantified using fluorometric methods. The library is loaded onto an Illumina flow cell. Each DNA fragment hybridizes to complementary oligonucleotides immobilized on the flow cell. Each DNA strand is sequenced base by base using fluorescently labeled nucleotides. The sequencer captures fluorescent signals as nucleotides are incorporated. This process is repeated for each cycle to build the sequence read. The sequencing run generates millions of short reads (typically 50-150 bp each). The output is in FASTQ format, containing sequence data and corresponding quality scores. Reads are aligned to a reference genome or transcriptome using STAR. The number of reads mapping to genes or transcripts is quantified using the software feature Counts. Statistical analysis is performed to identify differentially expressed genes or transcripts using the edgeR package.

2.6. Ribosome Profiling Library Preparation Protocol

This protocol [16] has been adapted with moderate modifications. Briefly, skin tissues were lysed with polysome buffer [20 mM Tris-Cl [pH7.4], 150 mM NaCl, 5 mM MgCl2, 1 mM DTT, 200µg/ml cyclohexamide, 1 % (vol/vol) Triton X-100 and 25U/ml Turbo DNase I (Ambion, Cat# AM2238)] and lysate were digested with RNAse I (Ambion, Cat# AM2294). For RNAseq library construction, total RNA was extracted and oligo-dT magnetic beads were used to select RNA containing poly-A tails. Selected RNA was fragmented, and 50-60 nucleotide-long fragments were selected following end-repairing, linker ligation, reverse transcription, and circular PCR. Ribosome-protected fragments (RPF) purified with PAGE/fragmented total RNA were dephosphorylated with T4-PNK (NEB, Cat# M0201S). Linkers were ligated to RPF/fragmented total RNA by T4 Rnl2tr K227Q ligase (NEB, Cat# M0351L), and the cDNA synthesis was performed with SuperScriptIII (Invitrogen, Cat# 18080- 093). The cDNA samples were circularized by CircLigase (Epicentre, Cat# CL4111). The DNA libraries were amplified with Phusion High-Fidelity DNA Polymerase (NEB, Cat# M0530S) and purified using the AMPure XP system (Beckman Coulter, Cat# A63880). The quantified libraries were sequenced with the Illumina sequencing platform following the manufacturer’s instructions.

2.7. Bioinformatics Analysis

All analytic pipelines were developed under Linux CentOS 6.6 environment. Due to A frenatus being a non-model species, it is challenging to study gene expression and differential expression without a reference genome or complete gene annotation. In this study we de novo assembled the A frenatus transcriptome using Trinity based on previous RNA-seq data [12]. The transcriptome is the reference for mapping reads and quantifying gene expression. We then annotated all transcripts using Trinotate. The linker sequences were trimmed from each sequence read by using Cutadapt and selected fragments between 27-40 nucleotides (nts) for Ribo-seq (25~65 nts for RNA-seq) with in-house developed scripts. We mapped the remains to A. frenatus transcriptome reference from Trinity assembly by RSEM [16]. We calculated the expression value and count for each gene at both the gene and isoform levels. For differential expression analysis, we used edgeR [17], and read counts as input to identify differentially expressed genes between groups with replicates. The log2 fold change of RNA-seq and Ribo-seq results was plotted to illustrate the transcription and translation status of each gene.

2.8. qRT-PCR

cDNA derived from 10 ng RNA was used for Quantitative PCR performed by the following conditions: 10 min at 95°C, and 40 cycles of 15 sec at 95°C, 1 min at 60°C using 2X Power SYBR Green PCR Master Mix (Applied Biosystems) and 200 nM of forward and reverse primers in 20 μL. Each assay was run in triplicate on an Applied Biosystems 7900HT Real-Time PCR system, and expression fold-changes were derived using the comparative CT method, with ef1a as the endogenous control and sample “AF-25D” as the calibrator. Table 1 lists the forward and reverse oligonucleotide primer sequences for quantitative real-time PCR (qPCR) gene expression analysis. The primer sets were used to amplify transcripts from various genes involved in cellular development, differentiation, signaling, and pigmentation in fish skin tissues. Primers were designed using Primer Express software and synthesized by Mission Biotech. Gene symbols correspond to those analyzed in Figure 5a and Figure 5b. Expression levels were normalized to ef1α (not listed in this table as the internal reference gene, and relative expression was calculated using the 2^^−ΔΔCt method [18].

3. Results

3.1. Morphological Development and Pigmentation Progression in A. frenatus

The morphological development and pigmentation progression of A frenatus from the larval to juvenile stages were documented through imaging (Figure 1). Newly hatched larvae (0 day post-hatch; dph) exhibited a wholly transparent body with minimal structural differentiation and no apparent pigmentation (Figure 1A). By 1 dph, yolk absorption had begun. Still, transparency persisted (Figure 1B). At 5 dph, the earliest signs of melanophore accumulation appeared, primarily along the dorsal and caudal regions (Figure 1C), which intensified by 10 dph, giving the body a faint grayish tone (Figure 1D). By 12 dph, the body had darkened further. The fin structures became more distinct (Figure 1E). At 15 dph, a brownish hue emerged along the trunk, and eye pigmentation increased (Figure 1F). The appearance of an orange coloration was first evident at 20 dph, beginning in the posterior body (Figure 1G), and became more pronounced throughout the fins and body by 25 dph (Figure 1H). At 40 dph, the first white vertical head bar—a hallmark of the species—was observed (Figure 1I). As development continued, the orange pigmentation became dominant, and the white bar was well defined by 90 dph (Figure 1J). By 140 dph, individuals exhibited a vibrant reddish-orange body with stabilized bar patterns (Figure 1K). By 240 dph, juveniles displayed the fully mature coloration typical of the species, characterized by deep orange pigmentation and a distinct white head bar (Figure 1L). These observations indicate a sequential and coordinated pattern of pigment cell differentiation and regional pigmentation during early development.

3.2. Histological Development of the Skin in A. frenatus

Histological examination of A frenatus skin at various developmental stages revealed progressive structural and pigmentary maturation from 15 to 240 days post-hatch (dph) (Figure 2). At 15 dph, the epidermis consisted of a thin, undifferentiated epithelial sheet, with scattered melanophores present primarily in the epidermal layer, and the dermis remaining rudimentary (Figure 2A). By 25 dph, epidermal stratification became apparent, with keratinocytes beginning to form organized layers. Melanophores increased in number and migrated toward the dermo-epidermal junction. At the same time, xanthophores and early iridophores appeared in the superficial dermis (Figure 2B). At 40 dph, the epidermis displayed a more defined basal layer adjacent to the dermis, and chromatophore density increased, with all three pigment cell types—melanophores, xanthophores, and iridophores—discernible (Figure 2C). By 90 dph, the epidermis had become well-stratified. Pigment cells were abundant and spatially distributed across the epidermis and upper dermis. Iridophores and xanthophores were more prominent, indicating increased chromatophore diversification (Figure 2D). At 140 dph, the epidermis showed full stratification and a thickened basal layer, and dermal pigment cells, particularly melanophores, appeared elongated and interspersed with reflective iridophores and yellow-hued xanthophores, suggesting the onset of organized pigmentation patterns (Figure 2E). At 240 dph, skin architecture was fully mature. The epidermis was composed of multiple keratinocyte layers anchored by basal cells at the dermal interface, while the dermis exhibited thick collagen bundles and densely packed pigment cells. Chromatophores were well differentiated and stably localized, reflecting the establishment of mature skin pigmentation (Figure 2F). These findings illustrate the coordinated development of skin architecture and chromatophore populations during juvenile growth in A. frenatus.

3.3. Scanning Electron Microscopy (SEM) of Skin Structure

Scanning electron microscopy (SEM) provided detailed visualization of the surface architecture and histological organization of the skin in juvenile A frenatus. The skin surface was characterized by overlapping scales bearing posteriorly oriented ctenii, partially covered by the epithelial layer (EPL). In some regions, the EPL was removed during sample preparation, revealing underlying annuli (ans) (Figure 3A). Damaged areas of the epidermis (EPD) exposed epithelial cells (epc) and dense arrays of microridges (MRs), with annular structures clearly visible (Figure 3B). A lateral view of the trunk skin showed the multi-layered arrangement of the integument, with EPD and EPL overlaying the basement membrane (BM) and underlying dermis (DM), which contained fibroblasts (fb) (Figure 3C). The apical surface of intact EPD exhibited pronounced MRs, and the apical openings of goblet cells (ogbc) were discernible between the ridges (Figure 3D), suggesting active mucus secretion. Stratified epc and capillaries (cp) derived from dermal vasculature were also observed in association with MRs in the upper EPD (Figure 3E). At the basal end, epithelial cells were closely adhered to the BM, indicating structural cohesion within the epidermal layers (Figure 3F). These morphological features reflect the specialized skin organization in A. frenatus, supporting its protective barrier function and potential roles in mucosal immunity.

3.4. Transcriptomic Shifts Between 0 and 14 dph

RNA-Seq analysis of fish larvae is a practical approach to identify and predict genes regulated in developmental processes of 0 dph compared with 14 dph in larvae of A frenatus. We sequenced 4 cDNA libraries, each pooled from 48 samples of larvae hatched at 0 dph and 14 dph using the Illumina HiSeq platform. The sequencing runs yielded 352.90 M reads with an average length of 80 bp. The Strand NGS software (version 2.1) was used with default parameters for pre-alignment and post-alignment quality control analysis. Notably, 100% of the raw reads remained in the dataset when using the Danio rerio (zebrafish) genome for contig alignments. Of these, 68.4 M reads (19.4%) aligned into contigs, whereas 284.5 M reads (80.6%) remained unaligned singletons because of the genetic distance between the two species. The unaligned reads were excluded from further analytical steps. The high-throughput sequencing performed for different libraries yielded similar numbers of reads, ranging from 75.21 to 95.96 million, and had the same average length. The number of reads for each condition of the larvae sample and the number of reads mapped to reference sequences are shown in Table 2.

Differentially expressed genes (DEGs) between 0 and 14 dph were calculated from the raw reads using the DESeq method and moderated t-test analysis. The significance level applied to the RNA-Seq data was corrected for multiple tests using Benjamin-Hochberg correction, allowing rigorous identification of DEGs in the larvae of 0 and 14 dah. We identified 4245 DEGs out of a total of 31967 genes, which were hierarchically clustered into two major and well-separated groups (Figure 4A). We analyzed overall gene functions and uncovered specific biological processes and terms among 4245 DEGs by using the GO enrichment analysis and found that a majority of DEGs were significantly associated with GO metabolic process (1761 genes, p = 3.47E-33) and developmental process (378 genes, p = 1.40E-16) within the biological process category (2938 genes, p = 1.45E-25). The genes significantly associated with metabolic process term were mainly involved in catalytic activity (1368 genes, p = 4.94E-41), biosynthetic process (440 genes, p = 8.26E-20), and hydrolase activity (673 genes, p = 1.08E-10) (Figure 4B). However, the genes associated with the cellular metabolic process and involved in energy derivation (4 genes, p = 0.0976) and the RNA biosynthetic process (128 genes, p = 0.0833) were not significantly enriched. Interestingly, all 27 DEGs involved in glycolysis were enriched considerably and downregulated at 25 days after hatching. To identify specific altered functions, we investigated significantly associated biological processes in the subsets of DEGs up-regulated (2328 genes) and down-regulated (1917 genes) in the d25 group. Among the up-regulated genes, key enriched processes included oxidoreductase activity (207 genes, p = 2.40E-44), proteolysis (121 genes, p = 6.19E-10), and intermediate filament formation (22 genes, p = 2.53E-10), highlighting an increased emphasis on metabolic and structural regulation. Additionally, the myosin complex (27 genes, p = 2.48E-09) and pigment metabolic processes (10 genes, p = 2.01E-04) were notably upregulated (Figure 4C, left panel).

Conversely, the down-regulated subset was significantly enriched in developmental processes (258 genes, p = 1.65E-26), consistent with suppressing differentiation-related functions. Specific processes included cell differentiation (93 genes, p = 4.30E-11), neurogenesis (51 genes, p = 8.04E-09), and ameboidal cell migration (19 genes, p = 2.27E-04), indicating a broad downregulation of growth and morphogenesis pathways (Figure 4C, right panel). Furthermore, glycolysis (27 genes, p = 1.49E-22) was significantly reduced, suggesting a metabolic shift during this developmental stage. The enrichment of eye development (31 genes, p = 5.73E-03), visual perception (15 genes, p = 6.13E-06) and detection of light stimuli (9 genes, p = 1.26E-05) in the down-regulated subset further suggests a suppression of sensory-related genes. These findings reinforce the distinct regulation of metabolic activation and structural remodeling in up-regulated genes, contrasted with the repression of differentiation and developmental pathways in down-regulated genes, underscoring critical functional shifts at this stage.

3.5. Validation of Pigmentation-Related Gene Expression by Quantitative Real-Time PCR

Quantitative real-time PCR (qPCR) was performed on candidate genes in skin tissues of A. frenatus to validate transcriptomic findings and assess stage-specific gene expression during pigmentation development at 25 and 40 dph. Gene selection was based on their putative roles in cell proliferation, differentiation, pigment synthesis, and ion regulation.

As shown in Figure 5A, expression levels of genes associated with cellular proliferation and differentiation were significantly elevated at 40 dph compared to 25 dph. These included myc, pcna, and p63, proliferation and epithelial renewal markers, and pax6 and gata3, which are linked to developmental and pigment lineage commitment. Cytoskeletal genes such as myo5a, myo7a, and myh7b also showed marked upregulation, suggesting increased intracellular transport and structural remodeling activity. The pigmentation-associated gene ostm1 displayed a substantial increase at 40 dph, aligning with active melanogenesis. In a separate cluster of genes (Figure 5B), significant upregulation was also observed in regulators of melanosome transport (rab38), and pigment patterning (qnr71, sufu, sox2). Notably, mchr1 exhibited the highest fold changes, implicating them in ion homeostasis and pigment cell signaling. The transient expression of timp3, which was nearly undetectable at 25 dph but robustly induced by 40 dph, suggests dynamic regulation of extracellular structure during chromatophore stabilization. These expression patterns were consistent with RNA-Seq results and support the hypothesis that pigmentation onset involves coordinated upregulation of genes governing proliferation, differentiation, vesicle transport, and extracellular matrix reorganization. Collectively, the qPCR results confirm the developmental timing of gene activation underlying the acquisition of skin coloration in juvenile A. frenatus.

3.6. Integrated Transcriptional and Translational Regulation During Pigmentation

To evaluate the coordination between transcriptional and translational regulation during early pigmentation in A. frenatus, we performed an integrated analysis of RNA-Seq and Ribo-Seq data at 25 and 40 dph. Figure 6 presents a scatter plot comparing log₂ fold changes in mRNA levels (x-axis) and ribosome-protected fragment (RPF) levels (y-axis), representing protein synthesis. Each gene was classified into one of nine regulatory categories based on the direction and magnitude of its transcriptional and translational changes (cutoff |log₂FC| > 1). A substantial proportion of genes fell into Group C, which comprises genes concordantly upregulated at both the mRNA and protein levels. This group includes myc, pcna, p63, myo5a, myo7a, rab38, and mchr1, all of which were also validated by quantitative real-time PCR (qPCR) as upregulated at 40 dph. The co-induction of these genes at both transcriptional and translational levels underscores a robust, unified gene expression program during chromatophore differentiation and pigment pattern formation. These genes span multiple functional categories: myc and pcna are markers of proliferative activity; myo5a, myo7a, and rab38 participate in melanosome transport; mchr1 is involved in pigment cell signaling and ion homeostasis. Their coordinated activation suggests a tightly coupled regulatory system that ensures both transcript abundance and protein synthesis are synchronized during this critical developmental window.

Figure 5. Quantitative real-time PCR (qPCR) validation of gene expression in skin tissues of A. frenatus at different post-hatch stages. (A) Relative expression of genes associated with development, proliferation, and stemness, including myc, myo5a, myo7a, ostm1, p63, pax6, and pcna, in skin samples collected at 25 dph, and 40 dph. (B) Expression of genes related to pigmentation, ion transport, and structural regulation, including qnr71, rab38, sufu, tm154, trpm7, and mchr1. Gene expression levels were normalized to ef1a as the internal control and are presented as fold change relative to AF-25D using the 2^{^−ΔΔCt} method. Data are shown as mean ± standard error (SE) based on three technical replicates. The observed expression trends are consistent with transcriptome profiling results and support stage-specific regulation of genes involved in skin development and pigmentation.

Figure 5. Quantitative real-time PCR (qPCR) validation of gene expression in skin tissues of A. frenatus at different post-hatch stages. (A) Relative expression of genes associated with development, proliferation, and stemness, including myc, myo5a, myo7a, ostm1, p63, pax6, and pcna, in skin samples collected at 25 dph, and 40 dph. (B) Expression of genes related to pigmentation, ion transport, and structural regulation, including qnr71, rab38, sufu, tm154, trpm7, and mchr1. Gene expression levels were normalized to ef1a as the internal control and are presented as fold change relative to AF-25D using the 2^{^−ΔΔCt} method. Data are shown as mean ± standard error (SE) based on three technical replicates. The observed expression trends are consistent with transcriptome profiling results and support stage-specific regulation of genes involved in skin development and pigmentation.

Figure 6. Integrated analysis of transcriptional and translational gene expression during skin pigmentation (25 vs. 40 dph) in A. frenatus. Scatter plots show log₂ fold changes in mRNA (x-axis) and protein (y-axis) levels. Genes are grouped into nine regulatory categories based on expression patterns, including concordant (Groups C, G) and discordant regulation (e.g., Groups A, I). Dashed lines indicate significance thresholds (log₂ fold change > 1). The analysis highlights widespread post-transcriptional regulation during skin coloration.

Figure 6. Integrated analysis of transcriptional and translational gene expression during skin pigmentation (25 vs. 40 dph) in A. frenatus. Scatter plots show log₂ fold changes in mRNA (x-axis) and protein (y-axis) levels. Genes are grouped into nine regulatory categories based on expression patterns, including concordant (Groups C, G) and discordant regulation (e.g., Groups A, I). Dashed lines indicate significance thresholds (log₂ fold change > 1). The analysis highlights widespread post-transcriptional regulation during skin coloration.

While other regulatory categories were represented across the transcriptome, the pigmentation-related genes primarily exhibited concordant regulation, affirming that the molecular control of early pigmentation in A. frenatus is driven by integrated transcriptional and translational upregulation. This finding also validates the importance of multi-omics approaches in accurately capturing the dynamics of gene expression. Overall, the integrative analysis revealed that pigmentation in A. frenatus is governed by a multi-tiered regulatory landscape, where post-transcriptional and translational mechanisms play dominant roles alongside classical transcriptional control. Numerous genes in discordant groups highlight the limitations of using transcript-level data alone to infer functional outcomes. It highlights the necessity of parallel proteomic profiling to fully understand the orchestration of developmental processes such as skin coloration.

4. Discussion

The present study comprehensively characterizes the histological, transcriptional, and translational processes underlying skin pigmentation development in A frenatus. By integrating morphological observations, transcriptome and translatome profiling, and quantitative gene expression validation, we reveal a multi-layered regulatory framework that orchestrates the emergence of this species’ characteristic pigmentation patterns.

4.1. Morphological and Histological Progression of Skin Development

The histological analysis demonstrated that skin pigmentation in A. frenatus follows a defined temporal sequence. Initially, larvae exhibit a transparent body lacking visible pigment cells. By 5 dph, melanophore precursors emerge along the dorsal and caudal regions. As development progresses, the appearance of the orange pigmentation and the distinct white head bar at 40 dph mark pivotal milestones in chromatophore organization and pigment patterning. Histologically, these changes parallel progressive epidermal stratification and dermal maturation. Melanophores first localize in the superficial epidermis but migrate toward the dermo-epidermal junction, differentiating into iridophores and xanthophores, thereby constructing a multi-layered pigment architecture. By 90–140 dph, the epidermis becomes fully stratified, and a densely pigmented dermis is apparent, reflecting mature tissue organization. These findings support previous models proposing a spatial and temporal orchestration of chromatophore lineage deployment in teleosts [19].

Recent research on brown trout and Japanese flounder provides complementary evidence for the sequential differentiation of pigment cells, contributing to species-specific coloration patterns [20,21]. Such layered chromatophore recruitment not only governs the intensity and stability of adult pigmentation but also reflects conserved morphogenetic programs across diverse teleost lineages.

4.2. Transcriptomic Shifts Associated with Early Pigmentation Development

RNA-Seq analysis comparing 0 and 14 dph revealed a substantial transcriptional reprogramming that accompanies early pigmentation development in A. frenatus. A total of 4,245 differentially expressed genes (DEGs) were identified, with enrichment analyses pointing to upregulated processes such as pigment biosynthesis, oxidoreductase activity, and energy metabolism. The elevated expression of genes involved in oxidoreductase function and proteolytic remodeling at 25 dph indicates a metabolic transition supporting the structural and biochemical requirements for pigment deposition. This metabolic shift mirrors observations in Japanese flounder and rainbow trout, where pigmentation onset coincides with increased transcription of genes involved in melanogenesis, oxidative stress regulation, and tissue remodeling [19,20].

Concurrently, genes associated with neurogenesis, glycolysis, and early cell differentiation were downregulated, suggesting a phase-specific suppression of general developmental programs in favor of epidermal specialization and commitment to the chromatophore lineage. The repression of glycolytic genes aligns with broader trends in pigmentation studies where energy pathways are reallocated toward anabolic biosynthesis and cellular differentiation [21]. These findings support a model in which early larval stages prioritize organogenesis, followed by a transcriptionally focused shift toward skin maturation and the establishment of pigment patterns. Comparative data from teleosts, such as Oncorhynchus mykiss and Paralichthys olivaceus, highlight conserved themes in these transitions, with the developmental timing of pigmentation tightly linked to the transcriptional regulation of pigmentation-specific gene networks [20,21].

4.3. Validation of Candidate Gene Expression Patterns by qPCR

The transition from larval to juvenile pigmentation in A. frenatus marks a critical developmental window characterized by intense differentiation of chromatophores and epidermal remodeling. Our qPCR validation of transcriptomic findings at 40 dph highlighted the upregulation of key genes associated with cell proliferation, epithelial stratification, pigment transport, and endocrine regulation, corresponding temporally with the emergence of species-specific pigmentation. Specifically, the elevated expression of myc and pcna underscores increased mitotic activity and epidermal expansion, a pattern consistent with pigmentation development in other teleost models [20]. The parallel rise in p63, a pivotal transcription factor for epithelial maintenance, aligns with histological evidence of epidermal stratification. Neural crest regulatory genes pax6 and gata3 were also significantly upregulated, suggesting conservation of developmental gene networks governing pigment cell fate across teleosts. Similar transcriptional activation patterns have been reported during early pigment cell migration in zebrafish and crucian carp [22].

Moreover, genes involved in pigment vesicle trafficking, such as myo5a, myo7a, and rab38, exhibited marked increases at 40 dph. These motor proteins facilitate melanosome transport along the cytoskeleton, which is crucial for uniform pigmentation [19]. Supporting this, recent functional studies in zebrafish confirm the roles of these genes in spatial pigment deposition and stripe formation [23]. Among the most prominently upregulated was mchr1, a G protein-coupled receptor implicated in melanocortin signaling, suggesting an endocrine regulatory axis during pigmentation. This is reinforced by recent findings showing thyroid hormones and melanocortin pathways modulate melanogenesis in marine teleosts and zebrafish under variable physiological and environmental conditions [24,25]. Other ion homeostasis regulators such as trpm7 may also contribute to vesicle acidification and pigment granule dynamics during melanogenesis. Together, the qPCR data corroborate transcriptomic trends and reinforce that pigment pattern formation in A. frenatus is driven by a tightly coordinated gene expression program integrating proliferative cues, neural crest lineage commitment, vesicle transport, and hormonal modulation.

4.4. Post-Transcriptional and Translational Regulation of Pigmentation Genes

Our integrative analysis of mRNA and ribosome-protected fragment (RPF) data revealed pervasive decoupling between transcript abundance and protein synthesis during pigmentation development in A. frenatus. Genes were categorized into nine distinct regulatory profiles based on their transcriptional and translational dynamics. While some (e.g., Groups C and G) exhibited concordant regulation, the majority showed discordance, indicating significant contributions from post-transcriptional mechanisms such as mRNA stability modulation, translational efficiency variation, or protein degradation. For instance, Group C genes such as muc5a and muc5b were strongly upregulated at both transcript (logFC > 5) and RPF levels (logFC > 8), reflecting active transcriptional and translational engagement during epidermal maturation. Conversely, Group G genes like pxn1 and rl27 demonstrated parallel downregulation, potentially signaling global repression during tissue differentiation. Of particular interest were genes in Group H, such as gblp and tar1, which showed stable transcript levels but sharply reduced translation, suggesting selective translational repression—a regulatory motif increasingly recognized in teleost development [26]. Group A genes, with high mRNA expression but low protein output, may reflect active proteolysis or mRNA sequestration. In contrast, Group I genes, where translation is upregulated without corresponding transcriptional change, suggest enhanced translational efficiency or stabilization of pre-existing mRNAs.

Similar regulatory dichotomies have been observed in golden-back crucian carp and zebrafish, where microRNAs such as miR-196d and translational repressors modulate melanogenesis and pigment cell behavior [25,26]. Such findings underscore the role of small RNAs and RNA-binding proteins (RBPs) in orchestrating pigment cell differentiation through post-transcriptional mechanisms. Emerging evidence from teleost models suggests that RBPs fine-tune pigment biosynthesis pathways by selectively binding to pigment-related mRNAs, thereby controlling their stability and translational readiness [24]. Additionally, environmental and endocrine factors have been shown to influence translation efficiency in pigmentation genes, suggesting that external cues integrate into the post-transcriptional regulatory landscape [19].

Our observations are aligned with recent multi-omics studies in zebrafish and flounder, which report widespread mRNA-protein mismatches during pigment cell development and tissue remodeling [19,21]. Collectively, these findings emphasize that transcriptional output alone is insufficient to predict functional protein levels, supporting the growing consensus that translational and post-translational layers are crucial determinants of pigmentation phenotypes in teleosts.

In conclusion, our study demonstrates that pigmentation development in A. frenatus is regulated through an intricate interplay of transcriptional, translational, and post-translational mechanisms. Decoupling mRNA and protein levels for key developmental genes underscores the limitations of transcriptomic data alone in developmental studies and highlights the value of integrative omics approaches. Future work should explore the roles of RNA-binding proteins, microRNAs, and translation factors that mediate these regulatory effects, and examine whether similar post-transcriptional strategies are conserved across other teleost lineages. Additionally, investigating environmental and hormonal modulators of these pathways may reveal how external stimuli interface with intrinsic genetic programs to shape pigment pattern diversity in coral reef ecosystems.

Our results demonstrate that pigmentation development in A. frenatus is a tightly coordinated process involving multiple regulatory layers. Morphological maturation of the skin and the differentiation of pigment cells are accompanied by dynamic transcriptional reprogramming. However, the discordance between transcript and protein levels reveals an additional level of translational and post-translational regulation that fine-tunes pigment production.