Fish were caught from the Minjiang River using brails and transferred to the indoor tanks (50 × 40 × 30 cm) with a water-circulating system at 26°C under a 13L:10D lighting cycle. The healthy individuals were selected as the brood stock, and their progenies were obtained through artificial fertilization. The larva was fed with artemia three times a day, while fry and adult fish were fed twice daily with commercial feed. Seven developmental stages were identified based on skin color changes during development: 3 dpf, 15 dpf, 1 month post fertilization (mpf), 3 mpf, 6 mpf, 9 mpf, and 12 mpf. At 3 dpf and 15 dpf, the body size of the larvae is too small to directly obtain the skin tissue. Therefore, larvae which are removed with heads, viscera and yolk sacs were sampled and denoted as ST-3d and ST-15d, respectively. Similarly, dorsal skins were sampled at 1mpf, 3 mpf, 6 mpf, 9 mpf, and 12 mpf, and denoted as ST-1m, ST-3m, ST-6m, ST-9m, and ST-12m, respectively. Fish were anesthetized with MS-222 (Sigma-Aldrich, USA) before tissue sampling. All experimental protocols were approved by Southwest University and conducted under the Institutional Animal Care and Use Committee’s protocols.

The skin, scale, dorsal fin, pectoral fin, ventral fin, and anal fin of adult fish (n=3) were immersed in 0.65% normal saline and were observed under a stereomicroscope (Nikon SMZ25, Japan) to identify types of chromatophore. To investigate the dynamic density changes of melanophores and iridophores, the dorsal skin and head skin (n=9) were observed and imaged under the Nikon SMZ25 stereomicroscope at different developmental stages, viz., ST-1m, ST-3m, ST-6m, ST-9m, and ST-12m. Then, the number of melanophores per unit area of the dorsal skin and the ratio of the iridophores area to the total head area were calculated using ImageJ, respectively. Because iridophores are not obvious at 1mpf, iridophores are calculated only at 3 mpf, 6 mpf, 9 mpf, and 12 mpf. Statistical analysis of the data was performed using SPSS 26.

The pigment contents of dorsal skins at ST-1m, ST-3m, ST-6m, ST-9m, and ST-12m stages were analyzed. For melanin detection, dorsal skin samples (n=3; weigh=0.1 g/sample) were placed in a tube, then added 0.9 mL of phosphate buffer (0.1mol/L; pH=7.0). This mixture was homogenized in an ice-water bath for 10 min. 50 μL of the supernatant was detected using a fish melanin ELISA kit (Ruixin, China). The absorbance (OD) value was measured at a wavelength of 450 nm using a microplate reader (Rayto RT-6100, USA) to calculate the melanin concentration.

The contents of pteridine, β-carotene, and guanine were detected using a high-performance liquid chromatography (HPLC) system (Agilent 1200, USA). For pteridine detection, dorsal skin samples (n=3; weigh=0.3 g/sample) were ground into a powder and transferred 0.1 g of powder to a brown tube, then added 0.5% potassium hydroxide solution (5mL). This mixture was homogenized at 4°C for 5 min and oscillated for 2 min. The pH value was adjusted to 7.5 with 0.1% phosphoric acid water, followed by sonication at 4°C for 30 min and the addition of distilled water to 8mL. After oscillating mixing, 1 mL of supernatant was filtered with a 0.45 μm filter membrane, and then 20 μL of supernatant was detected at a wavelength of 360 nm using the HPLC system.

To detect β-carotene, dorsal skin samples (n=3; weigh=0.3 g/sample) were ground into a powder and placed 0.1 g of powder in a brown tube. A 0.1% butylated hydroxytoluene (BHT) -ethanol solution (5mL) was added to the tube, and the mixture oscillated at 25°C for 4h in a dark place. The resulting supernatant was transferred to a fresh brown tube, and the process was repeated once. The supernatant was collected and combined with the previously obtained supernatant. The final supernatant was mixed with 0.1% BHT-ethanol to a total volume of 10 mL. After oscillating mixing, 0.2 mL of the supernatant was filtered through a 0.22 μm filter membrane and detected at a wavelength of 445 nm using an HPLC system.

To detect guanine, dorsal skin samples (n=3; weigh=0.3 g/sample) were ground into a powder and placed 0.1 g of powder in a brown tube. Trifluoroacetic acid (5 mL) and formic acid (5 mL) were sequentially added to the tube, and the mixture oscillated at 90°C for 12min. The cooled supernatant was transferred to a fresh tube, dried under nitrogen, and reconstituted in 1mL of mobile phase consisting of 0.01 mol/L phosphate buffered salt solution (pH = 6.2) and methanol. The resulting solution was filtered through a 0.22 μm filter membrane and detected at a wavelength of 210 nm using an HPLC system.

Total RNA was extracted from larvae (n=3) at two developmental stages (ST-3d and ST-15d) and dorsal skins (n=3) at five developmental stages (ST-1m, ST-3m, ST-6m, ST-9m, and ST-12m) using Trizol reagent (Takara Bio, Japan). RNA quality and quantity were assessed using a Nanodrop ND-1000 spectrophotometer (LabTech, USA) and 2100 Bioanalyzer (Agilent Technologies, USA) as described in a previous study (Yuan et al., 2020). Samples that met the criteria (total RNA concentration ≥100 ng/μL; RIN ≥ 8) were used for sequencing library construction using a paired-end sample preparation kit (Illumina Inc., USA), and the library was sequenced on an Illumina HiSeq 2000 sequencer.

Prior to bioinformatics analysis, the raw sequencing reads underwent quality filtering using fastp (Chen et al., 2018). The resulting clean reads were aligned to the reference genome of S. taeniatus using HISAT2 (Kim et al., 2019). Gene expression quantification was calculated using RSEM (Li and Dewey, 2011). Differentially expressed genes (DEGs) analysis was conducted using DESeq2 (Love et al., 2014), with a threshold value (P value < 0.01 and |log₂FC| > 1). Enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using R package clusterProfiler 4.0 (Wu et al., 2021), and results were visualized using ggplot2 and stringr packages. To explore gene co-expression patterns and their associations with phenotypic traits (i.e., chromatophores number and pigments content), weighted gene co-expression network analysis (WGCNA) was performed using R package WGCNA (Langfelder and Horvath, 2008). Gene expression filtering was applied to retain genes with FPKM values ≥ 1 in more than half of all samples. The power parameter was set to 8, and modules were identified using a minimum module size of 30 and an abline of 0.25. Based on the expression levels of all genes and the similarity of module eigengene (ME), the biological modules were then classified. Pearson correlation analysis was conducted to study the correlation coefficient of phenotypic traits (i.e., chromatophores number and pigments content) and modules according to the ME values for modules. Modules with the squared Pearson correlation coefficient (R²) > 0.55 and P < 0.01 were considered significant, and their gene functions were further explored using GO and KEGG enrichment analysis. Intramodular connectivity between genes and intermodular correlations were calculated using the WGCNA package. Based on enrichment analysis results, genes involved in pigmentation pathways were used to construct a gene co-expression network using Cytoscape (Shannon et al., 2003), and hub genes with high connectivity were identified.

Dorsal skins from differential development were used to extract total RNA, followed by cDNA was synthesized using the PrimeScript RT Reagent Kit (Takara Bio). Gene expressions were analyzed by real-time PCR (RT-PCR) using TB Green Premix Ex TaqTM II Kit (Takara Bio) and performed on Quant Studio 1 (Thermo Fisher, USA). Primers for target genes were designed based on the genomic and transcriptomic data of S. taeniatus ( Supplementary Table A.1 ). β-actin was employed as the reference gene, and the 2^-△△Ct method was utilized to determine relative gene expression levels.

Three chromatophores including melanophore, xanthophore, and iridophore were observed in the scale, fins, and skins of adult S. taeniatus using optical microscopy. Both stellate and punctate melanophores exist in S. taeniatus. Xanthophores are dispersed near melanophores and exhibit an irregular yellow shape. Iridophores, the smallest of the three, display a punctate shape and reflect light ( Figure 1 ). Melanophores first appeared on the head and base of the ventral fin of S. taeniatus at 3 dpf, followed by iridophores on the ventral trunk at 4dpf, and finally xanthophores on the head and base of the ventral fin at 5dpf ( Figure 2 ). These three chromatophores are throughout the body of S. taeniatus, including skin, head, eye, fin, and scale.

The yolk-sac larva of S. taeniatus is transparent, gradually developing body color in the postlarva stage. Although most regions of the body are not heavily pigmented, the organ systems remain visible and fish maintains translucence at 2 months post fertilization (mpf). At 3 mpf, the peritoneum becomes covered with chromatophores presenting a silver-black color, while a dark line composed of chromatophores appears below the lateral line at 9 mpf ( Figure 3 ). The body’s pigmented regions increase with development, but the caudal region of the body remains translucent until 12 mpf ( Figure 3 ). The density of melanophores showed a significant increase with individual development, particularly at 12 mpf. However, the density of melanophores in the dorsal skins was no significant change between 6 mpf and 9 mpf. The density of Iridophores in the head skins significantly increased at 3 mpf, 6 mpf, and 9 mpf, before decreasing at 12 mpf. In addition, there is no significant change between 3 mpf and 12 mpf ( Figure 4 ). Xanthophores are obscured by a large number of melanophores, so their density change is not counted in S. taeniatus.

Additionally, the pigments secreted from three chromatophores were analyzed and quantitated. The levels of melanin, pteridine, and β-carotenoid significantly increased during individual development, with peak values observed at 12 mpf ( Figure 4 ). In contrast, the guanine content showed an initial increase at 1 and 3 mpf, followed by a slight decrease at 6 and 9 mpf, and a subsequent increase at 12 mpf ( Figure 4 ).

A systematic understanding of the transcriptional level changes during color formation was achieved through RNA-seq analysis at various developmental stages (3 dpf, 15 dpf, 1 mpf, 3 mpf, 6 mpf, 9 mpf, and 12 mpf). The raw reads were deposited in the NCBI BioProject database under accession number PRJNA949530. Approximately 957.57 Mb of clean reads were obtained after filtering. The clean reads had an average Q20 and Q30 of 98.11% and 93.99%, respectively, and approximately 90.36% of them were mapped to the reference genome of S. taeniatus ( Supplementary Table A.2 ).

DEGs were identified among seven developmental stages through pairwise transcriptome comparison, and a total of 21 comparative results were obtained. The largest number of DEGs was observed between the ST-3d group and ST-9m group. Conversely, the smallest number of DEGs was observed between the ST-15d group and ST-1m group ( Supplementary Table A.3 ). GO and KEGG functional enrichment analyses for DEGs were conducted to elucidate the predominant pathways involved in pigmentation. The results showed that DEGs were significantly enriched (p < 0.05) in ten pathways, including retinol metabolism, thyroid hormone synthesis, purine metabolism, pyrimidine metabolism, thyroid hormone signaling pathway, tyrosine metabolism, melanogenesis, melanoma, and Wnt signaling pathway ( Figure 5 ; Supplementary Tables A.4, A.5 ). Among these ten pathways, four pathways with high frequency, including retinol metabolism, thyroid hormone synthesis, purine metabolism, and pyrimidine metabolism, may be closely related to pigmentation in S. taeniatus. Additionally, DEGs from ST-3d vs. ST-15d, ST-3d vs. ST-1m, and ST-15d vs. ST-1m were significantly enriched in tyrosine metabolism. The thyroid hormone signaling pathway was significantly enriched in ST-3d vs. ST-3m, ST-3d vs. ST-6m, ST-3m vs. ST-6m, ST-6m vs. ST-12m. The melanogenesis and Wnt signaling pathways were significantly enriched in ST-6m vs. ST-12m and ST-9m vs. ST-12m. However, melanoma was only significantly enriched in ST-3d vs. ST-9m ( Figure 5A ). The previously reported pigmentation genes from these ten pathways were chosen ( Supplementary Table A.6 ), and the dynamic expression changes of these genes at different developmental stages were presented in Figure 5B .

To validate the RNA-seq results, six differentially expressed pigmentation genes were randomly selected and analyzed by RT-PCR. The gene expression patterns detected using RT-PCR were highly consistent with RNA-seq data ( Supplementary Figure A.1 ).

A total of 19 modules were obtained by WGCNA based on 21,260 genes with relative expression levels > 1 ( Figure 6 ). Four modules were significantly associated with pigments and chromatophores. The darked module was found to be significantly related to melanin and guanine content, as well as the number of melanophores and iridophores. The thistle module was significantly associated with pteridine content and melanophore count, while the darkorange2 module was significantly related to the number of iridophores. The darkmagenta module was found to be significantly associated with melanin content and iridophore count ( Figure 6 ). Go and KEGG enrichment analyses of genes from these modules revealed significant enrichment of retinol metabolism, thyroid hormone synthesis, thyroid hormone signaling pathway, melanogenesis, melanoma, pyrimidine metabolism, and Wnt signaling pathway ( Figure 7 ; Supplementary Table A.7 ). Co-expression network analysis identified two hub genes (adcy8 and rdh8) in the darked module and six hub genes (mitfa, wnt9a, wnt10b, wnt11, adcy5, and edn1) in the darkmagenta module ( Figure 8 ), highlighting their potential roles in pigmentation. However, due to insufficient gene numbers, co-expression network analysis was not conducted for the thistle and darkorange2 modules.

This study aimed to characterize changes in chromatophores and pigments at various developmental stages of S. taeniatus. Chromatophores as a complex palette form the skin colors and color patterns in fish (Sköld et al., 2016). Three types of chromatophores, black melanophore, yellow xanthophore, and silvery or blue iridophore were identified in S. taeniatus. These chromatophores form three monolayers, with iridophores on the top layer, melanophores in the middle layer, and xanthophores on the bottom layer. The appearance of melanophores on the head and base of the ventral fin was observed at 3 dpf followed by iridophores at 4 dpf and xanthophores at 5 dpf. In comparison, zebrafish displays melanophores, xanthophores and iridophores around 30 h postfertilization (hpf), 42 hpf, and 48 hpf, respectively (Kelsh, 2004). In Nile tilapia, melanophores were first observed during embryonic stages, followed by xanthophores, iridophores, and erythrophores at 6 dpf, 8 dpf, and 50 dpf, respectively (Wang et al., 2021). In summary, it can be known that the emergence of chromatophores in fish is species-specific.

Chromatophores’ density and distribution determine the body colors and patterns (Olsson et al., 2013). In the case of S. taeniatus, the larval and juvenile stages exhibit a translucence color due to the low pigmentation in most regions. The body colors of S. taeniatus are not fully developed until 12 mpf, with high concentrations of pigments such as melanin, guanine, pteridine, and carotenoid at this stage. This suggests that S. taeniatus takes longer to develop its body colors compared to other teleosts, such as the juvenile crucian carp (Liu et al., 2016) and Nile tilapia (Wang et al., 2021), which already exhibit adult-like coloration. In nature, species often develop unique phenotypic traits in response to their environment. The translucence color in S. taeniatus may serve as a camouflage strategy under ambient light conditions, reducing the animal’s silhouette and allowing it to evade predators (Zylinski and Johnsen, 2011). However, further studies are required to confirm this hypothesis.

To elucidate the mechanisms of pigmentation in S. taeniatus, RNA-seq was conducted at seven developmental stages (3 dpf, 15 dpf, 1 mpf, 3 mpf, 6 mpf, 9 mpf, and 12 mpf). The analysis revealed significant enrichment of in the tyrosine metabolism, Wnt signaling pathway, thyroid hormone synthesis, and thyroid hormone signaling pathway. Previous studies have linked the tyrosine metabolism pathway to melanin-based pigmentation in fish, which comprises genes such as tyr, tyrp1, hpdb, hgd, dct, and fah (Gan et al., 2021; Wu et al., 2022). Among these genes, tyr is a key regulator that determines the rate of melanin production (Cooksey et al., 1997), and tyrp1 plays a crucial role in the catalysis of melanin synthesis (Burgoyne et al., 2015). Notably, in larval S. taeniatus, the expression of tyr and tyrp1 decreased at 15 dpf and 30 dpf as compared to 3 dpf, indicating early inhibition of melanin synthesis, consistent with the small number of melanophores observed in the larvae. However, in the later stage of development, particularly at 12 mpf, there was a significant increase in melanophore numbers. Our findings indicate that the Wnt signaling pathway plays a dominant role in regulating melanocyte differentiation and pigmentation in fish (Yamada et al., 2010; Lim and Nusse, 2013; Takeo et al., 2016). Specifically, we observed an increase in wnt5b expression at 12 mpf, while wnt10b expression decreased. In a hybrid grouper, wnt5b expression was downregulated in the albino type compared to the wild type (Zhou et al., 2022), whereas wnt10b was found to be higher in the albino type of Astronotus ocellatus (Wang Q. et al., 2022). These results suggest that wnt5b has a positive effect on melanocyte differentiation, while wnt10b has the opposite effect. Additionally, our study revealed that wnt7b and wnt11 were upregulated and downregulated in 12 mpf, respectively, indicating their involvement in the terminal differentiation of melanocytes in S. taeniatus. However, the role of wnt7b and wnt11 in pigmentation in other teleosts remains unknown, and further studies are needed to clarify their function. Furthermore, our WGCNA results identified wnt 9a, wnt10b, and wnt11 as central regulators of pigmentation, closely associated with other pigmentation-related genes such as mitfa, adcy5, and edn1. Nevertheless, the regulatory roles of the wnt gene family in other pigmentation genes require further investigation.

Mitf is a crucial activator of pigmentation genes, such as tyr, tyrp1, and mc1r, and plays a vital regulatory role in melanocyte development, proliferation, survival, and melanin synthesis (Vachtenheim and Borovanský, 2010). Our study revealed high expression levels of mitf at 6 and 9 mpf, suggesting its involvement in melanogenesis during the adult stage of S. taeniatus. A similar expression pattern was observed in the rb1 gene. Studies have shown that Rb1 stimulated melanogenesis via Pka/Creb/Mitf signaling in humans (Lin et al., 2014). However, it remains unclear whether an interaction exists between mitfa and rb1 in fish, and further research is needed to elucidate their relationship.

The thyroid hormone plays a critical role in regulating fish pigmentation by controlling pigment cell development and patterning (McMenamin et al., 2014). In zebrafish, the absence of thyroid hormone leads to an increase in melanophore numbers and a loss of visible xanthophore and carotenoid pigments in the xanthophore (Saunders et al., 2019). Similarly, in clownfish (Amphiprion percula), the thyroid hormone regulates iridophore stripe formation and promotes the expression of duox, a gene involved in iridophore development (Salis et al., 2021). Our WGCNA study found that the thyroid hormone synthesis pathway in S. taeniatus was closely related to melanophore numbers and melanin accumulation, while the thyroid hormone signaling pathway was associated with melanophore numbers and pteridine accumulation. However, the genes involved in the thyroid hormone pathway for pigmentation in fish remain poorly understood. Previous studies in zebrafish have shown that adcy3 and adcy5 play a role in melanosome dispersion and melanocyte stripe formation (Zhang et al., 2022). Our WGCNA results identified adcy5 and adcy8 as hub genes in pigmentation of S. taeniatus. Furthermore, we observed significant changes in the expression of adcy3, adcy5, and adcy8 during color formation, and these three genes were significantly enriched in thyroid hormone synthesis ( Supplementary Table A.4 ). Therefore, we suggest that adcy3, adcy5, and adcy8 may play an essential role in pigmentation regulated by thyroid hormone. Our findings provide insights into the regulatory mechanisms of thyroid hormone on fish pigmentation. Furthermore, the thyroid hormone serves as a crucial endocrine factor that contributes to the process of metamorphosis in teleosts, which involves irreversible changes in both morphology and physiology from larvae to juvenile stages (McMenamin and Parichy, 2013). During the metamorphic stage of fish, the levels of thyroxine (T4) and triiodothyronine (T3) increase, along with markers of gene expression, such as thyroglobulin (tg), thyrotropin (tshb), deiodinase 2 (dio2) (Campinho, 2019). In our study, we observed a significant upregulation of tg expression at 6 mpf and 9 mpf ( Supplementary Table A.4 ; gene ID: Sin0120320). This finding suggests that the thyroid hormone signaling pathway not only plays a crucial role in pigmentation but also contributes to the process of metamorphosis in S. taeniatus. However, considering that thyroid hormone levels must be tightly regulated within the body, further research is necessary to elucidate the simultaneous functioning of thyroid hormones in various physiological processes.

Fish have more complex mechanisms of pigmentation than mammals due to their possession of multiple pigment cell types, unlike mammals that only have one type (melanophore). However, compared to melanin-based pigmentation, fewer studies have focused on regulatory pathways underlying non-melanin-associated pigmentation. In our study, WGCNA results showed that the retinol metabolism pathway was also highly related to pigmentation. KEGG analyses showed that carotenoid-associated genes, such as beta-carotene oxygenase 1 (bco1), were significantly enriched in this pathway. Studies on salmon have shown that bco1 and bco2 expression are significantly associated with the red coloration of flesh (Helgeland et al., 2019; Lehnert et al., 2019). In S. taeniatus, bco1 was observed to be highly expressed at 1 mpf, indicating its involvement in carotenoid-based pigmentation during the larval stage. However, at the adult stage, despite the increased generation of xanthophores, the expression of bco1 significantly decreased in S. taeniatus. We suggest that adult S. taeniatus are more likely to absorb carotenoids from their foods. However, we did not observe any significant expression changes in the expression of bco2 concerning body coloration in S. taeniatus. Xanthophores are responsible for the synthesis and storage of both carotenoids and pteridines. Xdh is a key enzyme responsible for pteridine biosynthesis (Ziegler, 2003). Studies on Carassius auratus have identified xdh as a key driver of variation in skin color (Zhang et al., 2017). In our study, xdh was observed to have a high expression level at the early stage of S. taeniatus, coinciding with the accumulation of xanthophores. However, we did not find any reported genes involved in pteridine-based pigmentation with a high expression level during later developmental stages, despite the high levels of pteridine accumulation. Therefore, we speculate that some novel genes may participate in pteridine-based pigmentation during late stages, and further studies are needed to explore new functional genes.

Iridophores reflect light through platelets derived from guanine to generate iridescent colors. In our study, we observed the expression changes of genes in the purine metabolism pathway were observed accompanying body color formation in S. taeniatus. The purine nucleoside phosphorylase (pnp) gene family (i.e., pnp4a, pnp4b, pnp5a, pnp5b, and pnp6) occupied a major part of the purine metabolism pathway. However, these genes exhibited distinct gene expression patterns at different developmental stages of S. taeniatus. Pnp5a and pnp6 had a high expression level during larval/early stages, while pnp4a, pnp4b, and pnp5b were primarily expressed during late stages. This suggests that genes of pnp family may play different roles in the stage-specific pigmentation of S. taeniatus. Previous studies have shown that pnp4a knockout leads to a significant decrease in iridophores, indicating its important role in iridophores-based pigmentation (Kimura et al., 2017). Transcriptome analysis of a body color mutant showed that pnp5a is involved in the pigmentation of a carp-goldfish hybrid (Zhou et al., 2019). However, the roles of pnp4b, pnp5b, and pnp6 in pigmentation are unclear in fish. Further study is required to reveal the function of pnp genes in pigmentation and elucidate the relationship between these genes.