Transcriptional Differences of Coding and Non-Coding Genes Related to the Absence of Melanocyte in Skins of Bama Pig

Skin is the body’s largest organ, and the main function of skin is to protect underlying organs from possible external damage. Melanocytes play an important role in skin pigmentation. The Bama pig has a “two-end-black” phenotype with different coat colors across skin regions, e.g., white skin (without melanocytes) and black skin (with melanocytes), which could be a model to investigate skin-related disorders, specifically loss of melanocytes. Here, we generated expression profiles of mRNAs and long noncoding RNAs in Bama pig skins with different coat colors. In total, 14,900 mRNAs and 7549 lncRNAs were expressed. Overall, 2338 mRNAs/113 lncRNAs with FDR-adjusted p-value ≤ 0.05 were considered to be differentially expressed (DE) mRNAs/lncRNAs, with 1305 down-regulated mRNAs and 1033 up-regulated mRNAs in white skin with|log2(fold change)| > 1. The genes down-regulated in white skin were associated with pigmentation, melanocyte–keratinocyte interaction, and keratin, while up-regulated ones were mainly associated with cellular energy metabolisms. Furthermore, those DE lncRNAs were predicted to be implicated in pigmentation, keratin synthesis and cellular energy metabolism. In general, this study provides insight into the transcriptional difference involved in melanocyte-loss-induced keratinocyte changes and promotes the Bama pig as a biomedical model in skin research.


Introduction
Skin is the largest organ for mammals, exhibiting a complex heterogeneous and multilayered structure and containing various components and more than 10 types of cells [1]. This organ is the main barrier to the external environment and protects underlying organs from trauma and radiation damage.
Much research has been performed to investigate or demonstrate the structural components and physiological mechanisms of skin in humans and mice. For example, pigmentation, which depends on the production of eumelanin and pheomelanin in melanocytes, provides protection for the skin of humans and other animals, as well as directly affecting appearance [2]. Researchers have also put into a 1.5 mL tube with 1 mL of RNAlater (Life Technologies, Beijing, China). Eight tissue samples from three Bama pigs (heart, brain, liver, spleen, mesentery, longissimus dorsi, kidney, lung) were also collected. Finally, the tubes and tissues were immersed in liquid nitrogen and stored at −80 °C. All experimental and sample collection procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the College of Animal Science and Technology of Sichuan Agricultural University, Sichuan, China, under permit No. DKY-S20163629.  White and black skin were separately sampled from the pigs' back and buttocks, and the pie plot shows the comparisons of fibroblast, keratinocyte and melanocyte in two different skins after evaluation by CIBERSORT [13]. Except for fibroblast, keratinocyte and melanocyte varied significantly between the two groups using the Student's t-test (two-tailed). (*: p < 0.05).

Total RNA Extraction, Sequencing, and Read Mapping
A total of 3 µg of RNA (per skin sample) was extracted using Trizol reagent (Life Technologies, Beijing, China), in accordance with the manufacturer's instructions, as was RNA of eight tissues from 3 Bama pigs (heart, brain, liver, spleen, mesentery, longissimus dorsi, kidney, lung). The integrity of RNA was checked by gel electrophoresis on 1.0% agarose gel with GoldView staining using an Agilent 2100 bioanalyzer. A NanoDrop spectrophotometer (Thermo Scientific) was used to measure the RNA concentrations. The RNAs with a ratio of absorbance at 260/280 nm of over 1.8 were selected for further study.
Six samples from three different pigs (named B1, B2, B3, W1, W2, and W3) were selected for sequencing. Approximately 1 µg of total RNA (per sample) and oligo (dT) magnetic heads were used for enriching poly (A) RNAs. The resulting fragments were used as a template for reverse transcription. Random hexamer primers, buffer, dNTPs, DNA polymerase I, and RNase H were used for generating RNA-Seq complementary DNA (cDNA) libraries; next, RNA-seq was conducted following the manufacturer's standard procedures. High-quality strand-specific libraries were sequenced on the HiSeq X platform (Illumina, San Diego, CA, USA) and the bases were called using the software CASAVA v.1.8.2 (Illumina); then, 150-bp paired-end reads were obtained.
High-quality data were controlled by removing poly-N and low-quality reads from the raw data. As a consequence, a total of 44 GB of clean data was acquired. The Q 30 scores and GC content of the clean data were calculated. Clean data were mapped to the pig genome (Sus scrofa 11.1 of release 90 from Ensembl) using TopHat (version 2.1.0) [14]. (Table S1) The RNA-Seq data have been deposited in the NCBI (National Center for Biotechnology Information) Gene Expression Omnibus (GEO) and the accession number is GSE125517.

LncRNA Identification
Mapped reads were assembled by StringTie and then merged with Cuffmerge (part of Cufflinks version 2.2.1) [15]. Then, coding transcripts were filtered by the following steps: (1) remove transcripts of coding gene while comparing to annotated genome by Cuffcompare (part of Cufflinks); (2) comparing with the Pfam-27 database and trimming out transcripts with a p-value < 10 −4 by Hmmscan [16]; (3) comparison with uniref and the nr database and trimming out transcripts with a p-value < 10 −10 by BLASTX (https://blast.ncbi.nlm.nih.gov/); and (4) prediction and calculation of the coding potential of the remaining transcripts by CPC [17]. Transcripts without coding potential were retained [18].

Expression Analysis of mRNA and lncRNA
The mRNA and lncRNA expression level of fragments per kilobase per million mapped reads (FPKM) of each sample was calculated by StringTie (version 1.3.3) [19]; mRNAs with FPKM > 0.5 in at least one sample in at least one group were considered to be expressed and lncRNAs with FPKM > 0.1 in at least one sample in at least one group were considered to be expressed. Then Cuffdiff (part  [14] was applied to detect differentially expressed mRNAs and lncRNAs, and those mRNAs or lncRNAs with adjusted-p values < 0.05 and |log 2 (FC)| > 1 were considered to be differentially expressed.

Clustering and Principal Component Analysis
FPKM values of the six samples (B1, B2, B3, W1, W2, and W3) were used for principal component analysis (PCA) analysis, as well as clustering analysis.

Profiling Melanocyte Proportion with CIBERSORT
The RNA-seq data of melanocytes, keratinocytes and fibroblasts were downloaded from reference [20]. The CIBERSORT [13] was used to estimate the proportion of melanocytes, keratinocytes and fibroblasts within skin.

Construction of lncRNA-mRNA Interaction Network
The R package WGCNA was used to detect interactions between lncRNA and mRNA [21]. The Cytoscape (version 3.2.1) [22] was used to construct the lncRNA-mRNA interaction network.

Functional Enrichment Analysis
Gene ontology (GO) functional enrichment analysis and KEGG pathway functional enrichment analysis were performed with DE genes at the DAVID web server (http://david.abcc.ncifcrf.gov/). To predict the functions of the DE lncRNAs, the differentially expressed genes that were within 100 kb of the lncRNAs or showed a high correlation with the lncRNA were collected and performed functional enrichment analysis as well. The KEGG pathways or GO terms with Benjamini-corrected p-value < 0.05 were significant.
Pearson's correlation coefficients between the expression of DE lncRNAs and the DE mRNAs were calculated with Hmisc (an R package from https://cran.r-project.org/), with the aim of identifying correlations between expression of functional lncRNAs and DE mRNAs [23]. All the analysis strategy is shown in Figure S1.

Validation of Genes and lncRNAs by Real-Time PCR
Twenty-one genes (3 lncRNA included) (Table S2) were selected for the validation and another eight tissues of three Bama pigs were used to investigate the regulations between the lncRNA TCONS_00019024 and that of CYGB (n = 3). To investigate the correlation between expression of TCONS_00019024 and CYGB, RNAs from the skin of six Bama pigs were used for the experiment (n = 6). The genomic DNA in RNA samples was removed by gDNA Eraser (TaKaRa, Shanghai, China) at 42 • C for 5 min. Five micrograms of RNA was reverse-transcribed into cDNA using RT Reagent Kit (TaKaRa, Shanghai, China). The primers for the genes and lncRNAs were designed using Primer 5.0 and tested by NCBI Primer-Blast. The volume of the reaction mixture was 10 µL, with 1 µL of cDNA, 0.5 µL of primers, 5 µL of SYBR (TaKaRa, Shanghai, China), and 3 µL of RNA-free water. The following RT-PCR reaction was performed for all genes and lncRNAs: 95 • C for 3 min; followed by 40 cycles of 95 • C for 10 s and amplification at the optimal temperature for each sample for 30 s; 95 • C for 30 s; and then a melting curve analysis (65 • C to 95 • C). The expression of β-actin was used to correct the gene expression data. The 2 −∆∆CT method was used to analyze the RT-PCR data and calculate relative expression. If a gene was up-regulated in black skin, its expression relative to that in white skin was calculated. If a gene was up-regulated in white skin, its expression relative to that in black skin was also calculated. The t-test was used to test the significance of differences in gene expression.

Expression Profiles of mRNAs and lncRNAs
Approximately 49.03 million raw reads were generated for each sample, and after quality control, 48.81 million clean reads per each were obtained for further analysis. A total of 93-97% of the clean reads were mapped to the pig genome reference (Table S1). Finally, 14,900 mRNAs and 7549 lncRNA (including 82 annotated and 7467 novel lncRNA) were substantially expressed in our samples respectively (Tables S3 and S4). Then we compared the characteristics of lncRNA and mRNA in exon number, transcript length, expression level and coding potential. LncRNAs contain more transcripts with fewer exon numbers (median value = 2), shorter transcript lengths (median value = 588 nt), lower expression levels (mean FPKM value = 1.76) and lower coding potentials compared to mRNA (median of exon number = 8; median of transcript length = 2622 nt; mean FPKM value = 30.72) which were consistent with previous research ( Figure S2A-D) [24][25][26]. Great variations were found between white and black skin as the results shown by principal component analysis (PCA) (Figure 2A,B). The hierarchical cluster of expressed ( Figure S3A,B) or differentially expressed ( Figure S4A,B) mRNAs and lncRNA, and Pearson matrix correlation ( Figure S3C,D) based on expressed mRNA and lncRNA could clearly differentiate white and black skins. Further, as shown in Figure 1, the fibroblast is the major cell type in the skin, and there is a lack of melanocytes in the white skin, which is consistent with previous observations [5]. The above results suggested that our experiment was reliable.  Figure 1, the fibroblast is the major cell type in the skin, and there is a lack of melanocytes in the white skin, which is consistent with previous observations [5]. The above results suggested that our experiment was reliable.

Functional Enrichment Analysis of DE mRNA
A total of 2338 mRNAs were identified to be differentially expressed between two groups (Table  S5). Among these DE mRNAs, 1305 were down-regulated and 1033 up-regulated in white skin ( Figure 3A), including 239 genes down-regulated less than 0.25-fold and 295 genes up-regulated more than 4-fold (marked in Table S5).

Functional Enrichment Analysis of DE mRNA
A total of 2338 mRNAs were identified to be differentially expressed between two groups (Table S5). Among these DE mRNAs, 1305 were down-regulated and 1033 up-regulated in white skin ( Figure 3A), including 239 genes down-regulated less than 0.25-fold and 295 genes up-regulated more than 4-fold (marked in Table S5).

Functional Enrichment Analysis of DE lncRNAs
A total of 113 lncRNA were identified to be differentially expressed, among which 34 were down-regulated and 79 were up-regulated in white skin (Table S7, Figure 4A). Furthermore, 88 DE mRNAs were oriented nearby those DE lncRNAs with 100 kb. These DE mRNAs were found to be enriched in the PI3K-Alt signaling pathway (p = 4.0 × 10 −2 ) which was involved in pigmentation ( Figure 4B, Table S8) [29].

Quantitative Real-Time PCR Validation
A total of 21 genes were selected for the validation. The melanocyte-specific genes (TRPM1, TYRP1, PMEL, MLANA and DCT) were significantly highly expressed in black skin and 21 genes above were indeed differentially expressed ( Figure 5A-C).
RT-PCR indicated that the correlation coefficient between TCONS_00019024 and CYGB was 0.71 in the skin (R 2 = 0.50, p =3.5 × 10 −2 ) ( Figure 5D). We also investigated the relative expression of the TCONS_00019024/CYGB pair in eight other tissues. The results showed that TCONS_00019024 and CYGB presented different expression patterns in the tissues, implying that the correlation of TCONS_00019024 with CYGB may only occur in skin ( Figure 5E). In conclusion, the qRT-PCR experiment validated the discovery of RNA-Seq and supported the reliability of RNA-Seq.

Quantitative Real-Time PCR Validation
A total of 21 genes were selected for the validation. The melanocyte-specific genes (TRPM1, TYRP1, PMEL, MLANA and DCT) were significantly highly expressed in black skin and 21 genes above were indeed differentially expressed ( Figure 5A-C).
RT-PCR indicated that the correlation coefficient between TCONS_00019024 and CYGB was 0.71 in the skin (R 2 = 0.50, p =3.5 × 10 −2 ) ( Figure 5D). We also investigated the relative expression of the TCONS_00019024/CYGB pair in eight other tissues. The results showed that TCONS_00019024 and CYGB presented different expression patterns in the tissues, implying that the correlation of TCONS_00019024 with CYGB may only occur in skin ( Figure 5E). In conclusion, the qRT-PCR experiment validated the discovery of RNA-Seq and supported the reliability of RNA-Seq.

Discussion and Conclusions
Skin is both important and the largest organ of the body, and provides protection for underlying organs. Melanocyte inhabits the skin and its major function is the pigmentation of skin and hair. Loss of melanocyte can cause pigmentation deficiency. Mechanisms of melanocyte deficiency are still not fully understood and transcriptome profiles studies related to the lack of melanocytes are insufficient. Even the mechanism associated with the cross-talk between melanocytes and other types of cells in the skin is still not well-understood in the pig. Therefore, we performed transcriptomic analysis of Bama pig skin to reveal the mechanism potentially associated with the loss of melanocytes or interaction networks related to melanocytes in the pig.
Even with a great similarity between white and black skins, the Pearson correlation coefficient was 0.92 ( Figure S3C-D) which was higher than that between skin and other vascular tissues [57]. Furthermore, various differences such as skin thickness, PH [58,59], appearance, cellular energy metabolism, and melanocyte-keratinocyte interaction were found. Only based on these differences could we interpret the discrepancies in the synthesis of melanin, pigmentation, and coat color gene expression profiles.
Melanocytes and keratinocytes are important for their interaction in the skin, and the absence of melanocytes could result in the disappearance of this interaction and may impact on the function of keratinocytes. For example, melanosome transfer, based on calcium and plasma membrane, is an important process of keratinocyte-melanocyte interaction. [3,4]. Interestingly, many down-regulated genes in white skin, such as ITSN1, BTK, FYN and CCR1, were associated with calcium release of the plasma membrane ( Figure 3C, Table S5). These results could provide evidence for the disappearance of the melanocyte-keratinocyte interaction in white skin. Moreover, regarding keratinocytes, proliferation could be stimulated by melanocytes in the procedure of accelerating wound closure as well as fibroblasts [60][61][62]. In comparison with white skin, the black skin had more melanocytes. The black skin may be more prone to hypertrophic scars than the white skin. Previous wound healing experiments also demonstrated that the skin of the Duroc pig (with melanocytes) is more prone to hypertrophic scars than the skin of the Hampshire pig (with melanocyte deficiency) [63][64][65]. Therefore, the black skin of the Bama pig may be more suitable for hypertrophic scar research.
In this study, we identified two DE lncRNAs (TCONS_00077733, and TCONS_00060772) which might be involved in keratin synthesis and cellular energy metabolism. The black skin was different from the white skin in the expression of keratin genes and genes related to cellular energy metabolism. lncRNA was shown to play critical roles in the skin, such as melanin synthesis in goats [66], psoriasis in humans [26] and expression of keratin and promotion of melanogenesis in mice [67]. In addition, lncRNA H19 regulates Dsg1 expression and the consequent keratinocyte differentiation through acting as an endogenous "sponge" of miR-130b-3p in humans [68]. Thus, we deduced that TCONS_00077733 and TCONS_00060772 might be one of causes responsible for the difference between black and white skin, but more evidence is needed to support this assumption. Taken together, the results of this study provide evidence of the interactions involved in pigmentation. This research provides insight into the complex mechanisms associated with the "two-end-black" phenotype.
In conclusion, we had the opportunity to work on the same individual to systematically reveal transcriptome differences. Our results suggested that the loss of melanocytes could contribute to the expression of melanogenesis genes, and lack of melanocytes might be one of the causes of the white skin differing from the black skin in keratinocyte function and the development of hypertrophic scars. In addition, the present research implies that lncRNAs play roles in skin. Our research could promote the application of the Bama pig in research on melanocyte deficiency.