Comprehensive transcriptional analysis of pig facial skin development

Animals and tissue collection

In this study, 12 Chenghua sows from four different development periods, including 3-days-old, 90-days-old, 180-days-old, and 3-years-old, were selected as experimental animals, and three pigs from the same litter were considered as the biological replicates per development period. All pigs were divided into four experimental groups according to the development periods, which included the D3 group, D90 group, D180 group and the Y3 group. The piglets were weaned at the age of 28 ± 1 day. A starter diet containing of 18.0% crude protein, 7.0% crude ash and 1.32% lysine was administered to piglets from day 30 to day 45 after weaning. From day 46 to day 179, the pig’s diet consisted of 16.0% crude protein, 9.0% crude ash and 1.0% lysine. From the 180st day, pigs received a diet containing of 14.0% crude protein, 9.0% crude ash and 0.6% lysine. Pigs are allowed free access to food and water under the same conditions. Food and water were withheld from the pigs for 24 h before slaughter. After being transported to the site of slaughter, they were allowed to rest for 2 h and then were humanely slaughtering. To reduce pain, a 10-second sudden shock at 50 V and 90 Hz was used.

The 12 Chenghua sows in the experiment were sourced from the Chengdu Livestock and Poultry Genetic Resources Protection Center in Sichuan Province, China. All animal experimental procedures were approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University (permit number: 20220279).

The facial skin tissue of each pig and the fresh samples were flash-frozen in liquid nitrogen and then stored at −80 °C until RNA was extracted. TRIzol Reagent (Invitrogen, Waltham, MA, USA) was used to extract total RNA, which was subsequently treated with DNase and purified using an RNeasy Mini Kit (Qiagen, Valencia, CA, USA). The quality, concentration, and integrity of RNA were checked using a nanodrop photometer and an Agilent 2100 bioanalyzer.

RNA library construction and sequencing

The lncRNA library includes mRNA and lncRNA. According to the manufacturer’s information, the extracted total RNA was first removed using the MGIEasy rRNA Removal Kit. Then, it was submitted to RNA fragmentation and cDNA synthesis (second-strand cDNA synthesis with dUTP instead of dTTP), followed by end repair, addition of an A residue to the 3′  end and adapter ligation, PCR, circularization and generation of DNA Nano Ball (DNB), and finally sequencing on the DNBSEQ-G400 platform and 150 bp paired-end reads were obtained.

Once the Small RNA library was generated, RNAs of 18–30 nt in length were purified and separated using Polyacrylamide Gel Electrophoresis (PAGE) for 3′  and 5′  linkages. Reverse transcription extension with RT primers with Unique Molecular Identifier (UMI)was use to synthesize cDNA. PCR amplification of cDNA was performed with both 3′  and 5′  linkers linked to highly sensitive polymerase to amplify yield. PCR products in the 110–130 bp range were separated using PAGE. Library quantification, pooling cyclization, and quality inspection of the constructed library were performed. Libraries that passed the quality test were sequenced on the DNBSEQ-G400 and 50 bp single-end reads were obtained.

CircRNAs were constructed using DNase I to digest the DNA fragments present in the total RNA sample; then, they were purified and recycled, and ribosomal RNAs were removed from the total RNA sample using the Ribo-off method. The RNase R reaction system was prepared, linear RNA components were digested, reaction products were purified and recovered, RNA was fragmented at a certain temperature and ion environment, one-stranded cDNA was synthesized using fragmented RNA as a template, and two-stranded cDNA was synthesized with dUTP instead of dTTP. The ends of the double-stranded cDNA were repaired, an “A” residue was added to the 3′  end, the linker was ligated, and the two-stranded cDNA containing “U” residues was digested with UDGase and then submitted to PCR amplification. The constructed library was quality tested, the library products that passed the quality test were cycled, and the circular DNA molecules were copied through rolling rings to form DNA nano balls (DNBs), and finally sequenced by using the DNBSEQ-G400 platform, obtain 150bp paired-end reads.

Quality control and read alignment

After removing the reads containing the adaptor (adaptor pollution), low-quality reads and reads whose N content was greater than 5% from the raw sequencing data, the resulting clean reads were compared to the reference genome and transcriptome (GCF_000003025.6_Sscrofa11.1).

Expression quantification

We used Bowtie2 (Langmead & Salzberg, 2012) to map the lncRNA library clean reads to the reference sequence and then used RSEM (Li & Dewey, 2011) to calculate the expression levels of genes and transcripts. Similarly, we used Bowtie2 to align clean reads to the reference set and other small RNA databases and used UMI for relevant gene expression calculations. Meanwhile, we calculated the expression of circRNA based on the number of back-spliced reads compared to both ends of the circRNA and used two software programs, CIRI (https://sourceforge.net/projects/ciri) and find_circ (https://github.com/marvin-jens/find_circ), prediction. The final back-spliced read number is the average of the two results.

Differential expression analysis

The four developmental stages Chenghua sows had three biological replicates, respectively. The statistical power of this experimental design, calculated in RNASeqPower is 0.84 (https://doi.org/doi:10.18129/B9.bioc.RNASeqPower). A differential expression analysis of RNA was performed using the DEGseq R package. DEGseq (Wang et al., 2010) based on MA-plot (Yang et al., 2002) was used to calculate the differential expression. The p-value was adjusted to the q-value by Benjamini & Hochberg (1995) and Storey & Tibshirani (2003). A Q-value <0.005 and —log2 (fold change)—>1 were set as the thresholds for significantly differential expression.

According to the results of differential gene detection, the R package heatmap was used to perform hierarchical clustering analysis on the union set differential genes (https://cran.r-project.org/web/packages/pheatmap/). When multiple groups of DE miRNAs were clustered at the same time, we performed separate cluster analyses on the intergroup intersection and union of differentially expressed miRNAs.

Based on the GO and KEGG annotation results and the official classification, we performed functional classification and biological pathway analysis of genes derived from different circRNA sources, while enrichment analysis was performed using the hyperfunction in R software. The p-value was then adjusted to the q-value, and typically pathways with q-value <= 0.01 are considered to be significantly enriched.

Target gene prediction of lncRNAs and miRNAs

For differentially expressed lncRNAs, we predicted their target genes by the following steps: calculating the spearman and pearson correlation coefficients between lncRNAs and mRNAs and requiring spearman coefficient >=0.6 and pearson coefficient >=0.6. Then, it was considered the position relationship between lncRNAs and mRNAs. When the lncRNA was located within 20Kb upstream and downstream of the mRNA, it was considered to be cis-regulated. Beyond this range, it would analyze the binding energy of lncRNA and mRNA by RNAplex (v0.2) (Tafer & Hofacker, 2008), and when the binding energy is <-30 and it would be judged as tran-acting mode. RNAhybrid (Kruger & Rehmsmeier, 2006), miRanda (John et al., 2004) and TargetScan (Agarwal et al., 2015) were used to predict target genes of differentially expressed miRNAs.

RNA interaction network construction

The target genes of differentially expressed lncRNAs and circRNAs were predicted by miRanda and RNAhybrid software, and the target genes predicted by the two software programs were selected. Finally, Cytoscape 3.9.1 (http://cytoscape.org/) was used to map the corresponding RNA interaction network according to the predicted lncRNA−miRNA, miRNA−mRNA, and circRNA−miRNA.

Hematoxylin and Eosin (H&E) Staining

According to the manufacturer’s instructions, fixed tissues were embedded in paraffin (Servicebio, Wuhan, China) and cut into 3–4 µm thick sections using a microtome (Servicebio, Wuhan, China). Images were captured with a microscope (Nikon, Japan). All staining assays were performed in triplicate.

Validation Real-Time qPCR (qRT − PCR)

After RNA isolation, cDNA is synthesized by a reverse transcription kit (Takara, Chengdu, China). qRT −PCR was performed with SYBR-Green I nucleic acid dye in the CFX 96 real-time system (BIO-RAD, USA). Primer sequences for 2 mRNAs, 2 lncRNAs, 1 miRNA, and 1 circRNA were designed and synthesized by NCBI and Primer 5 (Table S1). GAPDH was used to normalize the expression levels of mRNAs, lncRNAs and circRNAs, and U6 was used to normalize the expression levels of miRNAs. At least three samples were analyzed for each developmental stage (D3, D90, D180, and Y3), and each sample was analyzed in three independent reactions. The results were statistically analyzed using 2−ΔΔ CT relative quantification.

Descriptive statistics and correlation analysis of skin development in pigs after birth

To understand the development of skin in pigs after birth, we performed RNA-seq on the facial skin of postnatal Chenghua pigs (D3, D90, D180, and Y3). H&E staining of the facial skin tissue of Chenghua pigs in four developmental stages showed that the skin thickness in the slices increased with increasing age (Fig. 1A). In this study, we constructed lncRNA and small RNA libraries. A total of 44,611 genes were detected in the lncRNA library (Table S2), and the average alignment rate of the sample com parison genome was 90.35% (Table S3). A total of 1,769 small RNAs were detected in the small RNA library (Table S4), and the average comparison rate of the sample comparison genome was 90.25% (Table S5). The expression of circRNA is shown in Table S6. The proportion of clean reads Q30 (or Q20) of the filtered RNA was greater than 92% (Tables S7–S9), and the base mass distribution of clean reads showed a low proportion of bases of low quality (Quality<20) (Fig. S1–S3), which indicates good sequencing quality.

Principal component analysis (PCA) was used for the identified mRNA and lncRNA transcripts, and the results showed that the D90 group relative aggregation, D3, and Y3 had a very low correlation, and the similarity of adjacent developmental stages was higher than that of nonadjacent developmental stages (Fig. 1B). According to the expression boxplot, the distribution of miRNA expression levels in the four periods was concentrated, especially at D180 (Fig. 1C). At the same time, the PCA results of miRNA also showed that the correlation between D3 and Y3 was lower than that between D3 and D90 and between D3 and D180, and the similarity of adjacent developmental stages was higher than that of nonadjacent developmental stages. The D90 group was clustered together, which was similar to the PCA results of mRNA and lncRNA (Fig. 1D).

After the quality control of transcriptome sequencing, the RNAs expression profiles at different time points were determined, and mRNA was confirmed to be expressed by 17,998, 17,957, 17,882, and 17,819 genes; miRNA was confirmed to be expressed by 1,447, 1,342, 1,297, and 1,294 genes; lncRNA was confirmed to be expressed by 16,392, 16,662, 16,454, and 16,108 genes; and circRNA was confirmed to be expressed by 6,260, 6,373, 6,243, and 6,043 genes at D3, D90, D180, and Y3, respectively (Fig. 1E).

Identification of differentially expressed mRNAs and noncoding RNAs (ncRNAs)

The molecular mechanism and pathway of pig skin growth and development are related to the expression abundance of specific RNAs. Statistics on the basal expression of mRNAs, lncRNAs, miRNAs and circRNAs at the four developmental stages showed that although there were tens of thousands of RNAs expressed in each stage, the proportion of high-expression genes that truly played a role in skin development was not high (Table 1). Interestingly, the genes with the highest expression for all four types of RNA at the four- time points were basically the same.

In this study, in skin tissue, a total of 242, 428, 625, 54, 323, and 138 upregulated and 220, 446, 453, 54, 114, and 41 downregulated mRNA differentially expressed genes (DEGs) (Fig. 2A); 38, 65, 39, 23, 18, and 20 upregulated and 41, 76, 78, 37, 55, and 41 downregulated lncRNA DEGs (Fig. 2B); 12, 19, 38, two, nine and five upregulated and 23, 21, 80, 2, 25, and 27 downregulated miRNA DEGs (Fig. 2C); and 59, 73, 85, 34, 55, and 43 upregulated and 47, 56, 69, 30, 46, and 31 downregulated circRNA DEGs (Fig. 2D) were detected in D3 vs. D90, D3 vs. D180, D3 vs. Y3, D90 vs. D180, D90 vs. Y3 and D180 vs. Y3 (—log2FC—>1, Q value <0.05). Interestingly, as the time interval between development increases, so does the number of differentially expressed RNAs. For example, D3 vs. Y3 had significantly more differentially expressed RNAs than D3 vs. D90.

The function of the differentially expressed mRNAs

The Upset diagram clearly shows the expression of the differentially expressed mRNAs (DE mRNAs) between the different comparison groups of the four developmental stages. Two DE mRNAs (immunoglobulin superfamily member 10, IGSF10 and Elastin, ELN) were expressed in all comparison groups (Fig. 3A), especially ELN, which was relatively high in each stage, and the expression of the two genes was continuously downregulated (Fig. 3B). The results show that these coexpressed DE mRNAs may have important roles in skin growth and development, especially due to their expression levels and the fact that they are co-expressed.

Heat shock protein 70.2 (HSP70.2) and heat shock protein family B (small) member 1 (HSPB1) were identified as DEGs with upregulated mRNA expression in D3 vs. D180 and D90 vs. Y3, and keratin 31 (KRT31), keratin 34 (KRT34), keratin 33A (KRT33A), collagen type XIV alpha 1 chain (COL14A1), collagen type XXI alpha 1 chain (COL21A1), collagen type VII alpha 1 chain (COL7A1) and ELN were identified as DEGs with downregulated mRNA expression in D3 vs. D180, D3 vs. Y3, D90 vs. D180 and D90 vs. Y3. As shown in the heatmap, the expression levels of these DE mRNAs were more significant (Fig. 3C). According to the KEGG pathway analysis, we found that these genes were mainly enriched in the following pathways: protein digestion and absorption, estrogen signaling pathway, and staphylococcus aureus infection (Fig. 3D). In addition, KEGG analysis of DEGs for D3 vs. D180 found that they were significantly enriched in the skin growth and development-related pathway such as peroxisome proliferators-activated receptors (PPAR) (Fig. 3E), and DEGs from D90 vs. D180 and D180 vs. Y3 were also enriched in this pathway (Fig. S4).

The function of target genes of differentially expressed noncoding RNAs (ncRNAs)

Gene Ontology (GO) analysis was used to analyze the main functions of the differentially expressed miRNAs (DE miRNAs). According to the GO database, we found that the terms enriched across groups of miRNAs had a high degree of overlap, and they were mainly enriched in the following biological processes: cellular process, biological regulation, and regulation of biological process and metabolic process; cellular components: cell, cell part, organelle and membrane; and molecular functions: binding, catalytic activity and molecular function regulator (Fig. 4A, Fig. S5). In addition, there was a high degree of overlap in KEGG pathways between the comparison groups, including cell growth and death, cellular community-eukaryotes, folding, sorting and degradation, development and regeneration, aging and other pathways (Fig. S6). This suggests that these terms may play an important role in the development of skin after birth.

To explore the function of lncRNA-associated mRNA, we selected the two genes (URS0000EF6FAF, URS0001974785) with the highest expression and 1 gene (URS0000EF1D34) with continuous upregulation from all upregulated lncRNAs and selected the top six genes with the highest expression from all downregulated lncRNAs, of which four genes (URS0001952842, URS0000EF04FE, URS0001977D95, and MEG3) were continuously downregulated, and one gene (URS0001961E18) showed no expression at the Y3 stage (Fig. 4B). We predicted the target genes of these nine lncRNAs and obtained 4 associated mRNAs. GO analysis found that the enriched terms were similar to those of miRNA (Fig. 4C).

The function of circRNA is related to the function of host linear transcripts, and we used GO and KEGG analyses to determine the host genes that the circRNAs differentially regulated. Through GO analysis, we found that each comparison group was mainly enriched in biological processes, such as cellular process, biological regulation and regulation of biological process, and developmental process, especially D3 vs. Y3, a comparison of highly separated developmental stages, for which 12 genes were enriched in developmental process (Fig. 4D, Fig. S7). In addition, KEGG analysis found that the differentially expressed circRNAs (DE circRNAs) had common differential signaling pathways, including cellular community-eukaryotes, cell growth and death, transport and catabolism, aging, and development (Fig. S8).

Construction of RNA interaction networks

RNA transcripts communicate through the ceRNA language, and lncRNAs act as sponges for miRNAs to regulate gene expression (Chen et al., 2019; Tay, Rinn & Pandolfi, 2014). In this study, we constructed the lncRNA–miRNA–mRNA coexpression network through Cytoscape, consisting of nine lncRNA nodes, 32 miRNA nodes and 107 mRNA nodes (Fig. 5A). LncRNAs, such as URS0001977D95 and URS0001961E18, target 28 miRNAs respectively, and miRNAs, such as novel-ssc-miR1-5p and novel-ssc-miR107-5p, are targeted by five lncRNAs. In addition, mRNAs, such as cyclin dependent kinase inhibitor 2B (CDKN2B) is associated with cell growth and death, and FERM domain containing 4B (FRMD4B) is associated with corpus callosum agenesis with facial anomalies and cerebellar ataxia.

Studies have shown that circRNAs can act as competitive endogenous RNAs (ceRNAs) to regulate miRNA function (Schorr & Mangone, 2021; Zhang et al., 2021), suggesting that circRNAs and their target miRNAs may be coexpressed in the development of skin tissues. Therefore, we used miRanda and RNAhybrid software to predict the target miRNAs of the circRNA. We identified a total of 13,584 circRNAs and 3,228 miRNAs with targeted binding relationships and constructed an interaction network diagram for the top 20 circRNAs with miRNAs with the greastet correlations (Fig. 5B). However, all these findings require further study.

Validation of RNA-seq data

We randomly selected two mRNAs, two lncRNAs, one miRNA, and one circRNA to validate the whole transcriptome sequencing data by real-time quantitative PCR (qRT–PCR). The q–PCR results were consistent with the RNA-seq data (Fig. 6, Table S10). The expression of mRNA, such as CTSK and ELN, showed significant differences over time in the four groups. In addition, lncRNA (MEG3 and URS0001952842), miRNA (ssc-miR-30b-5p), and circRNA (novel-circ-008475) also showed the same differences. Therefore, the q −PCR results verified the accuracy of the RNA-seq data.