RNA-seq transcriptomic and label-free quantitative proteomic analysis of M. javanica skin tissues
To identify the molecular mechanisms governing the development and differentiation of skin appendages, we harvested skin samples from M. javanica in biological replicates with 7 specimens including 2 embryos and 5 adult specimens at two developmental stages for two skin appendage types: embryo/adult-hair and embryo/adult-scale group (Fig. 1a). The 14 samples were divided into identical pools and subjected to label-free MS-based proteomics and RNA-seq-based transcriptomics. The adult-hair3 sample in the sequencing results was of poor data quality; thus, this sample was deleted from our transcriptome data.
The Illumina sequencing reads were deposited in the Short Read Archive as accession number PRJNA610466 for M. javanica. For the two types of pangolin skin, a total of 84,182,093 clean reads were obtained, corresponding to 90.54% of all high-quality reference genes of the M. javanica genome (NCBI BioProjects: PRJNA335369; Table 1). A total of 36,043 annotated RNAs were obtained, of which 21,302 were protein-encoding. A neural network graph based on self-organizing map (SOM) analysis20 revealed dynamic changes of gene expression occurring at the individual developmental stages (Fig. 1b). We evaluated biological reproducibility by comparing the biological replicates individually by principle component analysis (PCA) with all skin appendages (Fig. 1c). As expected, a similar expression pattern was observed for some of the genes in both hair- and scale-type samples.
Based on these results, a total of 8,854 new genes were discovered in the transcriptomic datasets (Supplementary Table 1), of which 7,800 were annotated by GO and KEGG database, 38 new genes were identified as being involved in KEGG pathways related to skin appendage development and differentiation, such as apoptosis, cell cycle, Hedgehog signalling pathway, etc. The new genes provide new information for studying the differentiation mechanisms of pangolins.
To complement the transcriptomic analyses, we carried out a systematic label-free quantitative proteomic analysis. Following MaxQuant analysis, peptide sequences were searched against the transcriptome reference dataset. In total, 2,904 proteins containing 26,472 peptides were successfully identified from 175,115 matched spectra, of which 2,643 proteins were shared by the transcriptome and proteome data.
DEGs related to abdominal hair development of M. javanica
To further understand the genetic basis of abdominal hair development of pangolin, we conducted DEGs analysis before (embryo-hair) and after (adult-hair) occurrence of abdominal hair. In total, 2,248 DEGs were identified between embryo-hair vs adult-hair group; among them, 862 genes were up-regulated in adult-hair group and the remaining 1,386 genes were up-regulated in embryo-hair group (Supplementary Table 2). Based on the annotation of the GO database, the up-regulated DEGs in embryo-hair group were significantly enriched in the component categories related to transporter activity, including transmembrane transporter activity, substrate-specific transporter activity, ion transport, channel activity, etc. However, the up-regulated DEGs in adult-hair group were significantly enriched in categories related to filaments, such as intermediate filaments, intermediate filament cytoskeleton, keratin filaments, etc. (Supplementary Table 3, 4). And its directed acyclic graph of enriched GO terms associated with the keratin filament term, genes related to the supramolecular complex, supramolecular polymer, supramolecular fiber, intermediate filament, intermediate filament cytoskeleton, polymeric cytoskeletal fiber, and cytoskeleton were significantly enriched. The associated genes KRT14, KRT25, KRT32, KRT71, and SYNM may play an important role in the development of abdominal hair in pangolin (Fig. 2a).
For the KEGG enrichment analysis, a total of 161 up-regulated and 117 down-regulated genes in embryo-hair group of 16 unique KEGG pathways related to the hair development process were identified. Among these genes, PI3K-Akt and the ECM-receptor interaction signalling pathway were significantly enriched, and several enriched pathways related to hair development such as focal adhesion, hedgehog, Wnt pathways, etc., were identified (Fig. 2b). It is noteworthy that in embryo-hair group, the up-regulated genes FGF10, VEGFC, FGF16, FGF18, FGFR4, THBS4, EREG, MCL1, and AREG do not only play an important regulatory role in the development of abdominal hair in pangolin, but they also participate in various processes including cellular apoptosis, proliferation, adhesion and attachment, and the inflammatory response. Furthermore, they are involved in epidermal growth, animal organ morphogenesis, and keratinocyte proliferation in GO annotation.
We conducted distinct expressed proteins (DEPs) analysis before and after hair development based on our proteome data. In total, 67 DEPs (embryo-hair vs adult-hair group) were identified, of which 34 proteins were shared with DEGs in the transcriptome (Supplementary Table 5). Among these DEPs, were proteins involved in cell attachment, controlling growth, and intermediate filaments in GO annotation, including 2 up-regulated genes in embryo-hair group: COL11A1 and POSTN, and 5 up-regulated genes in adult-hair group: FCGR2B, PTBP3, DSC1, CDSN, and CSTA.
DEGs related to dorsal scale development of M. javanica
To identify genes that activate the development of pangolin scales, we analysed DEGs before (embryo-scale) and after (adult-scale) the occurrence of dorsal scale. A total of 2,556 DEGs were identified between embryo-scale vs adult-scale group, of which 875 genes were up-regulated in adult-scale group, and the remaining 1,681 genes were up-regulated in embryo-scale group. In the GO enrichment analysis with all DEGs, 7 terms related to the extracellular region and metallopeptidase activity were significantly enriched. In addition, 15 terms related to development and growth were mainly enriched in growth factor activity, regulation of cell growth, multicellular organism development, anatomical structure development, etc. (Fig. 3a). In the KEGG enrichment analysis, genes associated with PI3-Akt, ECM-receptor interaction, Hedgehog, and focal adhesion signalling pathways were significantly enriched (Fig. 3b). In addition, there were 16 signalling pathways related to scale development (Table 2). Notably, FGF1, SGK1, IGFBP3, TEAD4, SERPINE1, and KRAS genes were highly expressed in adult-scale group. We speculate that these genes play an important role in activating the development of pangolin scales.
We conducted DEPs analysis before and after scale development (embryo-scale vs adult-scale group) based on our proteome data. In total, 54 DEPs were identified, of which 17 proteins were shared with DEGs in the transcriptome (Supplementary Table 6). According to the gene description, 3 proteins (COL16A1, RAB21, and MAPRE1) were related to cell spreading, alterations in cell morphology, and cell adhesion, respectively.
DEGs related to skin appendage differentiation of M. javanica
To explore genes controlling the differentiation of skin appendages in pangolin, we compared the gene expression profiles in skin appendages between adult-hair and adult-scale group. Differential expression analysis identified 1,825 up-regulated genes in adult-hair group and 2,486 up-regulated genes in adult-scale group (Supplementary Table 7). According to the gene description, we screened 8 genes related to cell differentiation (MYOD1, PPDPF, COPRS, MYADML2, GDAP1, MYADML2, MAL, and MALL), 7 genes related to cell proliferation (PPDPF, MTCP1, SIPA1L2, SIPA1L3, BTG1, BTG3, and MKI67), 8 genes related to cell apoptosis (AIFM1, AIFM2, BFAR, BCL2, CCAR1, DDIAS, LOC108396282, and LOC108384891), 5 genes related to cell development (NDE1, LMBR1, DRG1, MESDC2, and novel.887), and 17 genes related to keratin (KRT8, 14, 32, 36, 38, etc.). These genes might be indicative of pangolin transcriptome involvement in scale and hair development and differentiation.
Annotation of the M. javanica transcripts with the GO database revealed that the up-regulated DEGs in adult-hair group were significantly enriched in terms related to skin appendage differentiation such as cell adhesion, cell differentiation, biological adhesion, etc. However, the up-regulated DEGs in adult-scale group were significantly enriched in terms related to protein complex, protein catabolic process, transferase activity, etc. In one directed acyclic graph of GO enrichment, the terms cell differentiation, developmental process, anatomical structure development, and multicellular organism development were significantly enriched; the main DEGs involved in these processes included PPDPF, JAG1, CSTA, DLL1, NOTCH1, etc. (Fig. 4a). The KEGG enrichment analysis revealed a total of 35 up-regulated and 79 down-regulated genes in adult-hair group of 20 unique KEGG pathways related to hair development and the differentiation process, including ECM-receptor interaction, proteasome, FoxO, MAPK, NOTCH, P53, HIPPO pathway, etc. (Fig. 4b; Supplementary Table 8). It is noteworthy that the Notch signalling pathway and PPDPF, DLL1, JAG1, VANGL2, WNT9A, LOC108401680, NOTCH1, and NOTCH3 gene were found in both annotation results, and are closely related to the differentiation of keratinocytes and epidermal cells.
According to our proteome data, differential proteome analysis identified 91 distinct significant proteins, including 73 up-regulated proteins in adult-hair group and 18 up-regulated proteins in adult-scale tissues (Fig. 5a). In the GO enrichment analysis with all DEPs, the CDSN gene was enriched in the skin morphogenesis term, which is also involved in biological processes related to skin development and keratinocyte differentiation. JUP, DSC1, CDH1, and TLN2 genes were enriched in cell adhesion and adherens junction terms, which were also involved in cell migration and regulation of cell population proliferation (Fig. 5b). In the KEGG enrichment analysis with all DEPs, 41 proteins in 21 pathways were related to the development and differentiation of skin appendages (Fig. 5c). We also analysed the subcellular localisation of all distinct proteins and found that the proportion of extracellular proteins was the highest (26.5%), followed by cytoplasmic proteins (17.7%) (Fig. 5d).
The mRNA information obtained from the transcriptome was integrated with the protein information identified by the proteome to identify any corresponding relationships. A total of 1,144 IDs were identified by both RNA-seq and proteomics (Fig. 6a). There were 6 distinct proteins corresponding to transcriptome DEGs: CSTA, XDH, PSMC5, KLK10, VCP, and MYOZ1. The correlation between the different multiples of genes (proteins) identified by transcriptome and proteome analysis in the two groups was analysed (Fig. 6b). The Pearson correlation coefficient between mRNA and the corresponding protein was positive (Pearson = 0.461). Thus, we conclude that it is important to determine the protein level to understand phenotypic changes and to not rely solely on mRNA levels.
Immune-related genes in the epidermis of M. javanica
In this study, many immune-related genes were also identified in tissues of skin appendages from different stages of development in pangolin. Through GO enrichment analysis, the immune-related DEGs before and after hair growth (embryo-hair vs. adult-hair group), 34 up-regulated genes and 4 immune-related biological processes were identified as enriched in adult-hair group. These included defence response, inflammatory response, immune system, and immune response, while only 4 immune-related up-regulated genes and two terms were enriched in embryo-hair group (Fig. 7a). The results obtained from KEGG enrichment analysis were similar to those of the GO enrichment analysis. Compared with the embryonic samples (embryo-hair: 15 signalling pathways and 61 immune genes), more immune-related pathways and genes (adult-hair: 17 signalling pathways and 123 immune genes) were enriched in the adult-hair group (Fig. 7b).
Through GO enrichment analysis with all DEGs, before and after scale growth (embryo-scale vs. adult-scale group), we found that 17 up-regulated genes and 3 immune-related biological processes were enriched in adult-scale group. These included regulation of the immune system process, immune system, and immune response, while there was no enrichment of terms related to the immune response in embryo-scale group (Fig. 7c). The results obtained from KEGG enrichment analysis were similar to those of the GO enrichment analysis. Compared with the embryonic samples (embryo-scale: 15 signalling pathways and 61 immune genes), more immune-related pathways and genes (adult-scale: 17 signalling pathways and 123 immune genes) were enriched in the adult-scale tissue (Fig. 7d).
We also compared and analysed the DEGs related to immunity in the dorsal and abdominal skin samples (adult-hair vs. adult-scale group) of adult pangolins. After GO enrichment analysis with all DEGs, both tissue types were found to be enriched in 2 terms related to immunity (immune system process and immune response), while the number of genes related to immunity and enriched in scale-type tissues (4 genes) was higher than that in hair-type tissues (2 genes) (Fig. 8a). In addition, we observed similar KEGG enrichment analysis results. Scale-type tissues had more immune-related pathways and genes than that found in hair-type tissues (Fig. 8b). Scale-type tissues exhibited 16 enriched immune-related pathways and 47 immune-related genes. The number of terms and genes enriched in hair-type tissues was significantly fewer than in scale-type tissues.