Animal care and experimental protocols were approved for this study by the Guangxi University Ethical Committee under the approval document (GXU-2021-126). To investigate the adaptive mechanism for cave-dwelling in the Sinocyclocheilus radiation, 9 species were collected from caves in Yunnan, Guizhou and Guangxi of China during 2019 – 2021 as a part of an ongoing phylogenomic study. Cavefish were transported alive to the laboratory in oxygenated plastic bags in a cooler box to both keep the temperatures low (18 – 20°C) and to provide darkness. We let them settle down prior to initiation of the experiments. Following morphological observations (mentioned below), they were anesthetized using 30 mg/L MS-222 (ethyl 3-aminobenzoate methanesulfonate; Sigma-Aldrich). The skin tissue was taken from the right side near the back of each individual in approximately 0.8 × 1.0 cm under sterile conditions. Three biological repeats were taken for each species. The tissue extraction was done in a DNA/RNA-free clean room. Biopsied tissues were placed in RNAlater and stored in an ultra-low temperature freezer at -83°C prior to further analysis.

Given their rarity, sampling difficulty, and the need for a representation of the diversification, we enhanced taxon sampling. By adding skin transcriptomes for five additional Sinocyclocheilus species from previous studies and supplementing the morphological data with micrographs of Sinocyclocheilus qiubeiensis and Sinocyclocheilus tianlinensis (Yang et al., 2016; Li, 2019; Li et al. 2020). Overall, the dataset included species representing the four main clades in the context of Mao et al. (Mao et al., 2021), as well as the three main eye-types/habitats (normal-eyed/surface, SU; normal-eyed/stygophilic, SP; micro-eyed/stygophilic, SP; eyeless/stygobitic, SB), representative of the Sinocyclocheilus radiation ( Figure 1 ; Supplementary Table S1 ). The eye-type in this study serves as a practical, easily observable feature that corresponds with the habitat occupation of these species.

In order to generate transcription data for Sinocyclocheilus skin from the morphotypes, the total RNA was extracted from the 27 skin samples (9 species that we sampled) using Trizol methods; mRNA was captured from total RNA through Oligo (dT). To synthesize cDNA, fragmented mRNA was utilized as a template in the M-MuLV reverse transcriptase system. Random oligonucleotides acted as primers to initiate the synthesis of the first cDNA strand. RNaseH was then used to degrade the RNA strand, and dNTPs were added under the DNA polymerase I system to synthesize the second cDNA strand. The library was constructed using the NEBNext^® Ultra™ RNA Library Prep Kit (Illumina, USA). The RNA quality and concentration were measured using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) and Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). The effective concentration of the library was again accurately quantified by qPCR to ensure the quality of the library. An Illumina NovaSeq 6000 was used to sequence the library, producing paired-end reads of 150 bp each. The library construction and sequencing were conducted at Beijing Novogene Bioinformatics Technology Co., Ltd.

In addition to our data, RNA-seq data for 5 species were obtained from two public databases. We downloaded transcriptome datasets for 3 species from the public DRYAD database: Sinocyclocheilus oxycephalus, Sinocyclocheilus tianlinensis, and Sinocyclocheilus qiubeiensis (https://datadryad.org/stash/share/t5cZIXoVUgyhpzEP6z-GN6xjc5EU3TvPPEwbdIo7siI), and transcriptome datasets for 2 species: Sinocyclocheilus grahami (NCBI: SRR2960332) and Sinocyclocheilus anshuiensis (NCBI: SRR2960751) from the NCBI database. The skin transcriptomes from the 9 species from our study and the 5 additional species were evaluated with FastQC v0.11.9 (Andrews, 2010). To ensure the accuracy and reliability of our data, we implemented a rigorous filtering process. This involved removing reads containing adapters, poly-N sequences, and those with low-quality scores. The resulting high-quality, clean reads formed the basis for all downstream analyses. The clean reads are deposited in the NCBI database and will be made public upon acceptance (NCBI SRA accession number: SUB13474193).

We used Trinity v2.8.6 (Grabherr et al., 2011) to assemble clean reads for each species and extracted the longest transcripts (unigenes), to obtain single gene sequences. Subsequently, the CD-HIT v4.6.8 (Fu et al., 2012) was used with a 95% threshold to cluster sequences and eliminate redundancy in the final assembly. To predict the full open reading frames (ORFs) for each gene, we used TransDecoder (http://transdecoder.github.io/) default parameters. We then annotated the resulting protein sequences using BLASTP (Blast+ v2.6.0) with an E-value of 10^-5 and using the UniProt database (v2022_05) (Mahram and Herbordt, 2015). We retained only UniProt germplasm from the top 10 hits for each species and selected the final UniProt ID based on the number of species that used it as the best match. Using OrthoFinder v1.1.2. (Emms and Kelly, 2015), we performed gene clustering while retaining only those containing at least one transcript per species. These filtering steps resulted in 7,374 annotated orthologous groups for downstream analysis of expression differences.

For phylogenetic inference, we utilized single-copy gene families obtained from OrthoFinder, as previously described. We discarded gene families comprising sequences shorter than 200 amino acids. The retained amino acid sequence families were aligned with Muscle 5 (Edgar, 2022), with default parameters. To excise ambiguously aligned regions, we used Gblock v0.91b (Castresana, 2000). We conducted a maximum likelihood (ML) phylogenetic analysis on a refined set of 1,369 single-copy orthogroups using IQ-TREE 2. We selected the best-fit model, GTR+F+I+I+R3, which encompasses 72,820 distinct patterns, 77,156 parsimony-informative sites, 76,466 singleton sites, and 1,318,361 constant sites, to construct an unrooted phylogeny. To determine node support, we performed an ML bootstrap analysis using 1,000 pseudoreplicates over 102 iterations. (Minh et al., 2020). The resulting ML phylogenetic tree (best log likelihood score: -3,627,797.972) was visualized using Figtree (Rambaut, 2009). We explored two methods to root the tree. First, we used mid-point rooting, which places the root of the tree at the midpoint between the two most distant taxa, making the fewest assumptions about the evolutionary rates in different lineages. Second, we rooted the tree using a well-established early-merging Clade A species, S. guyangensis (Mao et al., 2021; Mao T. et al., 2022).

To understand the functional genes that may facilitate the adaptation of Sinocyclocheilus to its environment, orthologous of the skin of 14 species were tested for signals of positive selection. The dN, dS, and dN/dS values of the orthologous were calculated using the CodeML program of the PAML package (Yang, 2007). The CodeML parameter was set to “Runmode = 0, Model = 0”. Then, ω> 1 can be judged to have experienced positive selection pressure effects in this gene, and a total of 329 genes were screened. KEGG enrichment of PSGs was performed using KOBAS v2.0.12, and the PSGs in significantly enriched pathways were annotated based on the Evolutionary Genealogy of Genes: Nonsupervised Orthologous Groups (eggNog) database (Xie et al., 2011; Huerta-Cepas et al., 2016).

To study how orthologous genes are expressed in three different habitats - SB, SP, and SU - we conducted a mapping analysis of high-quality reads of each species against its corresponding representative orthogroups obtained earlier; for this, we used Bowtie 2 v2.2.9 for direct homology filtering (default settings) (Langmead and Salzberg, 2012), and for gene expression levels analysis we used RSEM v1.2.26 (default setting) (Li and Dewey, 2011). To analyze the gene expression profiles of each species, we utilized the ggbiplot package to conduct a principal component analysis (PCA) following the method outlined by Yeung and Ruzzo (2001). The ggplot2 package was used to statistically generate stacked bar graphs with clustered trees reflecting the similarity between species and habitat groups, as well as information on the expression profiles of the orthologous genes for each species (Thiergart et al., 2020). We then identified differential expression in the orthogroups between species using edgeR and adjusted the resulting P-values through the application of the Benjamini and Hochberg method, which effectively estimates the false discovery rate (FDR). We considered differentially expressed genes to be those with absolute log2 fold change (|log2FC|) greater than 1 and FDR less than 0.05. These thresholds were chosen to determine statistical significance. Finally, GOseq and KOBAS v2.0.12 were used to perform GO enrichment analysis and KEGG enrichment analysis of differentially expressed orthogroups (Kanehisa et al., 2008; Young et al., 2010). P-values were obtained using a hypergeometric test and the significance term with an adjusted P-value threshold of 0.05. Using these methods, we identified significantly enriched biological pathways, including important biochemical metabolic and signal transduction pathways across the orthogroups.

Using a digital camera (Canon EOS 6D Mark II AF-A) set at a fixed distance of 0.3m from the tank and LUX = 45 for ambient light, we captured images with the following settings: Shutter speed: 1/250s; F/20; ISO 200. Next, the fish were returned to their aquarium system (pH: 7.0-8.0; temperature: 19 ± 1°C; dissolved oxygen: 8.5 mg/L) devoid of light to maintain their natural skin characteristics. Subsequently, we examined the pigmentation and scales in samples preserved in 70% ethanol using a Leica M165FC stereomicroscope.

To compare the body color differences of the Sinocyclocheilus species, we used the Color Summarizer v0.8 (http://mkweb.bcgsc.ca/colorsummarizer/analyze) to digitize the average color values of the skin. In order to ascertain the average red-green-blue (RGB) color component of the image, we employed a systematic approach. Firstly, a set of reference colors that had the least deviation from the image colors were identified, which allowed to discern variation in color with greater precision. Specifically, 5 reference colors denoted by “k” were selected. Subsequently, 5 representative color ratios from each sample taken. Finally, the average RGB color values were calculated from these ratios to obtain an overall RGB color component of the image. This procedure provided an accurate and reliable estimation of the RGB color values. This process was repeated three times separately for the three samples collected to obtain a representative color for each Sinocyclocheilus species. Our analysis focused on the 5 most representative color ratios in order to calculate the proportion of darker blotches. For instance, in a particular species, mean color values of RGB (250, 250, 250) were observed in the white areas while the darker black areas had color values closer to RGB (0, 0, 0). We postulated that a higher concentration of melanin deposition is associated with a decrease in RGB values towards RGB (0, 0, 0). The area and number of lateral line scales were counted for each sample using ImageJ (Abràmoff et al., 2004). The scale sizes were classified into three categories: no-scales, small-scales, and large-scales, and the degree of scale cover was classified as full-cover, no-cover, and partial-cover.

To understand the molecular mechanisms underlying skin coloration in Sinocyclocheilus species, we generated high-quality transcriptome sequencing data from skin tissue samples of 9 distinct Sinocyclocheilus species ( Supplementary Table S2 ). Following this, a de novo assembly was performed for sequences from the 14 species. Among them, S. oxycephalus (205,993) has the highest number of transcripts and S. grahami (80,118) has the lowest. The average transcript length post-assembly varied between 856 and 1,158 bps, while the N50 length ranged from 1,621 to 2,414 bps. The OrthoFinder identified a total of 384,146 orthologous genes in the transcripts. Overall, orthogroups accounted for 94.3% of the genes. To facilitate comparative analysis among the various species, we focused on the 7,374 orthologous genes present in all examined species. These orthologous gene families were clustered, and 1,369 single-copy orthologues genes were obtained after multiple sequence alignment and low quality pairwise pruning. These orthologous genes were subjected to quantitative comparative analysis.

A total of 1,348 differentially expressed genes (DEGs) were identified on the skins of Sinocyclocheilus species by conducting a pairwise comparison across the three habitats ( Figure 2 ; Supplementary Table S3 ). When compared to SUs, SPs exhibited 200 up-regulated and 467 down-regulated genes; SBs displayed 292 up-regulated and 394 down-regulated genes. The majority of DEGs were up-regulated in SUs and the highest number of down-regulated DEGs was observed in SPs ( Figure 2A ). We also found more unique gene up-regulation in SP vs. SB, and that more of their shared genes were consistently down-regulated compared to SUs (SU vs. SP and SU vs. SB, Figure 2B ). This suggested a significant impact of cave habitats on gene expression patterns.

GO enrichment analysis of 1,348 DEGs in different eye-types/habitats identified 630, 565 and 799 significantly enriched GO terms in SU vs. SP, SP vs. SB and SU vs. SB, respectively (adjusted P<0.05) ( Supplementary Table S4 ). These DEGs were involved in functional responses mainly related to stimulus responses, catalytic activity, multi-organism process, and immune system process ( Figure 3A ). We found that the co-enriched terms in these three groups were energy metabolism related ( Supplementary Table S5 ). We also found two major categories of GO terms have been enriched in SU vs. cave dwellers (SP and SB). Stimulus response-related terms such as external biotic stimulus, bacterial, fungal, and viral reactions, and numerous immune-related terms. SU vs. SP was mainly related to cellular immunity, such as leukocyte mediated immunity (GO:0002444), regulation of macrophage derived foam cell differentiation (GO:0010743); while SU vs. SB increased the regulation of apoptosis, such as: positive regulation of MAPK cascade (GO:0043410), regulation of neuron apoptotic process (GO:0043523), regulation of epithelial cell apoptotic process (GO:1904035). Anyhow, this suggested that these observations may be due to the differences in the exogenous biological stimulation of surface water environment. Cave species seemed to have different immune strategies to these external stimuli.

Genes associated with changes in oxygen levels showed differences in expression across different habitats, for example, response to hypoxia (GO:0001666), response to decreased oxygen levels (GO:0036293), cellular response to hypoxia (GO:0071456), etc. Interestingly, only in SPs and SBs can respiratory electron transport chain (GO:0022904), mitochondrial respirasome (GO:0005746), respiratory electron transport chain (GO:0022904), and mitochondrial respirasome (GO:0005746), respiratory chain complex (GO:0098803) and response to pH (GO:0009268) be identified ( Supplementary Table S5 ). This may reflect the different cellular and mitochondrial respiration efficiency of cave dwellers.

The number of significantly enriched pathways in three groups was 30, 38, and 22, respectively, and the same pathways were mainly related to metabolisms, such as pentose phosphate pathway (ko00030), biosynthesis of amino acids (ko01130), Carbon metabolism (ko01200), PPAR signaling pathway (ko03320) and microbial metabolism in diverse environments (ko01120), etc. Secondly, immune response-related pathways were phagosome (ko04145) (adjusted P<0.05) ( Figure 3B ; Supplementary Table S6 ). This suggested that the strong influence of the habitat environment leads to differences in skin metabolism and skin microbial metabolism among different habitat populations. Compared to SPs, all DEGs (28) enriched in the PPAR signaling pathway and carbon metabolism pathway were up-regulated in SUs, followed by 16 DGEs in SBs. However, all DEGs (9) enriched oxidative phosphorylation (OXPHOS, ko00190) pathway and 5 DEGs enriched fatty acid degradation pathway (ko00071) were upregulated in SPs. This may indicate differences in the regulation of energy and mitochondrial metabolism between species in different habitats ( Supplementary Table S7 ). These pathways in different eye-types/habitats showed that there were common molecular mechanisms to understand their habitual differences and evidence of genetic variation and transcriptional plasticity or adaptation in response to environmental change.

In particular, microbial immune-related pathways such as complement and coagulation cascades (ko04610), intestinal immune network for IgA production (ko04672), and viral protein interaction with cytokine and cytokine receptors (ko04061) found in SU vs. cave dwellers (SP and SB) ( Supplementary Table S6 ). This may indicate differences in immune responses to microorganisms between SUs and cave dwellers. Moreover, hematopoietic cell lineage (ko04640) and ECM-receptor interaction (ko04512), and pathways associated with oxidative stress (OXPHOS, glutathione metabolism, cysteine and methionine metabolism, metabolism of xenobiotics by cytochrome P450, drug metabolism - cytochrome P450) exist only in SP vs. SB. This may indicate a widespread oxidative stress response of cave dwellers to stimuli of the cave environments.

We also found that phenylalanine, tyrosine and tryptophan biosynthesis pathway (ko00400), and especially the tyrosine metabolism pathway (ko00350) may be influenced by their habitats, and affected melanin differences ( Supplementary Table S6 ).

After analyzing the KEGG pathways of these 329 positive selection genes (PSGs), we found that genes under positive selection were most significantly enriched in six metabolic pathways ( Table 1 ). They were associated with viral infection, signaling, and cell necrosis and apoptosis. A total of 18 PSGs were identified in these enriched pathways, among which the tumor protein p63 regulatory 1 (OG0018161) was under the strongest selection pressure ( Supplementary Table S8 ). These genes may be involved in the adaptation process of Sinocyclocheilus to cave dwelling. We also found that Reticulocalbin 3, EF-hand calcium binding domain (OG0016528), S100 calcium binding protein U (OG0016522), and positive regulation of vitamin D (OG0016683) have all been implicated in calcium regulation.

We filtered the clusters that comprised a single sequence from each of the 14 transcriptomes and retrieved 1,369 putative single-copy orthologous genes. We concatenated and aligned these genes into a supermatrix 1,471,983 informative sites for the 14 taxa. The two orthologous gene-based maximum likelihood trees, rooted using both the mid-point method and S. guanyangensis as the outgroup, showed two similar trees, each with five well-supported clades ( Figure 4 ; Supplementary Figure S1 ). Three notable discrepancies with previous phylogenies were the positions of early diverging Sinocyclocheilus lineages: S. xunlensis, S. oxycephalus, and S. furcodorsalis. Our phylogenetic tree unambiguously showed that S. oxycephalus was a separate lineage, while S. xunlensis was recovered as the sister group of S. guilinensis; four species including S. furcodorsalis constituted a separate lineage.

The analysis of color variation in the 14 Sinocyclocheilus species highlighted a rich diversity in body coloration and patterns, with each species exhibiting unique colors, mainly in a combination of pinkish-white, gray, and yellow ( Figure 5 ; Supplementary Figure S2 ).

For every species, color attributes were assessed based on average skin color value and the proportion of dark blotches and distinct spots on the skin. Sinocyclocheilus tianlinensis had the highest skin RGB value (195, 181, 177), followed by S. anshuiensis (191, 179, 177) and S. furcodorsalis (188, 163, 155); S. grahami (138, 123, 111), S. oxycephalus (142, 116, 90), and S. qiubeiensis (142, 127, 104) displayed the lowest RGB values ( Supplementary Table S9 ). The RGB values for black blotches in these species were under (120, 108, 114) with lower values signifying darker colors, such as dark blotches and spots ( Figure 5 ; Supplementary Table S9 ).

Skin colors were associated with eye-types and habitats. Normal-eyed SUs had yellow or gray skin, while most SPs exhibited gray or pale gray tones. Unique spots were also indicative of eye-type and habitat groups; normal-eyed S. grahami (SU) had the highest proportion of black spots, succeeded by normal-eyed S. oxycephalus (SU) and normal-eyed S. qiubeiensis (SU). Eyeless SBs and micro-eyed S. xunlensis (SP), which displayed depigmented pink-white skin, had the lowest percentage of dark blotches ( Figure 5 ; Supplementary Table S9 ). However, microscopic observations showed that SBs and S. xunlensis (SP) still had numerous tiny black blotches dispersed across their skin ( Supplementary Figure S2 ); comparatively, eyeless SBs were pink with fewer pronounced dark blotches and distinct spots.

When examining color in the context of phylogeny, Clades I-V displayed various combinations of skin color and pigmentation relative to eye-types and habitats ( Figure 5 ; Supplementary Figure S2 ). Notably, in Clade III and Clade IV, which contain the largest proportion of exclusive cave dwellers, SBs were pink with fewer dark blotches. Hence, the Sinocyclocheilus genus adapted to cave environments and evolved with convergent skin coloration in different evolutionary clades.

Microscopic examination of scale characteristics showed a marked reduction in scale size and coverage in 11 Sinocyclocheilus species except for three species (S. mashanensis: mean ± SD: 3.74 ± 1.36, S. zhenfengensis: mean ± SD: 2.93 ± 0.95, S. guilinensis: mean ± SD: 2.00 ±0.62), which exhibited complete coverage by large scales, and the mean number of their lateral line scales are respectively 49 ± 3, 43 ± 4 and 50±8. S. oxycephalus (mean ± SD: 0.17 ± 0.11), S. qiubeiensis (mean ± SD: 0.61 ± 0.28), and S. grahami (mean ± SD: 0.72 ± 0.40) had smaller scales, also with more lateral line scales (mean ± SD: 70 ± 3, 78 ± 3 and 72 ± 2), but most of the body scales were buried in the skin and disappeared ( Figure 5 ; Supplementary Figure S2 ). We found that larger scales and fewer lateral line scales were mainly present in Clade III, IV while smaller scales and more lateral line scales species were mainly found in Clade V. The reduction of scale coverage can be found in all clades (I-V). This may indicate convergent evolution of scale reduction in Sinocyclocheilus.

Generally, SBs exhibited smaller and fewer scales than SPs. The scales of S. tianlinensis (SB) were absent altogether, and S. anshuiensis (SB) had only rare scales. These species in Clade III, Clade IV revealed a pattern of scales, characterized by a gradual reduction in scale size, decreased coverage and the number of lateral line scales, and eventual disappearance, which may be related to cave adaptation. Clade V species displayed variability in scale morphology and coverage, suggesting a chaotic pattern of scale evolution in Clade V unrelated to habitat. Here, the small scales were widely spaced ( Figure 5 ).

The PCA analysis and cluster analysis of the 14 Sinocyclocheilus species using expression data from 7,374 orthologous genes revealed that 13 species clustered together, separate from S. angustiporus ( Figure 6A ; Supplementary Table S10 ). In addition, micro-eyed SPs and eyeless SBs showed differences in gene expression compared with normal-eyed SUs ( Figure 6B ; Supplementary Table S10 ). This suggested that patterns of gene expression were similar between these species, but the effects of habitats can still be seen.

Interestingly, the second principal component (variance explained 16.55%) differentiated species with and without black blotches ( Figure 6A ). Moreover, S. tianeensis, S. cyphotergous and S. xunlensis with little black blotches were clustered in one branch, but expression patterns vary among species with full-cover large scales ( Figure 6B ). This suggested that changes in gene expression patterns may affect pigmentation similarly but not scales.

In our analysis of DEGs, we identified distinct patterns pertaining to metabolism, oxidative stress, and immune responses in various Sinocyclocheilus species. These patterns appear to be correlated with their specific habitats. Such results substantiate the presence of cave-environmental gradients within natural ecosystems. In line with our findings, several studies have pointed toward physiological and metabolic distinctions between surface dwelling and cave dwelling species (Yang et al., 2016; Stahl and Gross, 2017; Krishnan et al., 2020; Boggs and Gross, 2021; Medley et al., 2022).

These variances could be due to differences in habitat characteristics like light exposure, ambient oxygen levels, and nutrient resource availability. Furthermore, the effects of different experimental sampling methods on gene expression has also been shown (Tessler et al., 2020). The intricate interplay among these factors significantly intensifies the complexity of environmental gradients, making the prediction of gene expression responses more challenging (Passow et al., 2017).

While the analysis of coloration and large-scale gene expression patterns may provide some insights, they fall short of identifying the specific environmental factors driving convergent evolution in these habitats. A recent comprehensive genomic study encompassing 248 species of A. mexicanus from diverse habitats shed light on the underlying mechanisms for convergent evolution in cavefish (Moran et al., 2023). Likewise, changes in fish transcriptome levels have been linked to phenotypic responses to environmental extremes (Wang and Guo, 2019). However, solely relying on this approach provides insights that might be limited to the specific context of the study in question. A potentially more effective strategy would combine pinpointing and investigating differentially expressed genes (DEGs) with functional annotations based on a priori hypotheses. Such an integrated approach could provide deeper insights into the range of environmental stressors within a habitat and illuminate the evolutionary pathways leading to different phenotypes.

Our ML trees, based on 1,369 unique single-copy orthologous genes, showed some discordance when compared to recent mitochondrial gene (mt-DNA) and RAD-seq based studies (Zhao and Zhang, 2009; Mao et al., 2021; Mao T. et al., 2022; Wen et al., 2022; Xu et al., 2023). The robust node support in our phylogeny might be attributed to the abundance of genes under selection derived from transcriptome sequencing. However, this precision may have a downside: the unique genes associated with the transcriptome we used could potentially not be under neutral evolution, thereby causing the phylogeny of expressed genes important for the functions of the skin. Most notably, Clade I of our study contains species from three clades that employed mt-DNA/RAD-seq data. This was possibly due to shared skin traits. In addition, S. oxycephalus formed a distinct divergent lineage (II). The skin of these species had a distinctive gray color. And previous studies have shown that its skin is charcoal gray skin unlike other types (Li et al., 2020). Furthermore, we identified S. furcodorsalis, which was from clade B of previous studies placed in our Clade V, which contained predominantly species from clade D (which are mostly SUs) (Mao et al., 2021). This suggested that S. furcodorsalis may have a skin transcription profile similar to that of SUs. In fact, this species, despite being a cave-dweller, enigmatically had some of the most prominent scales. However, our phylogeny lacked high taxon sampling, which prevented us from making deeper inferences on transcriptional plasticity. However, this study highlights the significance and precision of transcriptomic data in understanding the evolutionary relationships among Sinocyclocheilus species.

It is known that the phenotypic adaptations of cave-dwelling in Sinocyclocheilus cavefish involve changes in eye types and color and pigmentation patterns (Meng et al., 2013; Li et al., 2020); our results from the context of skin morphology and gene expression confirm this. We found that eyeless SBs and some micro-eyed SPs without obvious black blotches showing similar expression patterns, even though they were in different clades ( Figure 5 ). Existing literature has established that changes in skin color, driven primarily by adaptation to diverse habitats, serve as a major contributing factor to differential gene expression patterns (Hinaux et al., 2013; Stahl and Gross, 2015; Gross et al., 2016; Li et al., 2020). Skin color adaptations in fish are known to occur rapidly when they are introduced to novel environments (Leclercq et al., 2010; Nilsson Sköld et al., 2013).

The relationship between fish skin color and melanin content was examined, focusing on key genes involved in phenylalanine, tyrosine and tryptophan biosynthesis pathway, and especially the tyrosine metabolism ( Figure 7 ). We found that the tyrosine metabolic pathway, influenced by habitat, played a role in determining melanosis in the skin of Sinocyclocheilus. Previous research has demonstrated that decreased phenylalanine-4-hydroxylase (PAH, OG: 0000390) activity can affect Try activity (Leandro et al., 2017), and inhibition of phenylalanine conversion to tyrosine will lead to reduced melanin production (Staudigl et al., 2011; Zhou et al., 2021). However, we found low expression of the gene encoding PAH in these species. This may be due to the lack of phenylalanine and protein catabolic inputs in the diet of cavefish (Borowsky, 2018). However, L-phenylalanine (L-Phe) is also an essential protein-producing amino acid, and it is important to avoid it being completely catabolic (Xiao et al., 2020). Therefore, in order to achieve the dual role of effective preservation and removal of excess L-Phe, this regulatory mechanism may be more favorable to the adaptation of their environment.

Our analysis revealed that SBs had the fewest black blotches and the lowest expression of tyrosine metabolism-related enzymes: aspartate aminotransferase (AST, OG0001823), fumarylacetoacetase (FAH, OG0006238; Figures 5 ; 7 ). This indicated that phenylalanine and tyrosine metabolism was inhibited, and the final melanin synthesis was affected. Consistent with previous research, melanin synthesis is triggered by sunlight, which can offer benefits such as camouflage, protection from ultraviolet radiation, and a role in social signaling (Rzepka et al., 2016). Thus, whereas melanin synthesis is inhibited in SBs that inhabit lightless environments.

SPs such as S. guanyangensis and S. cyphotergous had elevated AST activity. The way these species produce melanin may be related to the accumulation of homogentisic acid (HGA). HGA oxidation leads to melanin-like pigmentation (Giustarini et al., 2012). Interestingly, cultivation of these species in our laboratory revealed that SPs seemed to exhibit an enhanced ability to adapt rapidly to fluctuating light conditions. This may enable them to adapt to the changing light environment between cave and semi-enclosed cave environments. Therefore, AST may have an important role in fish skin differentiation and variation through the melanogenesis pathway. The observed pigmentation differences among Sinocyclocheilus species may represent adaptations to their specific environments. Reduced pigmentation could be advantageous in the dark cave environment, while increased pigmentation may be more beneficial for survival in semi-enclosed and surface environments (Romero and Green, 2005; Soares and Niemiller, 2020).

Degeneration of scales was found in all clades in our study ( Figure 4 ). It is not unexpected to find skin scale degradation prevalent in Sinocyclocheilus, as it is a trait seen in some other Carpiformes, unrelated to cave environments (Harris et al., 2008; Zhao and Zhang, 2009; Zhu et al., 2019). The scales of Sinocyclocheilus are known to be bony scales, and the posterior part of their scales consists of bone (Zhao and Zhang, 2009). We found that Reticulocalbin 3, EF-hand calcium binding domain (OG0016528), S100 calcium binding protein U (OG0016522), and Positive regulation of vitamin D receptor signaling pathway (OG0016683) were under positive selection, which may contribute to the calcium-bone homeostasis of fish skin (Schäfer and Heizmann, 1996; Bouillon and Suda, 2014; Girard et al., 2015). And PSGs are enriched in the MAPK signaling pathway, which is known to be associated with scale development, so the role of such genes involved in scale evolution is strong (Zhou Q. L. et al., 2022); there is also a strong association between scale regression and cave-dwelling. However, a deeper analysis is warranted to explore the connection between cave environments and scale features perse, i.e. a comparison between scaleless Sinocyclocheilus and scaleless other Carpiformes. Our research showcases the remarkable adaptability of cave-dwelling species to dynamic environments, offering insights into fish adaptation processes in response to anthropogenic alterations, such as habitat degradation and climate change. Sinocyclocheilus cavefish share few positive selections among their skin orthologous genes, and genes associated with the limited skin phenotype are more difficult to detect. The most likely plausible reason for this evolutionary phenomenon is the low rate of evolution of the cavefish genome and genome-wide relaxation of selection (Torres-Paz et al., 2018; Policarpo et al., 2020; Zhao et al., 2022). Our study demonstrates the remarkable adaptability of cave-dwelling species to dynamic environments and provides insights into the process of fish adaptation to anthropogenic changes such as habitat degradation and climate change.

PSGs in Sinocyclocheilus were enriched in herpes simplex virus 1 (HSV-1) infection, cytokine-cytokine receptor interactions, p53 signaling pathway, necrosis, and apoptosis related pathways. The current study has found that under prolonged environmental stress (e.g. hypoxia, prolonged hypothermia, starvation, etc.), facilitation and evolution of pathways such as p53 signaling pathway, MAPK signaling pathway, cytokine-cytokine receptor interaction, and apoptosis related pathways are prone to occur and echoes our study (Tong et al., 2021; Mao J. et al., 2022; Voskarides et al., 2022; Li et al., 2023). Our study highlighted that an important aspect of skin evolution is resistance to viral infection. These pathways promote coordination between cells and contribute to the clearance of pathogenic infections. Previous studies have shown that p53 plays a dual role in the replication of HSV-1 at different stages of infection (Sato and Tsurumi, 2013; Voskarides et al., 2022). Our study found that one of the p53 homologs: the Tumor protein p63 regulated 1 (OG0018161) was subject to the strongest positive selection (Aloni-Grinstein et al., 2018). Interestingly, it has been shown that HSV-1 manages to counteract this negative effect of p53 through the viral protein ICP22, which is able to bind p53 directly and eliminate its function (Sato and Tsurumi, 2013). Like many other viruses, herpesviruses their hosts co-evolved (Chawla et al., 2022). Fish skin is the first line of defense against viral and pathogenic microbial infections in the aquatic environment, and the adaptive evolution of immune and apoptosis-related genes may have facilitated the development of immune function in the skin of Sinocyclocheilus, enabling their ability to rapidly adapt to new environments and to rapidly occupy vacant ecological niches.