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). Overall, 85.6% of unigenes were successfully annotated with significant hits (e-value < 1E–10). Absolute unigenes were annotated using seven major public protein databases (Table 1).
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.
In addition, we annotated all unigenes using KEGG pathway analysis to understand their top-level function and role in the hedgehog biological system. Overall, 36,283 unigenes were mapped to 237 biological pathways clustered in five KO hierarchies in the functional ortholog system. Among these KO hierarchies, 15 unique KEGG pathways represented cellular processes and 30, 22, 98, and 72 pathways represented environmental information processing, genetic information process, metabolism, and organismal systems, respectively. The largest KO group was translation (5743 unigenes, 13%), followed by signal transduction (5551 unigenes, 12.6%), and endocrine system (3223 unigenes, 7.3%). Further, a small proportion of unigenes were annotated into categories related to skin appendage development, such as SHH (95 unigenes), p53 (161 unigenes), apoptosis (189 unigenes), and TGF-b (206 unigenes) signaling pathway (Fig. 2b).
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 combinations: EMH1 (Emb vs. HH1N), EMS1 (Emb vs. HS1N), 1N (HH1N vs. HS1N), 1Y (HH1Y vs. HS1Y), 5D (HH5 vs. HS5), 14D (HH14 vs. HS14), and 43D (HH43 vs. HS43).
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 were always identified: BP: GO:0042254 Ribosome biogenesis, GO:0006412 Translation; MF: GO:0003735 Structural constituent of ribosome; CC: GO:0005840 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.
Establishing gene modules underling the key developmental DEGs by 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 spine-differentiation 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). Analysis of the GO terms in each module revealed signaling activation and transmission-related genes (magenta module, C1), including genes for the regulation of cellular catabolism (ZNF365, TPM2, TPM1, and KRT2), cell activation (MUSTN1, TGFB3, GNDF, and PDGFC), and cell proliferation (MYH and DES), highly expressed at Stage I and actively preparing the follicular cells for hair growth. Genes enriched in the metabolic and biosynthetic processes (tan module, C2) exhibited increased expression at Stage II compared with that at C1 and included genes related to the primary metabolic process, cellular biosynthetic and organic substance process, and cell. During Stages III–IV of development, the expressed genes were associated with further development (hardening, thickening, etc.) of the spine, including organic substance transport, macromolecule/cellular/protein localization, and cell cycle. During that time, most of the genes in the green module (C3) were related to keratin (KRT1, KRT2, and keratin-associated proteins). GO terms in the brown module (C4) were involved in the primary metabolic and biosynthetic process, including cellular and organic substance metabolism (USP, TOP2, TGFB2, SSH, and TAF), location, gene expression, and transport (CEBPD, MYH, TCN2, and THBS2S), and cellular biosynthesis (MT1, FCN, ISYNA1, KBTBD5, and GALNT).
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 anti-apoptosis 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). Some key immune molecules were identified form transcripts of A. albiventris, such as ACTB_G1, AKT, COL1AS, GNAS, MHC2, P38, MYLPF, MAPK1_3 etc. The details of immune-related unigenes were listed in Supplementary Table 6. Among these annotated immune genes, 16 genes were also found in Manis javanica, such as ARPC2, SHC1, GNB3, CD5, MHC2, PLG etc.