Transcriptome analysis reveals the genetic basis underlying the development of skin appendages and immunity in hedgehog (Atelerix albiventris)

The expression of hair features is an evolutionary adaptation resulting from interactions between many organisms and their environment. Elucidation of the mechanisms that underlie the expression of such traits is a topic in evolutionary biology research. Therefore, we assessed the de novo transcriptome of Atelerix albiventris at three developmental stages and compared gene expression profiles between abdomen hair and dorsal spine tissues. We identified 328,576 unigenes in our transcriptome, among which 4,435 were differentially expressed between hair- and spine-type tissues. Dorsal and abdomen skin tissues 5 days after birth were compared and the resulting DEGs were mainly enriched in keratin filament, epithelium cell differentiation, and epidermis development based on GO enrichment analysis, and tight junction, p53, and cell cycle signaling pathways based on KEGG enrichment analysis. MBP8, SFN, Wnt1 and KRT1 gene may involve in the development of hedgehog skin and its appendages. Strikingly, DEGs in hair-type tissues were also significantly enriched in immune-related terms and pathways with hair-type tissues exhibiting more upregulated immune genes than spine-type tissues. Our study provided a list of potential genes involved in skin appendage development and differentiation in A. albiventris, and the candidate genes provided valuable information for further studies of skin appendages.


Results
Illumina sequencing and de novo assembly. To identify the transcriptome and molecular mechanisms governing skin appendage development and differentiation, we analyzed temporal changes in transcript abundance of A. albiventris (Fig. 1a). RNA sequencing generated 82.2 ± 6.0 G (mean ± SD) read pairs for 16 skin tissue samples (Supplementary Table 1). After quality trimming, 96.8 ± 0.4% of reads were retained, indicating a highquality dataset (> 90% reads with ≥ Q30). The de novo assembly using Trinity revealed 328,576 unigenes, which were used for subsequent analysis. We evaluated biological reproducibility by individually comparing biological replicates through principle component analysis (PCA) (Fig. 1b). As expected, similar gene expression patterns in both hair-and spine-type samples were observed. PCA analysis also demonstrated more separation between the two skin appendage types. A neural network graph based on self-organizing feature map (SOM) analysis revealed dynamic transcriptional changes at the individual stages of appendage development in A. albiventris (Fig. 1c).
Functional annotation and classification of unigenes. Use of the NR and Swiss-Prot databases yielded reliable protein annotations for approximately 42% of unigenes. Absolute unigenes were annotated using seven major public protein databases (Table 1), of which the greatest number of matches (68.3%) was obtained using the NT database. Of 328,576 unigenes assembled, 243,444 genes (74.1%) exhibited a positive match against at least one database. Thus, our gene annotation results from A. albiventris transcriptome were considered high quality.
In our study, 70,833 unigenes were assigned to 56 sub-categories of GO terms belonging to the following three main categories: biological process (BP), cellular component (CC), and molecular function (MF). These main categories included 20, 20, and 10 sub-term categories, respectively (Fig. 2). The most enriched GO terms were related to cellular and metabolic processes (BP), cell and organelle (CC), and binding and catalytic activity (MF). Importantly, we found some GO terms (level 4) in the BP category related to skin appendage development and differentiation, including (positive/negative) regulation of cell differentiation, regulation of cell proliferation, reproductive structure development, cell development, and ovarian follicle cell development.
Analysis of differentially expressed genes in skin appendage tissues. Analysis of DEGs in skin appendage tissues finally identified 4,435 DEGs in the three developmental stages, and 816 DEGs shared by the two types of tissues (Fig. 4a-c). As shown in the Venn diagram (Fig. 4d), more DEGs were found in tissues with spines (HH1Y and HS1Y) than in tissues without spines (HH1N and HS1N).
To validate expression patterns indicated by the transcriptome data, 10 differentiation-related DEGs were selected for RT-qPCR analysis including KRT1, LEF1, RSPO2, ARRB, TGFB2, GRB2, SFN, TCF7, Table 2). Expression trends determined by RT-qPCR significantly correlated with the RNA-Seq data ( Fig. 5a-c). Expression of these 10 genes at the three developmental stages was also analyzed, revealing expression patterns similar to those determined by RNA-sequencing. Expression of LEF1, TGFB2, SFN, and WIF1 at stages I-III increased significantly. However, expression of KRT1 and RSPO2 gradually decreased (Fig. 5d). Overall, both approaches confirmed the observed DEG trend patterns, indicating the accuracy of the transcriptome data and de novo RNA-Seq data.
Differentially expressed genes related to skin appendage development. In order to further understand the spine development mechanism of the hedgehog, we conducted DEG analysis before and after spine development. In total, 28 DEGs were identified between HS1N and HS1Y, among them, 11 genes were  Table 3). Based on GO and KEGG annotation, we found one downregulated gene of unknown open reading frame (LOC103122410) related to cell differentiation and multi-cellular organismal development in the BP category, and one downregulated keratinassociated protein like gene (LOC103118355) related to keratin filament. In order to identify genes that activate the development of hedgehog abdominal hair, we analyzed DEGs before and after the occurrence of abdominal hair. A total 82 DEGs were identified between HH1N and HH1Y, of which 62 genes were downregulated, and the remaining 20 genes were upregulated (Supplementary Table 4). Notably, LOC103122410 was also present in these DEGs; in addition, keratin gene KRT2 was highly expressed in HH1Y. We speculate that these genes play an important role in the development of hedgehog abdominal hair and spines.
Based on annotation of the A. albiventirs transcripts with the GO database, the DEGs were significantly enriched in keratin-related terms. Keratin filament is an important component of skin appendages; the main DEGs involved in this process included KRT1, KRT2, MYH, DES, BMP8, SHH, and several KRTAP genes ( Fig. 6a,b). In addition, CSTA was involved in the biological process of keratinocyte differentiation. In a directed acyclic graph diagram associated with this term, cell differentiation, epithelial cell differentiation, and epidermis development were enriched, and the main DEGs involved in this process included KRT1, KRT2, MBP8, FGF, Wnt10, and MYH (Fig. 6b). The only significantly enriched KEGG pathways among all DEGs comparing HH5 and HS5 tissues were tight junction and p53 signaling pathways ( Fig. 6c; Table 2). In addition, eight pathways related to skin appendage development and differentiation were enriched, including cell cycle, MAPK, Rap1, Hippo, VEGF, TGF-beta, and PPAR signaling pathways. It is noteworthy that SFN, BMP8, Wnt3, Wnt10, MYH, and SFN were found in both annotation results, suggesting they are closely related to differentiation of keratinocytes and epidermal cells.
In our data, keratin-related genes were highly abundant and also significantly overexpressed in the dorsal spine-type tissues. In the hedgehog transcriptome, more than 3,000 transcripts were annotated to 5 keratin genes However, in our transcriptome data, KRT1 and KRT10 were not co-expressed, KRT1 was significantly upregulated in spine-type tissue, while KRT10 was down-regulated. Through alignment of KRT1 gene sequences from 20 species, we found that the hedgehog KRT1 gene sequence has a deletion of 45 bp compared with other species (Fig. 6d). The two hedgehog species were isolated into one clade in the bayesian phylogenetic tree (Fig. 6e), suggesting that KRT1 may be one of the important gene involved in the development and differentiation of hedgehog skin appendages.
immune-related genes in epidermis of hedgehog. In this study, we also found many immune-related genes comparing HH5 and HS5. After KEGG enrichment analysis, we found no significant enrichment of immune-related pathways in spine-type tissues, and fewer immune-related pathways and genes than in hair-type tissues. Hair-type tissues exhibited 19 enriched immune-related pathways, of which cytokine-cytokine receptor interaction, intestinal immune network for IgA production and cell adhesion molecules (CAMs) pathway were significantly enriched, and including 6 immune-related genes: CD27, CCR7, CCL5, LTB, PIGR and MHC2 (Fig. 8a,b). In addition, we observed similar GO enrichment analysis results. The number of terms and genes enriched in hair-type tissues was significantly higher than in spine-type tissues, of which immune response and immune system process terms were enriched considerrably, and 6 genes related to immunity: MHC2, CCL5, END3, LTB, CXCL5 and C1QA (Fig. 8c,d).

Discussion
The genetic basis of morphological variation, both within and between species, provides a major topic in evolutionary biology. 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 21 . Specifically, mammals develop a coat containing many distinct types of hair. Such diversity is associated with molecular and signaling pathways that drive formation and induction in a specific spatial and temporal manner 22 . Reciprocal interactions between epithelial and mesenchymal tissues constitute a central mechanism that determines the location, size, and shape of organs 23 .
In the current study, we aimed to explore the genetic basis for hedgehog skin appendage differentiation and development and the resulting expression of the spine trait. We identified 328,576 unigenes in our transcriptome, all of which were annotated in 7 databases (Table 1). Taken together, our de novo assemblies revealed higher quality compared with previous studies 24,25 that formed the foundation of all our subsequent analyses. According to GO and KEGG annotation analysis, a total of 70,833 and 36,238 unigenes were mapped to 56 sub-categories and 232 biological pathways, respectively (Figs. 2, 3). We identified some of the key pathways involved in skin development; these annotations may provide a valuable resource for further understanding the specific functions and pathways in A. albiventris.
Newborn hedgehogs begin to develop hair/spines approximately 2 h after birth. We found 6 shared genes (APOE, COX2, COX3, FCN, RP-L11e, SH3GL) through analysis of DEGs before and after hair/spine development, indicating that these genes are essential regulatory genes in the development of skin appendages, whether hair or spine. Further, we compared the gene expression profiles of spine-and hair-type tissues to systematically assess key regulatory genes for spine and hair differentiation in A. albiventris. SFN, upregulated in spine-type tissues, is a regulator of mitotic translation that interacts with a variety of translation and initiation factors 26 , is enriched in the p53 signaling pathway, and has a role in keratinocyte differentiation and skin barrier establishment in the BP category (GO). SFN plays an important role in maintaining hair follicle development, especially affecting the formation of hair shaft structure 27 . Bu et al. 28 determined that SFN gene and protein were significantly highly expressed in the thicker, longer, and harder skin of wool, indicating that SFN was involved in the regulation of wool character development. In addition, we identified FGF as the upregulated DEG enriched in the MAPK and Rap1 signaling pathways. These pathways often act together by forming signaling loops during organogenesis 29 and induce the most fundamental biological processes, such as the formation of periodic patterns. Hebert et al.
(1994) demonstrated that FGF plays an important role in regulation of the hair cycle growth, and functions as an inhibitor of hair elongation by promoting progression from anagen stage, the growth phase of the hair follicle, to Scientific RepoRtS | (2020) 10:13920 | https://doi.org/10.1038/s41598-020-70844-y www.nature.com/scientificreports/ catagen stage, the apoptosis-induced regression phase 16,30 . Hence, we believe that SFN and FGF may be important regulators that affect spine development and differentiation in A. albiventris. At present, keratin and keratin associated proteins are the main structural proteins in model organisms 31,32 , which determine the basic properties of hair, such as fineness, length, hardness and luster, etc 33 , and in the process of hair formation, a series of keratin genes also determine the process of hair follicle cell differentiation 34 . In addition, many skin and hair-related diseases are also associated with mutations and expression of keratin family genes [35][36][37] . In this transcriptome study, it was found that there were significant differences in keratin gene expression in different skin regions at different development stages, such as KRT1, KRT2, KRT5, KRT10, KRT14, KRT6A, KRT75, KRATP9-2, KRTAP13-1, KRTAP19-2, etc. Parson et al. 38 and Rogers et al. 39 found that the differential expression of keratin and keratin associated proteins was related to the fineness and density of hair at different development stages. Khan et al. 40 found that the unique hair-related phenotypes, such as scales (armadillo) and spines (hedgehog), were correlated changes in the expression variation probably also influences hair diversification patterns. And the number and size of differences in gene family replication, loss and pseudogenization are very important for the evolutionary history of mammals 41,42 .
We also compared KRT1 gene sequences of 20 mammals with that of the hedgehog, in which we identified a 45-bp deletion in the latter. Keratins, the major structural proteins of epithelia, are a diverse group of cytoskeletal scaffolding proteins that form intermediate filament networks and provide structural support to keratinocytes that maintain skin integrity 43 . In general, KRT1 is highly conserved in mammals, however, molecular defects in keratin intermediate filament-related genes can cause keratinocyte and tissue-specific fragility, accounting for a large number of genetic disorders in skin and its appendages 44,45 . Therefore, whether the deletion phenomenon in hedgehog KRT1 has special significance for spine development and differentiation in the hedgehog warrants further research and verification.
Lastly, we analyzed immune-related genes in the skin transcriptome of A. albiventris. We found that hair-type skin had more immune-related genes than spine-type skin. Choo et al. found that interferon epsilon (IFNE), which was exclusively expressed in epithelial cells and was important to mucosal immunity, was pseudogenized in pangolins. They proposed that scale development provided protection against injuries or stress and reduced pangolin vulnerability to infection, thus protection by scales on the pangolin body compensated for the low immunity of this species to a certain extent 46 . From our current data, we speculate that there may be significant differences in the immune function of skin with different appendage types, and that the occurrence of spines may be an innovative physical protection for dorsal skin with relatively low immunity.  In the current study, we conducted a comprehensive transcriptome analysis of A. albiventris to explore the genetic basis for the growth of skin appendages. Transcriptome analysis provided a rich list of unigenes expressed in hair-and spine-type tissues at three different developmental stages, of which 328,576 unigenes and 3,598 DEGs were identified. Candidate genes were identified that are likely involved in the regulation of hair and spine growth and differentiation. The knowledge acquired in this study of the molecular and signaling pathways related to hair and spine expression greatly contributes to the current genetic resources for the hedgehog and mammalian species that harbor shaggy appendages, as well as traits that could be potentially exploited for curing skin diseases of other animals, even humans.

Methods
Biological samples. According to observations of the growth characteristics of hedgehogs, hairs and spines on the back and abdomen begin to appear approximately 2 h after birth; therefore, we obtained 8 hedgehogs at 3 different stages of appendage development from a commercial animal farm (Dongguan City, China). Hedgehog hair (HH) and spine (HS) tissues were collected from 2 specimens within 2 h of birth (Stage I, HH1N/HS1N), 3 specimens after 2 h but within the first birth day (Stage II, HH1Y/HS1Y), and 3 specimens 5 days after birth (Stage III, HH5/HS5). Prior to skin sampling, follow the standard operating procedures for euthanasia, conduct euthanasia for experimental animals, minimize the pain of small animals, do not affect the results of animal experiments, and shorten the time of death as far as possible. Sixteen skin tissue samples representing the two types of appendages (abdomen hair-type and dorsal spine-type) were rapidly excised, immediately snap-frozen on dry ice, and stored at − 80 °C until RNA extraction.
RNA extraction and RNA-seq. Total RNA from each tissue sample was extracted using the RNeasy Kit (Qiagen, Hilden, Germany). RNA purity was determined using a NanoPhotometer spectrophotometer (Implen, Inc., Westlake Village, CA, USA). All unigenes were annotated using 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 47,48 ; Gene Ontology (GO) database. Differentially expressed genes (DEGs) were identified by comparing gene expression levels between samples (or sample groups).  ethics approval and consent to participate. All animal procedures in the study were approved by the ethics committee for animal experiments at the Guangdong Institute of Applied Biological Resources (reference number G2ABR20170523) and followed basic principles. We confirm that all methods were performed in accordance with relevant guidelines and regulations.

Data availability
Data analyzed in the current study are included within the article and its supplementary material. All unigene sequences from A. albiventris have been deposited in the GenBank Sequence Read Archive (SRA) under accession number PRJNA561241 for SUB6195278. We have uploaded supplemental material to figshare via the GSA Portal. Supplementary Table 1 contains the summary of sequencing statistics for the transcriptomes; Supplementary Table 2 contains the summary of primer information used in real-time PCR analysis; Supplementary Table 3 contains the DEGs between HS1Y and HS1N; Supplementary Table 4 contains the DEGs between HH1Y and HH1N; Supplementary Table 5 contains the DEGs between HH5 and HS5; Supplementary Table 6 contains the species name and NCBI serial numbers of 19 additional species for phylogenetic analysis.