Scale Development-Related Genes Identified by Transcriptome Analysis

Scales, as key structures of fish skin, play an important role in physiological function. The study of fish scale development mechanisms provides a basis for exploring the molecular-level developmental differences between scaled and non-scaled fishes. In this study, alizarin red staining was used to divide the different stages of zebrafish (Danio rerio) scale development. Four developmental stages, namely stage I (~17 dpf, scales have not started to grow), stage II (~33 dpf, the point at which scales start to grow), stage III (~41 dpf, the period in which the scales almost cover the whole body), and stage IV (~3 mpf, scales cover the whole body), were determined and used for subsequent transcriptome analysis. WGCNA (weighted correlation network analysis) and DEG (differentially expressed gene) analysis were used for screening the key genes. Based on the comparison between stage II and stage I, 54 hub-genes were identified by WGCNA analysis. Key genes including the Scpp family (Scpp7, Scpp6, Scpp5, and Scpp8), the Fgf family (Fgfr1b and Fgfr3), Tcf7, Wnt10b, Runx2b, and Il2rb were identified by DEG analysis, which indicated that these genes played important roles in the key nodes of scale development signal pathways. Combined with this analysis, the TGF-β, Wnt/β-catenin, and FGF signaling pathways were suggested to be the most important signal pathways for scales starting to grow. This study laid a foundation for exploring the scale development mechanism of other fishes. The scale development candidate genes identified in the current study will facilitate functional gene identifications in the future.


Introduction
Scales are the structures that cover the surface of birds' legs, reptiles' bodies, and fishes' skin. Different from birds and reptiles, fish scales are derived from the dermis. Fish scales, like fin rays, belong to the fish exoskeleton and generally refer to all the hard, flat bones in the fish skin [1,2], which protect the fish body. Fish scales account for 2% to 3% of the body weight and play important roles in protecting fish from external microorganisms, reducing frictional damage caused by particles in water, and maintaining fish body shape [3,4].
Studies of fish scales mainly focus on morphological comparison and observation, classification, structural composition, and the genes involved in scale development [5][6][7]. Recently, researchers have come to believe that the starting position of scale development is mainly affected by epigenetic factors [8]. At the gene level, the related signal pathways and genes mainly include the EDA/EDAR, Shh, FGFs, BMP, and Wnt/β-catenin signal pathways. [9][10][11]. For example, many studies showed that the Wnt signaling pathway plays a key role in the development of fish scales [12]. However, fish scale development is a complex process regulated by a variety of signals and there is no clear regulatory network reported.
Zebrafish (Danio rerio), a small tropical freshwater fish, has the advantages of short maturation cycle and easy reproduction, which makes zebrafish an important experimental model for the study of fish scale development [13,14]. Although no major breakthroughs have been made regarding the molecular mechanism, research on zebrafish scale development began decades ago. Several gene mutants (such as Fls, nkt, spd, and rs-3) resulted in fish scale development defects in zebrafish [15]. These mutants provide a convenient and relevant basis for fish scale development studies.
With the rapid development of high-throughput sequencing technology, the process of obtaining gene-level information has become possible [16]. Transcriptome analysis is a useful technology and has been widely used in developmental biology studies. Many studies regarding skin formation patterns and epidermal structural development have been carried out on mammals and birds through transcriptome analysis [17]. However, no transcriptome analysis regarding fish scale development has been reported.
In this study, the developmental process and signal transduction pathways of zebrafish scale development were analyzed by transcriptome sequencing. For key genes' screening, WGCNA (weighted correlation network analysis) analysis was carried out. WGCNA is an R package for weighted gene co-expression network analysis, which can simplify a large amount of complex gene data and quickly identify core genes related to target traits, and has been widely used in animal and plant research [18]. The purpose of this study was to determine the key genes and pathways in the development of fish scales. These results will increase our understanding of the molecular mechanism of fish scale development.

Zebrafish Maintenance
Zebrafish (AB type) were maintained at 28 • C in a closed water ultrafiltration purification system at our zebrafish feeding laboratory (Shanghai Haisheng Biological Experimental equipment Co. Ltd., Shanghai, China). The animals used in the present study were cultured and euthanized following the terms of use of animals approved by the Institutional Animal Care and Use Committee at the Shanghai Ocean University (Shanghai, China; SHOU-DW-20171022). Zebrafish embryonic developmental stages were determined in days post-fertilization (dpf) or months post-fertilization (mpf). According to the days after embryo fertilization, after staining observation, we identified four stages of scale development, namely stage I (~17 dpf, scales have not started to grow), stage II (~33 dpf, the point at which scales start to grow), stage III (~41 dpf, the period in which the scales almost cover the whole body), and stage IV (~3 mpf, the period in which the scales completely cover the fish body).

Alizarin Red Staining and Observation
Zebrafish tissue samples were fixed with 4% paraformaldehyde (Sangon, Shanghai, China) overnight and rinsed with ddH 2 O for 20 min the next day. Equal volumes of 1% KOH and 1% H 2 O 2 (Sangon, Shanghai) were mixed to make a bleaching solution and fish samples were immersed and exposed to the air for at least 10 min. After the bleaching solution was removed, protease solution (100 mL system: 65 mL ddH 2 O, 35 mL saturated sodium borate, and 1 g trypsin (Sangon, Shanghai)) was then used to remove the excess muscle tissue until the fish body became transparent. After digestion, the protease digestion solution was discarded and the fish samples were rinsed with 1% KOH several times; then soaked in 1% KOH solution; and 1% alizarin red solution (1 g alizarin red (Yuanye, Shanghai, China) dissolved in 100 mL 0.5% KOH) was gradually added. When the fish body surface turned light purple, the fish bodies were observed under a microscope (Nikon SMZ1500, Tokyo, Japan).

RNA Isolation and Sequencing
Total RNA was extracted from zebrafish skin (with scales) collected at four developmental stages with the Trizol reagent (Ambion, Austin, TX, USA). The sampling site was at the caudal peduncle shown in the alizarin red staining diagram ( Figure 1). Only high-quality samples (OD260/280 ≥ 1.8, OD260/230 ≥ 1.8) were used to build a sequencing library.

RNA Isolation and Sequencing
Total RNA was extracted from zebrafish skin (with scales) collected at four developmental stages with the Trizol reagent (Ambion, Austin, TX, USA). The sampling site was at the caudal peduncle shown in the alizarin red staining diagram ( Figure 1). Only high-quality samples (OD260/280 ≥ 1.8, OD260/230 ≥ 1.8) were used to build a sequencing library. The library was established according to the instructions of the VAHTS Stranded mRNA-seq Library Prep Kit for Illumina V2 (Vazyme, Nanjing, China). MRNA with a poly-A tail was enriched using magnetic beads with oligo dT. The combined mRNA was eluted and broken into fragments with buffer and the first chain of cDNA was formed by reverse transcription with random primers. RNase H digestion solution was used to digest the mRNA strand in the hybrid double-strand, synthesize the complementary strand of cDNA, realize the synthesis of the second strand of cDNA, and conduct end repair and dA-tailing. After the linker was connected, the second chain containing U was digested with the UDG enzyme and then PCR primer mix 3 in the kit was used to amplify and enrich the library. The sequencing work was carried out by Meiji Biological Company (Shanghai, China) using an illumina 4000.

Weighted Gene Co-Expression Analysis (WGCNA) of the Sequencing Data
WGCNA is a freely accessible R package for constructing the weighted gene coexpression network [19]. After cluster analysis, two outlier samples were deleted and the soft threshold was determined to be 12. On this basis, the co-expression network was constructed and the relationship between sample characteristics and modules was analyzed. Finally, it was determined that the turquoise module (turquoise, cor = 0.74, p < 0.01) was the module with the highest correlation with scale development. Then, correlation analysis between the module and gene importance was conducted. Finally, Cytoscape (Version 3.6.1) was used to draw the gene network diagram and Cytohubba's EPC algorithm was used to identify hub-genes [20]. The library was established according to the instructions of the VAHTS Stranded mRNA-seq Library Prep Kit for Illumina V2 (Vazyme, Nanjing, China). MRNA with a poly-A tail was enriched using magnetic beads with oligo dT. The combined mRNA was eluted and broken into fragments with buffer and the first chain of cDNA was formed by reverse transcription with random primers. RNase H digestion solution was used to digest the mRNA strand in the hybrid double-strand, synthesize the complementary strand of cDNA, realize the synthesis of the second strand of cDNA, and conduct end repair and dA-tailing. After the linker was connected, the second chain containing U was digested with the UDG enzyme and then PCR primer mix 3 in the kit was used to amplify and enrich the library. The sequencing work was carried out by Meiji Biological Company (Shanghai, China) using an illumina 4000.

Weighted Gene Co-Expression Analysis (WGCNA) of the Sequencing Data
WGCNA is a freely accessible R package for constructing the weighted gene coexpression network [19]. After cluster analysis, two outlier samples were deleted and the soft threshold was determined to be 12. On this basis, the co-expression network was constructed and the relationship between sample characteristics and modules was analyzed. Finally, it was determined that the turquoise module (turquoise, cor = 0.74, p < 0.01) was the module with the highest correlation with scale development. Then, correlation analysis between the module and gene importance was conducted. Finally, Cytoscape (Version 3.6.1) was used to draw the gene network diagram and Cytohubba's EPC algorithm was used to identify hub-genes [20].

Differentially Expressed Gene Analysis and Functional Enrichment
In order to study the molecular mechanism of zebrafish scale development, we analyzed the number and biological function of DEGs in six groups. Firstly, we mapped all clean reads onto the reference gene sequence using HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts 2) [21]. Significant DEGs were detected by EdgeR software (Version 3.26.8) and both |log2 (fold change)| ≥ 1 and p ≤ 0.05 were used as the filtering thresholds [22]. The gene expression level of each sample was then calculated to determine the fragments per kilobase of the exon model per million mapped fragments (FPKM) using cufflinks software [23] with the default settings.
In addition, to further understand the function of DEGs at different scale developmental stages, we selected genes with differential expression |log2 (fold change) | ≥ 1 and p ≤ 0.05, as well as used the clusterprofiler package in R software for enrichment analysis, and GO terms and KEGG pathways with p ≤ 0.05 were regarded as significant [24].

Quantitative Real-Time PCR Validation of RNA-Seq Data
To verify the accuracy of the transcriptome sequencing results, 11 genes were selected to perform qRT-PCR analysis. First, we extracted RNA from the caudal peduncle of zebrafish at four different stages. The first strand of cDNA was synthesized according to the instructions of the Goldenstar RT6 cDNA Synthesis Kit, version 2 (Tsingke Biological Technology Co., Ltd., Wuhan, China). The specific process was as follows: First, we removed the possible residual DNA in a 10 µL reaction system: 1 µg total RNA, 1 µL gRNA remover, and 1 µL gRNA remover buffer, and added RNase-free water to 10 µL. We incubated at 42 • C for 2 min and 60 • C for 5 min. Then, we added the following ingredients: 4 µL 5 × Goldenstar reaction buffer, 1 µL Goldenstar RT6 (200 U/µL), 1 µL oligo (DT) 17, 1 µL dNTP mix (10 mM), 1 µL DTT, and 2 µL RNase-free water. We incubated at 55 • C for 40 min and 85 • C for 5 min. Then, the synthesized cDNA was used as a template and the primers used are as shown in Table 1. The total qPCR mixture reaction volume of 20 µL contained 10 µL 2 × T5 Fast qPCR Mix (SYBR Green I (Tsingke Biological Technology Co., Ltd., Wuhan, China)), 0.8 µL of each primer (10 µM), 2 µL cDNA template, and 6.4 µL RNase-free water. The quantitative analysis of genes was carried out by using a CFX96 fluorescence quantitative PCR instrument (BioRad, Hercules, CA, USA). The specific qPCR procedure was as follows: pre-denaturation at 95 • C for 1 min, amplification at 95 • C for 10 s, amplification at 58 • C for 10 s, and amplification at 72 • C for 15 s, for 40 cycles in total. The relative expression of each of the ten genes was analyzed using the comparative cycle threshold (2 −∆∆CT ) method (∆CT = CT target gene − CT reference gene , ∆∆CT = ∆CT treatment − ∆CT control ) [25]. Statistical analysis and Pearson correlation analysis were performed with SPSS (Version 25.0, IBM, Armonk, NY, USA). The significance of the test results was tested by t-test. Pearson correlation analyses were used to assess the strength of the relationship between RNA-seq data and RT-qPCR measurements.

Scale Development in Zebrafish
We divided the zebrafish scale development process into four stages by alizarin red staining. Stage I (~17 dpf) is the stage when zebrafish do not grow scales (see Figure 1A,B). Stage II (~33 dpf) is the point at which scales start to grow ( Figure 1C). The third stage (~41 dpf) is when scales extend forward along the midline with the dorsal axis and also vertically up and down until it covers the whole fish body ( Figure 1D). Stage IV (~3 mpf) is the period in which the scales completely cover the fish body.

Gene Expression Quantification and Analysis of DEGs
Gene expression levels for the four fish scale developmental stages are shown and visualized by Venn diagrams (Figure 2A) clearly showing the numbers of genes that were specifically expressed at each stage and the number of genes that were co-expressed between them. The results showed that 20,345, 16,608, 14,577, and 7674 genes were expressed in the four stages, of which 7533 genes were common among the four stages. It can be seen from Figure 2B that the data distribution of each group was relatively scattered; the minimum value is zero and the median is close to upper quartile.

WGCNA Analysis and Hub-Genes Screening
In total, FPKM values of the 35,116 genes for the four developmental stages were obtained. Among them, 7223 genes with FPKM ≥ 1 were selected for cluster analysis by the Fastcluster in the WGCNA software package. In construction of the WGCNA network, it is important to determine the soft threshold. We analyzed the network topology with a threshold power of 1 to 20 and determined both the scale independence and mean connectivity as the relative balance of WGCNA. As shown in Figure 3A, the optimal threshold β value was 12, indicating that after reaching the high threshold of 12, the curve began to flatten and the network topology was almost connected ( Figure 3A).
The "one-step method" was used to construct the network. The difference coefficient between genes was calculated to obtain the genetic system cluster tree. After the gene modules were determined, the eigenvalues of each module were calculated in turn. Cluster analysis was performed on the modules. In total, 7233 genes were divided into 13 modules, among which 10 genes did not belong to any module ( Figure 3B).
In order to explore genes that were highly correlated with the module and phenotype, we conducted a heat map of module-associated phenotypes ( Figure 3C). Among the four scale developmental stages, we mainly focused on stage II, the starting point of scale growth (~33 dpf, zf33). It was clearly shown that compared with other modules, the turquoise module was highly correlated with stage II (zf33), which indicated that genes in the turquoise module might play an important role in zebrafish scale development ( Figure 3C).

Gene Expression Quantification and Analysis of DEGs
Gene expression levels for the four fish scale developmental stages are shown and visualized by Venn diagrams (Figure 2A) clearly showing the numbers of genes that were specifically expressed at each stage and the number of genes that were co-expressed between them. The results showed that 20,345, 16,608, 14,577, and 7674 genes were expressed in the four stages, of which 7533 genes were common among the four stages. It can be seen from Figure 2B that the data distribution of each group was relatively scattered; the minimum value is zero and the median is close to upper quartile.      Then, the connectivity of genes in the turquoise module was calculated. Genes with gene significance > 0.87 and intra-module connectivity > 0.8 were identified as the key genes (hub-genes). In total, 54 hub-genes with high connectivity were identified (see Table S1). In addition, the top 10 key genes were obtained by Cytoscape software(Version 3.6.1, UCSD et al, Califonia, USA)package CytoHubba(Version1.5.1, UCSD et al, Califonia, USA) using the EPC algorithm ( Figure 3D).

GO and KEGG Pathway Enrichment Analysis
In order to explore important pathways and biological processes involved in scale development, the gene expression in the four developmental stages was compared and both GO and KEGG enrichment were carried out. As shown in Figure 4A, compared with stage I, stage II mainly involves metabolic processes, such as ATP metabolism and translation, indicating that a large amount of energy is required in the process of scale formation. In addition, the corresponding enrichment pathways are mainly energy metabolism pathways, such as carbon metabolism, ribosome, and other pathways ( Figure 4B).

Gene-Act Network
Compared with stage I, 320 genes (log2FC ≥ 3.0, p ≤ 0.01) in stage II (the starting point of scale growth) were selected to explore the relationship between DEGs. Among them, 48 DEGs were identified for the linkage map constructing of the gene interaction network by Cytoscape software (Figure 4C). Several gene families, such as Fgf (Fgfr1b and Fgfr3) and Scpp (Scpp5, Scpp7, Scpp6, and Scpp8), were included ( Figure 4C). It was clearly shown that Scpp6, Scpp5, Fgfr1b, Runx2b, Il2rb, and Tcf7 had complex interactive relationships Compared with stage I, stage III is mainly involved in the biological process of energy metabolism and translation ( Figure S1A), and the corresponding KEGG pathway was also Fishes 2022, 7, 64 8 of 13 related to energy metabolism, such as the TCA cycle and glucose metabolism ( Figure S1B). Compared with stage I, stage IV is mainly involved in a series of biological processes of skeletal muscle growth ( Figure S1C) and both focal adhesion and the ECM receptor in the enriched pathway are mainly involved in the process of cell adhesion ( Figure S1D).

Gene-Act Network
Compared with stage I, 320 genes (log2FC ≥ 3.0, p ≤ 0.01) in stage II (the starting point of scale growth) were selected to explore the relationship between DEGs. Among them, 48 DEGs were identified for the linkage map constructing of the gene interaction network by Cytoscape software ( Figure 4C). Several gene families, such as Fgf (Fgfr1b and Fgfr3) and Scpp (Scpp5, Scpp7, Scpp6, and Scpp8), were included ( Figure 4C). It was clearly shown that Scpp6, Scpp5, Fgfr1b, Runx2b, Il2rb, and Tcf7 had complex interactive relationships with other genes, which indicated the important roles that those genes played ( Figure 4C).

Key Genes and Pathways Related to Scale Development
The 10 hub-genes obtained by WGCNA analysis were highly expressed in stage II of scale development ( Figure 5A). Through traditional transcriptome analysis methods, this study found that the genes of these families were also upregulated in stage II of scale development, which included the Scpp family (Scpp6, Scpp7, Scpp5, and Scpp8), Wnt (Wnt10b), the Fgf family (Fgfr1b and Fgfr3), Tcf7, Il2rb, and Timp4.1 ( Figure 5B). In order to further understand the process of scale formation, after further analysis of these genes, we found that they are mainly in three pathways: the Wnt, TGF-β, and FGF signaling pathways ( Figure 5C). For further details about these genes, see Table S2. found that they are mainly in three pathways: the Wnt, TGF-β, and FGF signaling pathways ( Figure 5C). For further details about these genes, see Table S2.

qPCR Analysis Verifications
In order to further evaluate the reliability of transcriptome data, we selected 11 genes for qPCR detection. The correlation between RT-qPCR and RNA-seq results was analyzed by calculating the Pearson correlation coefficient (R 2 ). The expression data of 11 selected genes showed the consistency between RNA-seq and RT-qPCR. The R 2 value (>0.7) of the correlation between RT-qPCR and RNA-seq data confirmed the reliability of differential expression analysis in this study. Although the R 2 values of two genes (Runx2b and Col10a1a) were less than 0.7, their trend of expression was the same (Figure 6).

qPCR Analysis Verifications
In order to further evaluate the reliability of transcriptome data, we selected 11 genes for qPCR detection. The correlation between RT-qPCR and RNA-seq results was analyzed by calculating the Pearson correlation coefficient (R 2 ). The expression data of 11 selected genes showed the consistency between RNA-seq and RT-qPCR. The R 2 value (>0.7) of the correlation between RT-qPCR and RNA-seq data confirmed the reliability of differential expression analysis in this study. Although the R 2 values of two genes (Runx2b and Col10a1a) were less than 0.7, their trend of expression was the same (Figure 6).
In order to further evaluate the reliability of transcriptome data, we selected 11 genes for qPCR detection. The correlation between RT-qPCR and RNA-seq results was analyzed by calculating the Pearson correlation coefficient (R 2 ). The expression data of 11 selected genes showed the consistency between RNA-seq and RT-qPCR. The R 2 value (>0.7) of the correlation between RT-qPCR and RNA-seq data confirmed the reliability of differential expression analysis in this study. Although the R 2 values of two genes (Runx2b and Col10a1a) were less than 0.7, their trend of expression was the same (Figure 6). Figure 6. The comparison of the selected gene expression patterns in RT-qPCR and RNA-seq showed a high correlation between the two methods. The calculation and mapping of the log 2 ratio of expression changes in the other three stages was relative to the first stage (the log 2 ratio of stage I was set to 0). The blue line shows the RT-qPCR results and the orange line shows the RNA sequence results. The Pearson correlation coefficient (R 2 ) was used to measure the correlation between RT-qPCR and RNA-seq results. The results were analyzed by t-test. * indicates significant correlation (p < 0.05) while ** indicates highly significant correlation (p < 0.01).

Discussion
In fishes, scales are produced by directed differentiation from fish skin primitive stem cells [26]. Relevant research data show that the skin has already been formed at the beginning of the scale developmental process [27]. In this study, we used zebrafish as a model to study scale development by using alizarin red staining and we divided the developmental process into four stages: stage I (17 dpf), stage II (~33 dpf), stage III (~41 dpf), and stage IV (~3 mpf) (Figure 1). Studies have shown that first scale appearance in zebrafish is related to fish size. Zebrafish with a body length of 8.1-8.2 mm at 30 dpf had the first scale at the caudal peduncle [1]. The time of stage II of scale development (~33 dpf) was similar to that reported in the literature. The molecular mechanisms of fish scales starting to grow attract the interest of researchers. In order to find the key genes of scale development, we used WGCNA and differentially expressed gene analyses to perform gene identification.
According to the WGCNA analysis results, the study first identified hub-genes related to scale development, which included nusap1, tbx1, smad6a, and mynn ( Figure 3). These genes have been reported to be involved in skin or bone development and have been confirmed to participate in many known skeletal development-related pathways. For example, Tbx1 is connected to several major signaling systems, such as FGF, WNT, and SHH, and is involved in the cell proliferation and regulation of cell shape and cell dynamics. It is expressed only in hair follicles in healthy human skin [28][29][30]. Id1 is a small molecule activator of BMP signaling [31] and is positively correlated with bone morphogenetic protein (BMP) signal metabolism during fin repair [32]. Bmp participates in the TGF-β signaling pathway and is a key regulatory gene in the pathway. Tgif1 is in the TGFβ1/Smad2/3 signaling pathway [33]. All these observations indicate that these hub-genes play both direct and indirect roles at key nodes in the scale development pathway.
Since stage II is the point at which scales start to grow and stage I is the stage zebrafish do not grow scales, the transcriptome comparison between II and I is very important. It is helpful for us to identify important upregulated genes involved in scale development. Therefore, we focused on the gene interaction map of groups (II vs. I) in which we found that the Scpp family (Scpp7, Scpp6, Scpp5, and Scpp8), the Fgf family (Fgfr1b and Fgfr3), Tcf7, Wnt10b, Runx2b, and Il2rb are closely connected, indicating that these genes play a key role in scale development ( Figure 4C). Of course, some of the selected genes are consistent with existing literature reports. For instance, Z Liu et al. [34] further analyzed the transcriptome of scaled and non-scaled fishes, and found that two genes that play a leading role among the upregulated genes are the apolipoprotein and Scpp genes. In addition, Scpp has been reported to be involved in skeletal muscle development and tooth tissue calcification [35][36][37][38][39][40]. Wntl0b participates in the regulation of animal hair follicle growth and regeneration hair follicles are highly expressed during embryonic development, which promotes hair growth and participates in the maintenance of the hair follicle cycle [41]. FGFR (fibroblast growth factor receptor) is a member of the immunoglobulin superfamily and the gene itself also belongs to a polygenic family [42]. The FGFR gene is considered to be an important developmental gene which can regulate the differentiation and proliferation of a variety of cells [43]. Studies have shown that Fgfr1a and Fgfr20a mutations will lead to changes in the scale size of zebrafish during scale development [44,45]. The Fgfr1b and FGFR3 identified in this study belong to the FGFR family and are also highly expressed at the stage of scale formation, which is consistent with the results of existing reports. The TCF/lef transcription factor family can regulate the Wnt signaling pathway, including the Tcf7 gene identified in our study [46,47]. Ducy et al. [48] found that the Runx2 gene regulates bone formation by controlling osteoblast activity.
In the comparison of the RT-qPCR and RNA-seq results of the screened genes, we conducted correlation analysis and the t-test was used to test whether the correlation was significant. The results showed that the R 2 of most genes was greater than 0.7 and two genes (EDAR and Bmp2a) had a very significant correlation (p < 0.01); the correlation of one gene (Scpp7) was significant (p < 0.05) but the correlation of the others was not significant (p > 0.05). However, while the R 2 value of the RT-qPCR and RNA-seq of the genes was less than 0.7 or their correlation was not significant, their expression trend was consistent, indicating the reliability of the transcriptome results. This result is consistent with the results presented in many literature reports, that is, the expression trend of RT-qPCR and RNA-seq is consistent [49][50][51].
Scale development is an extremely important and complex process. Therefore, studying the mechanism of scale development and its related signaling pathways can help understand animal evolution, growth, development, and reproduction. It involves many signaling pathways and genes, such as the SHH, FGFs, BMP, EDA/EDAR, and Wnt/βcatenin signaling pathways, as well as transcription factors including Twist2 [9][10][11]52]. These factors control the occurrence and development of scales. However, in our current study, the Shh, EDA/EDAR, BMP, twist2, and ApoE genes were not identified during scale development, although they are still important. Based on the results of WGCNA analysis and traditional transcriptome analysis, we mainly found that the Wnt/β-catenin, TGF-β, and FGF signaling pathways play a very important role in the process of scale development.
At present, only a few genes related to scale development have been functionally verified in zebrafish and most genes have not been functionally verified to prove their impact on scale development. Moreover, with the development of gene editing technology, it is possible to use gene knockout to assess gene function. In future work, our focus is to identify the genes related to scale development using gene knockout.

Conclusions
In this study, four different stages of zebrafish scale development were determined by alizarin red staining and the skin transcriptome was studied regarding these four stages. Using WGCNA analysis, we found 54 hub-genes highly related to scale development. In addition, combined with transcriptome analysis, we found the TGF-β, Wnt/β-Catenin and FGF signaling pathways, and the highly expressed Scpp family (Scpp7, Scpp6, Scpp5, and Scpp8), the Fgf family (Fgfr1b and Fgfr3), Tcf7, Wnt10b, Runx2b, and Il2rb genes to be involved in scale development. In conclusion, this study provides a powerful reference for further study on the molecular regulation mechanism of gene regulation, morphogenesis, and the function of related genes in fish scale formation.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/fishes7020064/s1, Figure S1: Compared with stage I, the biological process and pathways enriched in stage III and stage IV; Table S1: WGCNA analysis of 54 highly connected genes in the turquoise module; and Table S2: According to the results of WGCNA and transcriptome analysis, the FPKM values of 21 candidate genes related to zebrafish scale development were identified.