Comparative transcriptomics reveals the molecular genetic basis of pigmentation loss in Sinocyclocheilus cavefishes

Abstract Cave‐dwelling animals evolve distinct troglomorphic traits, such as loss of eyes, skin pigmentation, and augmentation of senses following long‐term adaptation to perpetual darkness. However, the molecular genetic mechanisms underlying these phenotypic variations remain unclear. In this study, we conducted comparative histology and comparative transcriptomics study of the skin of eight Sinocyclocheilus species (Cypriniformes: Cyprinidae) that included surface‐ and cave‐dwelling species. We analyzed four surface and four cavefish species by using next‐generation sequencing, and a total of 802,798,907 clean reads were generated and assembled into 505,495,009 transcripts, which contributed to 1,037,334 unigenes. Bioinformatic comparisons revealed 10,629 and 6,442 significantly differentially expressed unigenes between four different surface‐cave fish groups. Further, tens of differentially expressed genes (DEGs) potentially related to skin pigmentation were identified. Most of these DEGs (including GNAQ, PKA, NRAS, and p38) are downregulated in cavefish species. They are involved in key signaling pathways of pigment synthesis, such as the melanogenesis, Wnt, and MAPK pathways. This trend of downregulation was confirmed through qPCR experiments. This study will deepen our understanding of the formation of troglomorphic traits in cavefishes.

to perform critical functions in pigmentation in many fish species (Parichy, 2006).
The wild freshwater teleost genus Sinocyclocheilus (Cypriniformes: Cyprinidae) is a treasured endemic fish genus, distributed in the karst regions of the east and northwest Yungui Plateau, Guangxi, southwestern China. This genus is composed of more than 55 know species (Meng et al., 2013). Owing to distinct differences in terms of phenotype and habitat, both morphotypes (cave and surface) exist within one genus. High species diversity and skin phenotypic variation make Sinocyclocheilus particularly suitable for studying the molecular genetic mechanisms underlying the evolution of differential pigmentation. Previous studies on three Sinocyclocheilus species (Sinocyclocheilus graham [surface fish], Sinocyclocheilus rhinocerous [cavefish], and Sinocyclocheilus anshuiensis [cavefish]) have shown that the expression of oca2, Tyr, Tyrp1 (tyrosinase-related protein1), and Mpv17 (mitochondrial inner membrane protein 17) was lower in the skin of Sinocyclocheilus cavefish than that in the surface fish (Yang et al., 2016). However, considering the high species diversity and habitat differences, it is necessary to study more Sinocyclocheilus fish species to understand additional evolutionary mechanisms underlying pigment loss.
In this study, we conducted comparative histology and Illumina sequencing technology, (c) select pigmentation-related differentially expressed genes (DEGs) by comparing four cave-surface fish groups based on skin transcriptome profiles, and, finally, (d) to identify and validate key candidate genes linked to the differential skin pigmentation between cave and surface fish species through pathway and function annotation of DEGs and qPCR. Our results are expected to provide novel insights into the transcriptional regulation of pigment-related genes in adaptive evolution, leading to a better understanding of the molecular regulation mechanisms of troglomorphic traits in cavefishes.

| Sample collection and photograph with a digital camera
The eight Sinocyclocheilus species were captured in Yunnan Province and Guangxi Zhuang Autonomous Region, China, with three biological replicates for each species (Table 1). Pictures were taken of fresh fish immediately after capture. After general anesthesia with 30 mg/L of MS-222 anesthetic (3-aminobenzoic acid ethyl ester methanesulfonate; Sigma-Aldrich), all lateral skin tissues were surgically excised, divided into two sections, collected into two sterile tubes for comparative histological and transcriptomics (tissue immersed into RNAlater) analysis, and placed in liquid nitrogen. After euthanizing the fish, other tissues (muscles) and organs (heart, liver, brain, etc.) were collected for other studies. The samples were stored at −80°C in the laboratory. Since
Sections were imaged using an optical microscope (OLYMPUS BX 51; Olympus).

| RNA quantification and quality
RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN), and RNA concentration was measured using a Qubit ® RNA Assay Kit in a Qubit ® 2.0 Fluorometer (Life Technologies). RNA samples were run on 1% agarose gels, and their integrity was checked using an RNA Nano 6000 Assay Kit and the Agilent Bioanalyzer 2100 system (Agilent Technologies).

| Clustering and sequencing
Clone clusters were generated on Illumina cBot Cluster Generation System, using TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer's instructions. After cluster generation, highthroughput sequencing of library preparations was performed on Illumina Hiseq 2000, and paired-end reads were generated.

| Data analysis
The raw reads were cleaned by removing reads containing adapter, reads containing poly-N, and low-quality reads through in-house Perl scripts. Furthermore, to ensure the high quality of the data used for downstream analyses, the quality (Q20, Q30, GC-content, and sequence duplication level of the clean data) of the clean data was determined. Clean and high-quality transcriptome from the eight fish species were assembled using Trinity (min_kmer_cov: 2; other parameters: default values) (Grabherr et al., 2011). The whole sequence data were submitted to the DRYAD database (https://datad r yad.org/stash/ share/ t5cZI XoVUg yhpzE P6z-GN6xj c5EU3 TvPPE wbdIo 7siI).

| Gene function annotation
Gene functions of unigenes were annotated based on the homology searches of the following major public databases: NCBI nonredun-

| Differential expression analysis
According to different pairing modes, we investigated more than 30 combinations of pair-group among the eight species, and four species pairs (M vs. K, Q vs. X, ji vs. T, J vs. D; Table 1) were randomly selected for further description. The target genes were identified according to the following three criteria: (a) genes with expression log 2 |Foldchange| ≥ 2 and false discovery rates (FDR) ≤0.001 were considered to be differentially expressed. (b) The functions and pathways of these DEGs were involved in the formation of fish skin pigmentation. (c) The regulation trends (upregulation or downregulation) of these DEGs must be consistent across at least three of the four groups. Differential expression analysis of all four groups was performed using the DESeq R package (1.10.1) (Love et al., 2014). DESeq analyses count data from high-throughput sequencing to determine differential expression using a model based on a negative binomial distribution. The resulting P-values were adjusted using the Benjamini-Hochberg procedure to factor in FDR. DEGs were determined by comparing levels of gene expression within each group after using Bowtie and RSEM to calculate fragments per kilobase of exon per million mapped reads (FPKM) (Langmead et al., 2009;Li & Dewey, 2011). Genes with expression log 2 |Foldchange| ≥ 2 and FDR ≤ 0.001 were considered to be differentially expressed.

| Validation of DEGs by qPCR
RNA-seq results of a total of 9 DEGs were validated by qPCR; β-actin (housekeeping gene) was used as an internal reference. qPCR analysis using the same RNA samples as for the transcriptome profiling was performed using qTOWER2.2 (Analytik Jena, Jena, Germany); the primer sequences are listed in (Table S1). The amplification conditions were as follows: initial denaturation at 95°C for 3 min, followed by 39 cycles of denaturation at 95°C for 10 s and extension at 58°C for 30 s. Finally, melting curves were generating from 60 to 95°C. All reactions were performed with three technical and three biological replicates. Relative gene expression was calculated with the 2 −ΔΔCT method using qPCRsoft 3.2 software (The value of threshold cycle Ct was used to calculate the DEG expression fold change. ΔCt = Ct target genes − Ct β-actin , ΔΔCt = ΔCt control − ΔCt Indicated condition, log 2 (FC) = 2 −(ΔΔCt) .

| Weighted gene coexpression network analysis (WGCNA)
We used the R package WGCNA for weighted correlation network analysis. By using gene expression data from the skin tissue of cave and surface fish, we constructed a coexpression network to find additional important genes associated with each skin color phenotype. First, we constructed a gene-gene similarity network (Pearson's correlation) for all unigenes. In this analysis, we filtered unigenes with expression quantity <1 to calculate the optimal power value, and transformed these values into an adjacency matrix under the soft power (beta = 7). Next, all unigenes were hierarchically clustered, and the network was divided into modules.
Gene modules corresponding to the branches cutoff of the gene tree were color-coded (networkType: signed; soft power: 7; min-ModuleSize: 50; minKMEtoStay: 0.3; mergeCutHeight: 0.20). The "gray" module contained unigenes that could not be associated with any expression patterns. To find the core genes (hub genes) that are significantly associated with the phenotypic trait, we focused on the most highly correlated modules (R 2 >|.5|, p < .005) and the top three most correlated genes were selected for further analysis.

| Color observations in Sinocyclocheilus skin
The eight Sinocyclocheilus species were photographed in the living state with a digital camera (Figure 1), and skin color was noted.
Based on the body color, the four surface fish species were divided into three different subtypes within the same type of habitat. Sinocyclocheilus maculatus skin has a large number of black flecks and appears golden yellow, whereas the skins of S. qiubeiensis and S. oxycephalus exhibit yellow coloration but with fewer flecks.
Distinct from the other three surface fish species, S. jii has a charcoal gray body color with little black spots. The four cave-dwelling Sinocyclocheilus species were divided into two color subtypes:

| Unigene functional annotation
A total of 131,766 unigenes from the eight species of Sinocyclocheilus were annotated using several databases (including Nr, GO, COG, F I G U R E 1 Digital photographs of eight Sinocyclocheilus species in living state and KOG, Table S4). The BLASTx similarity analysis of the unigenes against the NCBI Nr database found matches for most of the uni-

| DEG analysis
DEGs were identified in the four fish groups (M vs. K, Q vs. X, ji vs. T, J vs. D). 4,332 upregulated and 6,297 downregulated) followed by Q vs.
X (9,635 DEGs: 3,628 upregulated and 6,007 downregulated) ( Figure 4b). In every surface fish and cavefish group compared, the proportion of downregulated genes was higher than that of upregulated genes (Table S5)

| DEG functional and pathway analysis in the databases
GO, KEGG, and COG annotations were performed to identify the functions and biochemical pathways of DEGs to further filter the candidate genes that may be responsible for differences in skin pigment (Table S5). DEGs in all four groups were classified into three functional categories (biological process, cellular component, and molecular function) in the GO database. Most of the DEGs of the four groups were annotated to seven GO subcategories, namely "cellular process," "single-organism process," "metabolic process,"  (Table S6).
To further study the function and BPs of DEGs, and identify target genes annotation of DEGs (four groups) was performed using KEGG pathways ( Figure S5). Under the category of the biochemical pathway in the KEGG classification, DEGs in all four comparison groups were associated with cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. Furthermore, the KEGG subcategories with a high proportion of DEGs in the four groups were endocytosis, protein processing in the endoplasmic reticulum, ribosome, and MAPK signaling pathway. Interestingly, some pathways showed differences among the four groups, such as the cardiac muscle contraction pathway, which was annotated in all groups except for the Q vs. X group, and the RIG-l-like receptor signaling pathway and the fatty acid metabolism, which were only observed in the ji vs. T and M vs. K groups, respectively. In this study, DEGs in all four groups were partially annotated to some pigment-related pathways such as the melanogenesis, MAPK signaling, and Wnt signaling pathways.

| Target DEGs selection and expression trend validation
After the functional analysis and pathway annotation of DEGs in the four groups, we selected target genes according to the restrictions stated above. Five DEGs were found in the melanogenesis  Note: "Inf" and "−Inf" indicate that the gene is not expressed in the cavefish or surface fish. "<|2.0|" indicate that the log 2 |FC| <2 (no significant difference) pathway, and the expression of these genes in all surface-cave pairs was checked. Furthermore, seven DEGs involved in Wnt signaling, MAPK signaling, and apoptosis pathways, which are putatively involved in skin pigmentation, were selected from a total of 35,101
The expression of a total of 9 selected DEGs was verified by qPCR to confirm the accuracy of the RNA-Seq data ( Figure 5).

| Weighted gene coexpression network analysis (WGCNA)
After filtering of genes with low levels of expression, the modules were divided based on clustering genes with similar expression patterns and calculating the correlation between the eigenvalue and traits of each module. The coexpression network was used to link gene expression to five skin color phenotypes quantitatively. As a result, all analyzed genes resulted in a total of 35 color modules. We found that the color subtypes of each cavefish and surface fish have different highly correlated modules ( Figure 6). Here, we focused on the five most correlated modules: "steel blue," "dark green," "grey60," "green-yellow," and "black" (R 2 > |.5|; p < .005) to determine hub genes. The top three hub genes based on eigengene connectivity in each module may play an important role in each trait. The details and functions of these genes are presented in Table 3.

| D ISCUSS I ON
For cave-dwelling fish species, the loss of evolutionarily undesirable traits is pivotal to control energy cost and may be one of the reasons for the convergent evolution of regressive traits; however, the regulatory mechanism underlying this process is complex (Protas et al., 2007). Studies on skin pigments in A. mexicanus have reported that different genes determine the degrees of skin pigment loss (brown and albino) (Gross et al., 2009;Protas et al., 2006). Besides, the mutations in these genes may be pleiotropic. Helena et al. found that the albinism of A. mexicanus caused by loss-of-function mutations in the oca2 gene is a by-product of catecholamine-related behaviors (anesthesia resistance), which is related to the capacity to be responsive to stimuli in the caves (reduced sleep) (Helena et al., 2018), which indicates that the cause of reduced pigmentation of cavefish may be more complicated than earlier believed, and there may be additional, yet unidentified pigment-related genes.
The link between loss of pigmentation and the underlying genetic mechanisms in cave-dwelling fishes can be better understood by investigating the pattern of transcriptional regulation of genes that may cause loss of skin pigmentation in Sinocyclocheilus cavefish.
In this study, we selected eight representative (skin color and habitat type) Sinocyclocheilus species, which included four surface species (with functional eyes and pigmentation) and four cave species (with reduced or absent eyes and pigmentation) to evaluate the phenotypic differences in skin pigmentation between the two ecological types of Sinocyclocheilus species. We aimed at identifying and confirming key candidate genes linked to the diverse skin pigmentation between cave and surface fish species through comparative transcriptomics analysis.

| Relationship between skin color and habitat of Sinocyclocheilus species
Skin color can be influenced by environmental interactions (Leclercq et al., 2010). During the long-term sampling process and looking at lots of other of species within this genus, we found that the skin color of Sinocyclocheilus species varies widely among species of the surface also cave-dwelling fish, but remains consistent within populations. The eight representative Sinocyclocheilus species selected in this study covered almost all the color types we have observed in the wild. Our comparative histological analysis of the two ecotypes of Sinocyclocheilus species clearly indicated the differences in distribution patterns, quantity, and melanosome intracellular arrays of melanocytes in the skin between cavefish species and surface fish species, suggesting that different habitats may be one of the crucial reasons for this difference. Melanocytes synthesize melanin to protect skin by absorbing ultraviolet radiation (UV) (Kim et al., 2018).
Previous studies on fish, such as whitefish larvae (Coregonus lavaretus) (Winberg, 2000), red sea bream (Pagrosomus major) (Kumai, 2005), and rainbow trout (O. mykiss; Little, 1995) showed that the concentration of skin melanin and eumelanin precursor increase in response to UV exposure to adapt to the changing environment. Melanin can absorb specific wavelengths of light to protect the skin against UV radiation, resulting in the dark gray skin of zebrafish (Kelsh, 2004).
However, troglobites live in completely enclosed underground habitats characterized by permanent darkness, and complete lack of autochthonous photosynthesis, leading to limited food resources (Bussotti et al., 2018). When adapting to these extreme environments, the accumulation of melanin in cavefish skin is a waste of limited resources (Culver & Pipan, 2009). Interestingly, among the four cavefish species, digital photography showed that the skin color of S. rhinocerous was the darkest but still transparent, and the stereoscopic microscopic photographs and paraffin section also confirmed that the number of melanocytes in the S. rhinocerous was greater than that in the other cavefish species. Li et al. (2000) studied the cave habitats at sampling sites of S. rhinocerous and found that this F I G U R E 6 Coexpression modules conducted by WGCNA. (a) hierarchical cluster tree constructed by WGCNA shows coexpression modules (the gray module represents genes that are not assigned to specific modules). Each branch in the tree point connects to a gene. Genes were assigned to network modules by dynamic tree cut after the dynamic tree cut algorithm has been used to identify all modules. minModuleSize:30 (lowest number of genes in each module) and Merged dynamics (the modules with similar expression patterns (80%) were then merged.); (b) Trait-module associated heat map. Column: traits; row: modules. The number in each grid represents the correlations between the module and corresponding traits (correlation coefficients) and p-value (in parentheses).
Color of the grid indicates the correlation (the deeper the color, the stronger the correlation) species occupies funnel-shaped doline caves, which partially receive solar radiation; in contrast, other cavefishes, such as Sinocyclocheilus altishoulderus, live in completely enclosed underground rivers. This half-open environment may partially explain its unique skin color.
Based on comparative histology, we divided the Sinocyclocheilus cavefish and surface fish body colors into two and three subtypes, respectively. The coexpression network of WGCNA was used to link gene expression with five skin color phenotypes quantitatively. The "steel blue" module and the "dark green" module were positively correlated with the cavefish gray and pink color traits, respectively, indicating that the two modules might play an important role in cavefish pigmentation. We found hub genes among the top three hub genes, C1q-like (complement C1q) and FN1 (fibronectin 1), in the "gray-steel blue" module whose expression was related to solar ultraviolet radiation (UVR). Mei et al. found that C1q-like is involved in UV-induced apoptosis in zebrafish. After the UV exposure, the transcripts of C1q-like were upregulated 3-4 fold (Mei et al., 2008). Fibronectin is a globular glycoprotein ubiquitous in the dermal extracellular matrix (ECM), and exposure to UV irradiation modulates FN1 expression, thereby enhancing ECM degradation (Hibbert et al., 2017;Parkinson et al., 2015). Furthermore, in the "pink-dark green" module, PLEC (plectin) is considered the top hub gene. Previous studies of human skin color have shown that the attenuated expression of PLEC leads to increased melanosome uptake by keratinocytes and skin hyperpigmentation after UVA exposure (Coelho et al., 2015). This finding may reconfirm that the differences in body color between cavefish in this study are related to the UV intensity among different habitats.
Wild Sinocyclocheilus species establish complex interactions with their habitats, which may require the development of different mechanisms or regulatory genes for each skin color subtype of cavefish such as in S. rhinocerous and S. altishoulderus. Furthermore, some surface fish or cavefish may occasionally change their habitat through the underground river for predation (Chen et al., 2019), which may be part of the reason for the diversity of body colors between the same type of fish. Complex wild habitats may have diverse influencing factors such as seasonal changes and diet, which may affect the formation of host skin pigmentation. Nevertheless, further research is needed to reveal the exact mechanisms underlying these variations.

| DEGs related to loss of pigmentation of Sinocyclocheilus fish
The distribution of the genus Sinocyclocheilus is very narrow, and the habitat of each species is merely a single cave or one surface waterbody (Lunghi et al., 2019). Furthermore, a large number of variables (age, geographical location, and species) may affect gene transcription. As such, in the grouping, these surface-cave fish pairs were combined only based on the "habitat type" to find the DEGs of each group, before comparisons were made. We believe that the candidate genes selected in this way could truly reflect the relationship between skin pigmentation and the two habitat types (cave and surface waterbody), rather than the effect of other variables.  (Jiang et al., 2014) and stockpiled by the melanosomes in melanophores (Slominski et al., 2004). Previous studies have shown that the expression of melanin is regulated by multiple genes in the melanogenesis pathway, such as Mc1r, MITF, GNAQ (α melanocyte-stimulating hormone).
Accordingly, in this study, we first focused on the melanogenesis pathway for candidate DEGs screening and found five genes: GNAQ, AC (adenylate cyclase), NRAS (NRAS proto-oncogene, GTPase), SFRP2 (secreted frizzled-related protein 2), and PKA (protein kinase) with significant differential expression, and further verified the expression trends of these five genes in all surface-cave pairs.
In the melanogenesis pathway, Mc1r on the cell membrane is activated by α-MSH, which results in the stimulation of Gnαq, thereby activating AC to increase the production of cAMP and PKA activation. PKA upregulates the expression of MITF, leading to a rise in melanin synthesis (Buscà & Ballotti, 2000;García-Borrón et al., 2005;Oscar & Van, 2016;Van Raamsdonk et al., 2009), eventually leading to darker skin coloration through the intracellular dispersal of membrane-bound pigment granules (melanosomes) within the melanophore (García-Borrón et al., 2005;Hoekstra, 2006;Lin & Fisher, 2007;Yamaguchi et al., 2007). Notably, we found three genes, GNAQ, AC, and PKA, that showed a significant downregulation trend in most cavefish skins.
GNAQ encodes Gnαq, a subunit of a trimeric G protein complex that binds to the endothelin B receptor in melanocytes (Raamsdonk et al., 2004), and plays an important role in the regulation of pigment formation ( Van Raamsdonk et al., 2009). We found that GNAQ was significantly downregulated in more than 75% of surface-cave pairs (12/16 pairs). Recent research found that the expression of GNAQ in black mouse skin was significantly higher than that in the white mouse skin (Yin et al., 2015). In fish species, GNAQ showed a significant association with skin pigmentation in three spine sticklebacks (Gasterosteus aculeatus; Greenwood et al., 2011). Furthermore, Gessi et al. (2013) found that a point mutation in GNAQ may be related to the development of primary melanocytic tumors in humans.
They identified six cases harboring mutations in codon 209 of the GNAQ gene. As members of a cascade of regulatory genes in the melanogenesis pathway, GNAQ, AC, and PKA control the process of melanin production by regulating the expression of upstream or downstream genes (Bennett & Lamoreux, 2003). This suggests that the above-mentioned genes, which are downregulated in cavefish species, are likely to be closely related to the loss of pigmentation in Sinocyclocheilus cavefish.
In addition to the melanogenesis pathway, genes (Tyr, TYRP1, P38, etc.) in Wnt signaling pathway (Fujimura et al., 2009), MAPK signaling pathway and apoptosis pathway (Squarzoni et al., 2011), have been indicated to be involved in the regulation of pigmentation. Tyrosinase, encoded by the Tyr gene, is the rate-limiting enzyme in melanogenesis and is an important regulatory factor in the synthesis of melanin. Previous studies reported that loss or downregulation of Tyr gene leads to an albino phenotype in zebrafish (Page-McCaw et al., 2004), duck (Li et al., 2012), mouse (Ito & Wakamatsu, 2011) and rabbit (Song et al., 2017 transcriptomics. We infer that these genes may participate in the partial regulation of skin pigmentation, and the downregulation of these genes may lead to the pigmentation loss in Sinocyclocheilus cavefish species. Furthermore, we found that the trend of gene expression was different in each group (Table 2), and some genes even no significant difference in some group may due to different mechanisms of pigment loss, which could be responsible for these different species. Therefore, future studies should include gene knockout or overexpression studies to validate the function of these genes.
In this study, we found that the results of qPCR of some pairs may not perfectly support RNA-seq but they used the same RNA. In order to minimize RNA degradation, in this study, we sent the skin tissue to the sequencing company directly to extract the RNA, and then sent the RNA back to our laboratory to conduct qPCR, we think this may cause the degradation of a small amount of RNA during the sample delivery. Furthermore, since the two experiments are not conduct at the same laboratory, the final results may have some differences.
In conclusion, to further elucidate the underlying genetic mechanism of skin pigment loss in Sinocyclocheilus cavefish species, we conducted a comparative histological, WGCNA, and comparative transcriptomics analyses of the skin of Sinocyclocheilus fishes dwelling on surface and in caves, and found differences in the distribution, quantity, and morphology of skin melanocytes. A total of 35,101 DEGs were found in four surface-cave groups. Through GO, COG, and KEGG pathway annotations, we identified 11 candidate genes which showed significant differential expression in Sinocyclocheilus cave and surface fish species analyzed in this study, that may participate in the regulation of skin pigmentation. Most of the DEGs were downregulated in Sinocyclocheilus cavefishes as validated by qPCR. However, these genes require further functional validation in Sinocyclocheilus. This study provides a strong foundation for better understanding of the molecular genetic mechanisms underlying troglomorphic traits in cavefish species.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

E TH I C S A PPROVA L A N D CO N S E NT TO PA RTI CI PATE
The Sinocyclocheilus fish species used in this study were caught from wild water bodies, and no specific permissions were required. All experiments were conducted after review and approval from the local Ethical Committee at Yunnan University in accordance with China's local and global ethical policies (Grant No: ynucae 20190056), and all the procedures were approved and assisted by the local government.

CO N S E NT FO R PU B LI C ATI O N
Not applicable.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data generated or analyzed during this study are included in this published article and its supplementary information files. The raw reads produced in this study were deposited in the DRYAD database (https://datad ryad.org/stash/ share/ t5cZI XoVUg yhpzE P6z-GN6xj c5EU3 TvPPE wbdIo 7siI).
Differential expression of Gnαq gene in white and black skin tissues of mice. Journal of Zhejiang University (Agriculture & Life Sciences), 41(4), 414-420.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.