The procedures for animal housing, care, and collection of skin samples were approved by the Animal Ethics Committee (2017[050]) at Shanxi Agriculture University (Taigu, Shanxi, China; Approval number: 2023008) and performed according to the Committee guidelines.

The CDK5-knockdown mouse models were prepared according to previously described instructions [12]. Three healthy wild-type and three CDK5-knockdown mice were randomly selected, and three pieces of skin (2×3 cm) from their backs were collected under local anesthesia and immediately stored in liquid nitrogen. Total RNA was extracted from the samples using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA integrity was evaluated by gel electrophoresis, and the concentration was measured by absorbance at optical density (OD_260/280) using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

After total RNA extraction and DNase I treatment, magnetic beads bound to oligo-dT were used to isolate the mRNAs. These were mixed with a fragmentation buffer to break the mRNAs into short fragments. cDNA was synthesized from the mRNA fragment templates. Short fragments were purified and resolved using the Qiagen EB buffer (Qiagen, Hilden, Germany) for end preparation and single-nucleotide adenine addition. These fragments were then ligated to adapters and suitably sized fragments were selected for polymerase chain reaction (PCR) amplification using agarose gel electrophoresis. During these quality control steps, an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and an ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) were used for mRNA quantification and qualification. Finally, the library was sequenced using an Illumina HiSeq 2000 (Illumina, San Diego, CA, USA) to generate raw reads, which were filtered into clean reads by removing adaptors and low-quality reads, and then aligned to the reference sequences. These alignment data were used to calculate the distribution of reads on the reference genes and mapping ratios for further analyses, including gene ontology (GO) enrichment analysis, Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis, and transcription factors identification.

UniGene assembly was performed using the Trinity software (http://www.genomics.cn). BLASTX alignment (e-values <0.00001) between the unigenes and several protein databases (nr, Swiss-Prot, KEGG, and clusters of orthologous genes) was performed, with the best alignment results used to determine the sequence direction of the unigenes. GO functional annotation was based on the nr annotation, and Blast2GO (http://www.blast2go.com) was used to assign GO annotations. WEGO (http://wego.genomics.org.cn/cgibin/wego/index.pl) was used to perform the GO functional classifications of all unigenes.

A rigorous algorithm was used to identify genes in the skin that are differentially expressed between CDK5-knockdown and wild-type mice. The algorithm consisted of a threshold false discovery rate value of ≤0.001 and a reads per kilobase of transcript per million mapped reads ratio of ≥2. Differentially expressed genes were mapped to each term in the GO database (http://www.geneontology.org/) and the number of genes assigned to each GO term was calculated. The calculated p-values were Bonferroni corrected and corrected p-values ≤0.05 were determined to be significant. GO terms that fulfilled this condition were defined as significantly enriched. The differentially expressed genes was also mapped to terms on the KEGG pathway database (http://www.genome.jp/kegg/pathway.html) to reveal the biological pathways that differed between CDK5-knockdown and wild-type mice.

To validate the sequencing data, 10 genes were randomly selected from a list of differentially expressed genes using the quantitative real-time PCR (qRT-PCR). Total RNA from six samples used for RNA sequencing was used for this analysis. One microgram of DNase-treated RNA was converted to cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China). The cDNA was used for qRT-PCR quantification using mRNA-specific primers (Table 1). β-Actin was used as an endogenous control. The qRT-PCR was performed in triplicate using a Stratagene Mx3005P system (Stratagene, La Jolla, CA, USA). Each 10 μL PCR reaction volume included 5 μL of SYBR Premix Ex TaqTM II (TaKaRa, China), 0.2 μL of specific forward primer, 0.2 μL of reverse primer, 0.2 μL of ROX reference dye, 1 μL of 10-fold diluted cDNA, and 3.4 μL of water. The cycling parameters were 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, 56°C or 58°C for 30 s, and 72°C for 15 s. Melting curve analysis was performed at each amplification step. The mRNA abundance of each target gene was quantified using the comparative threshold cycle (CT) method and normalized to β-actin [13]. Differences in the mRNA abundance of genes were determined using analysis of variance.

To identify the RNAs expressed in mouse skin, wild-type and CDK5-knockdown skin mRNA libraries were generated using high-throughput sequencing. A total of 98,801,386 and 96,503,448 raw reads were obtained from the wild-type and CDK5-knockdown mice, respectively. After filtering low-quality reads, 86,182,330 (87.23%) and 86,398,588 (89.53%) clean reads remained in the wild-type and CDK5-knockdown datasets, respectively. Using an algorithm based on a previously described method [14], 8,002 known genes were identified as differentially expressed between the skin of CDK5-knockdown and wild-type mice. Of these, 4,344 were downregulated (≤2-fold) and 3,658 were upregulated (≥2-fold) in CDK5-knockdown mice compared to wild-type (Figure 1; Supplementary Table S1).

Genes differentially expressed between the skin of CDK5-knockdown and wild-type mice were grouped into 49 classes based on their putative functions. For subsequent GO analysis, the differentially expressed genes were characterized at the biological process (48.98%), cellular component (32.65%), and molecular function (18.37%) levels (Figure 2). A total of 318 unknown genes were also identified as differentially expressed, of which 171 were downregulated (≤2-fold) and 147 were upregulated (≥2-fold) in the skin of CDK5-knockdown mice compared to the skin of wild-type mice (Supplementary Table S2).

More than 150 mutations affecting hair and skin pigmentation have been identified in mice. They are well dispersed throughout the mouse genome and are found at more than 50 distinct genetic loci [15]. Known hair color genes are typically classified into six functional groups: melanocyte development, melanosome components and precursors, melanosome construction and protein routing, melanosome transport, eumelanin and pheomelanin production, and other wider systemic effects [14]. In this study we found that 80 of these hair color genes are expressed in mouse skin. Our analysis revealed that 61 of these genes were expressed at lower levels in the skin of CDK5-knockdown mice than in wild-type mice, and 19 genes were expressed at higher levels (Figure 3; Supplementary Table S3). We also found that all the genes involved in encoding the components of melanosomes and their precursors were expressed at lower levels in the skin of CDK5-knockdown mice than in the skin of wild-type mice (Supplementary Table S4). Most genes involved in the production of eumelanin and pheomelanin also showed lower expression levels in the CDK5-knockdown mice. Among the hair color genes with lower expression in the skin of CDK5-knockdown mice, tyrosinase-related protein 1 (TYRP1) showed the largest decrease in expression between CDK5-knockdown mice and wild-type mice, followed by paired box3 (PAX3), dopachrome tautomerase (DCT), adenomatous polyposis coli (Apc), and lysosomal trafficking regulator (Lyst) genes (Table 2).

Of the 8,004 genes that were differentially expressed between the skin of CDK5-knockdown mice and wild-type mice, 234 had specific KEGG annotations. The top 20 pathways that were most highly expressed in the skin of mice included the cGMP-PKG signaling pathway, calcium signaling pathway, and MAPK signaling pathway (Figure 4). The genes identified in the mouse skin transcriptome related to pigmentation and melanogenesis and their relative level of differential expression between the CDK5-knockdown and wild-type mice are shown in Table 3 and Figure 5.

To validate the transcriptome sequencing results, we randomly selected 10 genes from a population of genes known to be involved in hair color for validation using qRT-PCR. Transcriptome sequencing analysis revealed that these genes were also differentially expressed between the CDK5-knockdown and wild-type mice. The results of qRT-PCR analysis showed that the mRNA expression of the selected genes was consistent with the transcriptome sequencing data (Figure 6). Among the differentially expressed genes, TYRP1 showed the greatest difference in expression between CDK5-knockdown and wild-type mice, again supporting the sequencing data.