Developmental transcriptome and genes related to skin appendage differentiation in hedgehog (Atelerix albiventris)

Hedgehog spines are skin appendages that have evolved as a result of the interaction of their skin with the environment. However, such a differentiation mechanism during skin development leads to a high skin appendage diversity, the origins of which are still not fully understood. Spine-skin and hair-skin offer a natural model for the analysis of the genomic basis for the evolution of epidermal appendage formation. We performed a comparative transcriptomic analysis of hedgehog ( Atelerix albiventris ) at multiple developmental stages, and tried to explore the genetic basis for differentiation and development and the resulting expression of the spine-trait. A total 15,158 differentially expressed genes (DEGs) were identified. We determined the gene modules and programs corresponding to the various phenotypic traits at different developmental stages by WGCNA analysis. Objective analysis of gene module expression revealed that HIPPO, TGFB, MAPK and Wnt signaling pathways regulate the activation and cell proliferation and differentiation at the skin-appendage development stage. Further, the key genes encoding keratin, FGF, TEAD, and other proteins regulate molecular localization and the cell cycle for hair development and differentiation. Finally, we found a number of highly expressed immune genes in the skin, suggesting that hedgehog spines, unlike pangolin scales, have evolved independently to protect against predators rather than compensate for low autoimmune immunity. and were from seven comparison combinations from the transcriptome datasets. Candidate genes related to five key signaling pathways were identified, which are likely involved in the regulation of the growth and differentiation of the hair and spine. The knowledge of molecular and signaling pathways greatly contributes to the current genetic resources for hedgehog and certain mammalian species that harbor shaggy appendages as well as traits that could potentially be exploited commercially for curing skin diseases of other animals, even human.

signal [7][8][9]. Another secreted protein in the follicular placode that plays a major part in epithelial-mesenchymal signaling is Sonic hedgehog (SHH) [10,11]. SHH depends on WNT signaling, and this protein is required for the proliferation of follicular epithelium and the development of the dermal condensate into the dermal papilla [4]. In addition, many key genes related to skin appendage development have been identified, such as TGF-α receptor (epidermal growth factor receptor) and the transcription factor ETS2, which regulate the shape of hair follicles and are responsible for hair follicle architecture and wavy hair [12,13], and HOXC13 and KRT75, which control the differentiation of the hair shaft [14][15][16][17]. Further, EDA and EDAR interact with members of the bone morphogenetic protein (BMP) family, some of which inhibit follicle development to establish follicle patterning [2,[18][19][20].
The abovementioned studies primarily focused on mammalian hairs and scales and little is known about changes in gene expression during the development of spine appendages in mammals. Therefore, here, we aimed to perform a detailed analysis of the transcriptome of the skin of abdomen hair and dorsum spines of hedgehogs (Atelerix albiventris) by comparing groups of transcripts differentially expressed at different hedgehog developmental stages. We identified key candidate genes related to the development and differentiation of skin appendages. Their importance should be verified in future studies. Nonetheless, the presented novel information will be widely applicable in many fields, e.g., skin disease, skin immune genes, also providing insight into the molecular mechanisms of skin appendage development in general.

Biological samples
Ten hedgehogs at six developmental stages (Stage 0: Emb, the embryo; Stage I:  Table 1). For the skin appendage collection, the animals were first anaesthetized by using diethyl ether and then killed via cervical dislocation. Two embryos and 20 skin tissues representing the two types of appendages (abdomen spine-type and dorsal hair-type) were rapidly excised, immediately snap-frozen on dry ice, and stored at -80 °C until RNA extraction.

RNA extraction, library preparation, and Illumina sequencing
Total RNA from each tissue sample was extracted using the RNeasy Kit (Qiagen, Germany). RNA purity was determined using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNANano 6000 assay kit and the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). RNA concentration (ng/ml) was determined using the Qubit® RNA assay kit and Qubi® 2.0 Flurometer (Life Technologies, CA, USA).
Qualified RNAs were used for cDNA library construction and sequencing. All the sequencing stages were conducted at Novogene sequencing company (Beijing, China).

De novo transcriptome assembly and functional annotation
Sequencing reads from 22 samples were first evaluated using FastQC to assess the initial read quality. Clean reads were obtained by removing reads containing adapters or poly-N, and low-quality reads from raw data. Simultaneously, Q20, Q30, GC content, and sequence duplication level of the clean data were calculated. All downstream analyses employed clean data of high quality. The trimmed paired-end reads were resynchronized (https://github.com/enormandeau/Scripts) and assembled de novo using the Trinity assembler [21].
All unigenes were annotated in the course of analysis with the basic local alignment search tool (BLASTX), considering hits with e-values of 1E−5 against seven databases: NCBI non-redundant protein sequences (Nr); NCBI non-redundant nucleotide sequences (Nt); Protein family (Pfam) database; Clusters of Orthologous Groups of Proteins (KOG/COG) database; Swiss-Prot (a manually annotated and reviewed protein sequence database); KEGG Ortholog (KO) database; Gene Ontology (GO) database. The best-aligned results were used to determine the sequence direction of the transcripts. In case of discrepancies of analyses using different databases, the sequence direction of the transcript was determined according to the prioritization order of NR > Swiss-Prot > KEGG > COG.

Differential expression analysis
Differentially expressed genes (DEGs) were identified by comparing gene expression levels between samples (or sample groups). Bowtie2 [22] was used to align the clean reads with all unigenes, and RSEM [23] was used to calculated the gene expression levels for each sample. Differential gene expression was analyzed using For samples with biological replicates, differential expression analysis of two conditions or groups was performed using the DESeq R package (1.10.1). DESeq [24] provides statistical approaches for determining differential expression in digital gene expression datasets using a model based on negative binomial distribution.
The resulting P-values were adjusted using the Benjamini and Hochberg's approach to control the false discovery rate. Genes with an adjusted P-value < 0.05 identified by DESeq [25] were deemed to be differentially expressed. For samples without biological replicates, prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by using edgeR program package [26] and a single scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq R package. The P-values were adjusted using the q-value [27]; q-value < 0.005 and |log2(foldchange)| > 1 were set as the threshold for significantly differential expression.
Further, KOBAS software [28] was used to test the statistical enrichment of DEGs in KEGG pathways. GO enrichment analysis of DEGs was implemented by the GO seq R package, based on Wallenius non-central hyper-geometric distribution [29], adjusting for gene length bias.

Quantitative real-time reverse-transcription PCR (RT-qPCR)
RNA was extracted from 22 skin appendage tissues using TRIzol reagent. Then, cDNA was synthesized by using the Toyobo reverse-transcription kit. Eleven genes related to hair and/or spine development were selected for analysis based on the functional annotation data, and fluorescent quantitative PCR primers were designed accordingly (Supplementary Table 2). RT-qPCR was performed in a 20-ml reaction volume, with four technical replicates for each sample, using the TransStart® Top

Transcriptomic atlas for developing hedgehog skin appendages
To identify the transcriptome and molecular mechanisms governing the development and differentiation of skin appendages, we analyzed the temporal changes in the transcript abundance of A. albiventris (Fig. 1a). RNA sequencing generated 81.4 ± 6.6 G (mean ± SD) read pairs for 22 skin tissues (Supplementary Table 3). After quality trimming, 96.8% ± 0.6% of reads were retained, indicating a high-quality dataset (>90% reads with ≥Q30). The de novo assembly using Trinity revealed 32,8576 unigenes, which were used for subsequent analysis. The main length distribution of 94.15% of unigenes was in the range of 0.2-2 kb (Fig. 1b).
A neural network graph based on SOM analysis revealed dynamic transcriptional changes occurring at the individual developmental stages of A. albiventris (Fig. 1c).
Among the 56 annotated GO terms, cellular process (38,648 unigenes) was the most common annotation in the biological process (BP) category and metabolic process (31,933 unigenes) and single-organism process (28,291 unigenes) were the next most abundant GO terms. In addition, a low proportion of unigenes were annotated as terms related to skin appendage development, such as cellular component organization (9844 unigenes), developmental process (2127 unigenes), cell killing (79 unigenes), and cell aggregation (1 unigene) (Fig. 2a). In the molecular function (MF) and cellular component (CC) categories, binding (31,320 unigenes) and cell (25,455 unigenes) were the most common annotation terms. Details of GO classification of the A. albiventris unigenes are provided in Supplementary Table 4.

Characterization of the developmental cycle of skin appendages in hedgehog
To systematically assess the key regulatory genes related to the development and differentiation of spine and hair in A. albiventris, we compared the gene expression profiles of spine-type and hair-type tissues in A. albiventris at six developmental stages (Stages 0-V). Highly expressed DEGs and those with type-specific expression in the skin appendage were screened using the following seven comparison Overall, to identify the core transcriptome, we screwed 22,587 shared unigenes expressed at each developmental stage. GO enrichment analysis of these shared unigenes revealed no molecular function change in the most enriched terms during the six sequential developmental stages of hair-type tissue (Fig. 3a). The following Ribosome. However, in the spine-type tissue, except for BP, the most enriched terms in the MF and CC categories changed during stages III-IV (Fig. 3b). We also screwed and annotated the top 100 genes based on the FPKM values in tissues at each stage. In general, protein-coding genes related to skin appendage development (e.g., RPL13A and RPLP1) and keratin (e.g. KRT1 and KRT2) were highly dynamic and their expression was time-dependent throughout a developmental stage. This suggested their widespread roles in organ development. Considering the transformation of biological functions during stages III-IV, the corresponding regulatory genes with highest frequency changed from ribosomal protein-coding genes (L13a and lateral stalk subunit P1) to keratin-encoding genes (types I and II) ( Fig. 3c).
We also identified a set of 15,158 DEGs based on comparison combinations of 1N, 1Y, 5D, 14D, and 43D, and one intersection with 163 DEGs (Fig. 4a-f). According to the KEGG enrichment analysis, we screwed 27 KEGG pathways related to skin appendage development and then analyzed DEG clustering patterns in these 27 KEGG pathways for different comparison combinations (Fig. 5). As shown in the heatmap, three clusters of gene expression patterns were apparent: Cluster 1, including the groups EMH1 and EMS1; Cluster 2, including the groups 1N, 1Y, and 5D; and Cluster 3, including the groups 14D and 43D.

WGCNA analysis
We first imported hair-type tissue, spine-type tissue, and embryo sample data to R and executed WGCNA analysis with an emphasis on genes that underwent spinedifferentiation changes. WGCNA identified 23 modules of co-expressed genes (Fig.   6a). The brown, tan, magenta, green, and purple modules were the top five modules related to the spine-type trait, and the red, light cyan, cyan, light green, and pink modules were the top five modules related to the hair-type trait (Fig. 6b).
GO analysis of the top eight modules revealed several key biological processes that took place during spine and hair development (Fig. 6c) We also performed KEGG analysis of modules related to the two skin appendage traits. We screwed five key signaling pathways involved in spine development and differentiation, namely, the HIPPO, TGFB, MAPK, cell cycle, and Wnt signaling pathways. Some key genes that promote cell reproduction, differentiation, and antiapoptosis were significantly enriched in these pathways, e.g., FGF, EGF , PAK1 , TEAD, CDK6 , and TCF (Fig. 7a, Table 2). Further, five key signaling pathways were identified as related to hair development, i.e., TGFB, p53, apoptosis, cell cycle, and tight junction signaling pathways (Fig. 7b). Genes that regulate apoptosis, arrest cell development, and reduce cell numbers were enriched in these pathways, e.g., CCND3, CCNB , CD1 , ACTB_G1 , and TUBA.

Verification of important DEGs related to hair and/or spine development and differentiation-related molecules using RT-qPCR
To validate the expression patterns indicated by the transcriptome data, 11 differentiation-related DEGs were selected for RT-qPCR analysis. These were KRT2, LEF1, RSPO2 , ARRB , TGFB2 , GRB2 , SFN , TCF7 , CTNNB1 , HOXC13 , and WIF1 . Indeed, the expression trends determined by RT-qPCR significantly correlated with the RNA-Seq data (Fig. 8a). The expression of these 11 genes at five developmental stages was also analyzed, revealing expression patterns similar to those determined by RNA-Seq experiments (Fig. 8b). The expression of LEF1, TGFB2, SFN, and WIF1 at stages I-III significantly increased. However, the expression of KRT1 and RSPO2 gradually decreased. Overall, both approaches confirmed the observed DEG trend patterns, indicating the accuracy of the transcriptome data and de novo RNA-Seq data.

Immune-related genes in skin epidermis of hedgehog and pangolin
In the KEGG classification, we also identified 4,687 unigenes were involved in immunity related pathways. The 16 immune system related pathways were extracted, Platelet activation, Leukocyte transendothelial migration and Chemokine signaling pathway were the top three categories with higher proportions (Fig. 9).

Discussion
In mammals, the development of skin appendages, such as hair, tooth, and scale, involves complex interactions between the epidermis and the underlying mesenchyme as part of an established hierarchical morphogenetic process [30].
Specifically, mammals develop a coat containing many distinct types of hair. Such diversity is associated with molecular and signaling pathways that drive the formation and induction in a specific spatial and temporal manner [31]. Reciprocal interactions between the epithelial and mesenchymal tissues constitute a central mechanism that determines the location, size, and shapes of organs [32].
In the current study, we aimed to explore the genetic basis for hedgehog differentiation and development and the resulting expression of the spine-trait.
First, according to the comparison of the number of up-and down-regulated genes in 5 developmental stages, we found that spines development requires more genes to participate in expression and regulation. Especially from stage I to stage II, the number of DEGs increased by 30%.
Second, according to the core gene analysis, the six developmental stage in the current study indeed covered the main skin development process throughout the hedgehog life cycle. We found some key nodes in the development of hedgehog skin appendages. (1) Spines developed earlier than hairs. At stage II, the hedgehog began to grow spines on its back, but no hair on its abdomen, indicating that the 2- Third, WGCNA analysis identified several key signaling pathways related to the development of the appendage trait (Fig. 7). The first was Wnt signaling pathway which was involved in the establishment of the first dermal signal [33]. Our DEGs associated with two phenotypes both were annotated in the Wnt pathway, however, genes associated with the spine-type tissue were significantly enriched in this pathway, while hair-type related genes were not. This indicated that the activation of Wnt signaling pathway was required at the early stage of hair and/or spine development in hedgehog, similar to other mammals. As well as three b-catenin genes (CTNNB1, ICAT, and JUP), genes for FZD, LEF/TCF, and Wnt family members were upregulated in the spine-type tissue (Fig. 8). Sustained b-catenin activity leads to increased dermal thickness and proliferation of fibroblasts, that could lead to increased hair follicle substrate and dermal condensate area during development, accelerating hair follicle differentiation [4,5]. Otherwise, a null mutation in a gene encoding the Wnt effector protein LEF1 in mouse leads to a failure of formation of the vibrissae and two-thirds of the body hair follicles [34,35]. This is consistent with the findings for A. albiventris. In addition, we identified is associated with a wide range of pathologies in human. It has been reported that some diseases related to human skin appendages are regulated by the Wnt signaling pathway, e.g., tooth agenesis [36], odontoonychodermal dysplasia [37], palmoplantar hyperkeratosis with squamous cell carcinoma of skin and sex reversal [38], acne inversa [39], focal dermal hypoplasia [40], and nonsyndromic congenital nail disorder [41].
Further, DEGs related to the two phenotypes were also significantly enriched in the TGF-b signaling pathway, affecting follicle fate direction and distribution pattern.
TGFB, E2F4 , DCN , PAK1 , and ID4 were here identified as DEGs upregulated in the spine-type tissue. We hypothesize that these genes are involved in skin appendage development and follicle size. Jackowska et al. (2012) observed increased TGFB1, TGFB2, and TGFB3 protein levels in oocytes in large follicles compared with those in small follicles [42]. This is consistent with the hedgehog data from the current study. By contrast, members of the BMP family of secreted signaling molecules appear to act as inhibitors of follicle formation. In the current study, several BMP and Noggin genes were identified as DEGs, with opposite trends of expression in the hair-type and spine-type tissues. Noggin is a secretory molecule that may inhibit BMP activity and is involved in the development of hair follicle. In support of this hypothesis, mouse lacking Noggin has fewer hair follicles than wild type and retarded follicular development [43]. Further, some diseases related to human skin appendages are regulated by the TGF-b signaling pathway, e.g., the Loeys-Dietz syndrome [44] and Axenfeld-Rieger syndrome [45].
In addition, we identified some DEGs enriched in the MAPK and cell cycle signaling pathways. These pathways often act together by forming signaling loops during organogenesis [46]. They induce the most fundamental biological processes, such as the formation of periodic patterns. In the current study, FGF, EGF , CDK6 , and SFN were upregulated in the spine-type tissue. Hebert et al. (1994) demonstrated that the FGF gene plays an important role in the regulation of hair cycle growth, and functions as an inhibitor of hair elongation by promoting progression from the anagen stage, the growth phase of the hair follicle, to the catagen stage, the apoptosis-induced regression phase [16]. Hence, we believe that these genes may be important regulators that affect the development and differentiation of spine.
We identified UCP1 as a DEG in the Wnt signaling pathway. It encodes a mitochondrial protein responsible for thermogenic respiration that participates in non-shivering adaptive thermogenesis in response to temperature variation and, more generally, in the regulation of energy balance. UCP1 was significantly upregulated in the hair-type tissue, suggesting that the spiny back of hedgehog is more sensitive to cold than abdomen with hair. We also compared the KRT1 gene sequences of some mammals with those of the European hedgehog, and identified an approximately 40-bp deletion in the KRT1 gene in the latter. Keratin gene is highly conserved in mammals. Whether this deletion has a special significance for the development of hedgehog skin appendages requires further study and analysis.
At last, we found some immune-related genes in the skin epidermis transcriptome of hedgehog. By observing hedgehog captive population, we found that there was no obvious immune deficiency in this species, which was consistent with the normal expression of the detected immune genes. Choo et al. (2016) found that interferon epsilon (IFNE) which expressed in epithelial cell and mucosal immunity on pangolins skin has been pseudogeneized, so they propose that scale development was an innovation that provided protection against injuries or stress and reduced pangolin vulnerability to infection, the protection of pangolin scales on the body can compensate for the low immunity of this species to a certain extent. Interestingly, 16 of the immune genes we found on hedeghog were also annotated on pangolin.
However, neither IFNE gene nor the 16 genes we detected showed pseudogenetization in hedgehogs, which indicated that the compensatory effect of skin appendages on immunity did not occur in hedgehogs. So, hedgehog's ratchet may have evolved independently of the immune system, and its function is more physical defense.

Conclusions
In the current study, we conducted a comprehensive transcriptome analysis of A.

Data availability
The data analyzed in the current study are included within the article and its additional files. All unigene sequences from A. albiventris have been deposited in the GenBank Sequence Read Archive (SRA) under the accession number PRJNA561241 for SUB6195278.