Comprehensive Profiles of mRNAs and miRNAs Reveal Molecular Characteristics of Multiple Organ Physiologies and Development in Pigs

The pig (Sus scrofa) is not only an important livestock animal but also widely used as a biomedical model. However, the understanding of the molecular characteristics of organs and of the developmental skeletal muscle of the pig is severely limited. Here, we performed a comprehensive transcriptome profiling of mRNAs and miRNAs across nine tissues and three skeletal muscle developmental stages in the Guizhou miniature pig. The reproductive organs (ovary and testis) had greater transcriptome complexity and activity than other tissues, and the highest transcriptome similarity was between skeletal muscle and heart (R = 0.79). We identified 1,819 mRNAs and 96 miRNAs to be tissue-specific in nine organs. Testis had the largest number of tissue-specific mRNAs (992) and miRNAs (40). Only 15 genes and two miRNAs were specifically expressed in skeletal muscle and fat, respectively. During postnatal skeletal muscle development, the mRNAs associated with focal adhesion, Notch signaling, protein digestion, and absorption pathways were up-regulated from D0 to D30 and then down-regulated from D30 and D240, while genes with opposing expression patterns were significantly enriched in the oxidative phosphorylation and proteasome pathways. The miRNAs mainly regulated genes associated with insulin, Wnt, fatty acid biosynthesis, Notch, MAPK, TGF-beta, insulin secretion, ECM–receptor interaction, focal adhesion, and calcium signaling pathways. We also identified 37 new miRNA–mRNA interaction pairs involved in skeletal muscle development. Overall, our data not only provide a rich resource for understanding pig organ physiology and development but also aid the study of the molecular functions of mRNA and miRNA in mammals.


INTRODUCTION
Intensive transcriptome sequencing is increasingly being used to study the mechanisms of organ physiology and development, with a goal of understanding the genetic expressions of tissuespecific diseases (Lage et al., 2008;Koh et al., 2014). Since the expressions, and thus the functions, of many genes vary between tissue types, the study of tissue-specific genetic expression describes those transcriptional variations and provides insights into the underlying genetic mechanisms in each specific tissue type. Thus, obviously the construction of comprehensive tissue-specific transcriptome profiles for both humans and model organisms is of great importance. Indeed, tissue-specific expression of both protein-coding and non-coding RNAs in both humans and mice has been well studied (Roux et al., 2012;Szabo et al., 2015;Zeng et al., 2016;Iwakiri et al., 2017) and the results have been used to elucidate organ physiologies and diseases. For example, the identification of tissue-specific genes was used to build a tissue-specific gene database for human cancers (Kim et al., 2018). However, while many studies have accumulated both mouse and human tissue-specific transcriptome data (Zhang et al., 2015a), RNA-seq transcriptome analyses across tissues and developmental stages of other mammals are still relatively scarce.
Wild Sus scrofa (pig) was domesticated approximately 9,000 years ago and has become one of humankind's most important livestock animals (Giuffra et al., 2000). Because of its similarity with humans in body size, lifespan, anatomy, and other distinct physiological characteristics, the pig has become a model in many disciplines of biomedical research, such as pharmacology, obesity, oncology, cardiology, and many more (Groenen et al., 2012). While pig mRNA and miRNA profiles have been obtained for several tissues, such as adipose, liver, skeletal muscle (Hou et al., 2012;Tang et al., 2015), and testis (Zhang et al., 2015b), those analyses focused mainly on a single organ or developmental stage. Systematic studies of miRNA and mRNA spatiotemporal expression patterns and their interactions in multiple organs and developmental stages are still needed to further understand their physiological functions and development, which would aid both biomedical research and pig husbandry. Interactions between miRNAs and their target mRNAs play important roles in regulating various biological processes (Bai et al., 2015), so identification of those interactions provides insights into the different mechanisms at work in each tissue type (Neville et al., 2011;Wang et al., 2012b).
In this study, we used high-throughput transcriptome sequencing to comprehensively explore S. scrofa mRNA and miRNA profiles in nine different tissues and three developmental stages of skeletal muscles. First, we systematically analyzed expression characteristics of protein coding genes (PCGs) and miRNAs in nine tissues, identifying the tissue-specific and -associated PCGs/miRNAs and then exploring the interactions of miRNAs with their target mRNAs in nine organs. Finally, we detected a set of miRNA-mRNA interaction pairs potentially associated with postnatal skeletal muscle development. Overall, this study provides a comprehensive profile of mRNAs and miRNAs in multiple organs and developmental stages in the pig and provides meaningful insights into tissue-specific metabolic regulation at the RNA level.

Animals and Organ Collection
In this study, we collected nine tissues (fat, heart, kidney, liver, lung, skeletal muscle, ovary, spleen, and testis) from the Guizhou miniature pig, one of the most primitive pig breeds in China and a source of high quality meat, at 240 days of age and two additional skeletal muscle samples at postnatal 0 and 30 days. These samples were collected from three biological individuals at each postnatal date (0, 30, and 240 days). All samples were rapidly isolated and immediately frozen in liquid nitrogen. All animal procedures were performed according to protocols approved by the Biological Studies Animal Care and Use Committee in Beijing Province, China.

Isolation of Total RNA and Construction of RNA-seq Libraries
We extracted total RNA from various tissues at least three times, mixing RNA samples from each tissue type into one group per type, and small RNA library for each tissue group was produced. Polyacrylamide electrophoresis gel was used to purify the fragments of 18-30 nt, and then these fragments were ligated to adaptors on both 5′ and 3′ ends. After reversetranscription amplification, the PCR products in length of 90-bp were isolated from 4 % agarose gels, and then sequenced on the Illumina HiSeq 2500 platform. The RNA-seq data for mRNA were deposited in the Gene Expression Omnibus (accession codes GSE73763) as our previous reports Liang et al., 2017), and the reliability of transcriptome data has was verified by qRT-PCR in the study by Tang et al. (2017).

RNA-seq Data Analysis
First, using custom scripts, we trimmed adapters from all RNA sequencing data and then mapped the processed reads from each sample to the S. scrofa reference genome (v10.2) using TopHat2 (Kim et al., 2013) (v2.0.12) with fr-frststrand and the following parameters: mate-inner-dist 20, mate-std-dev 50, microexon-search segment-length 25, and segment-mismatches 2. Alignment results from each sample were then processed using Cufflinks (v2.2.1) with known annotations for transcript assembly (fr-firststrand and min-frags-per-transfrag 3) and then the consensus transcriptome was merged. Finally, we used HTseq-count (v0.6.1) required strand-specific counting (Anders et al., 2015) to quantify genes and transcripts. We then calculated the reads per kilobase million (RPKM), counted on read pairs in cases of paired ends, for the PCGs.

miRNA-seq Data Analysis
Raw sequencing reads were obtained after removing reads without a 3′ primer or insert tag, reads with polyA or 5′ primer contaminants, and reads shorter than 18 nt or longer than 32 bp. Then, the clean reads that could be annotated and aligned to rRNAs, snoRNAs, and tRNAs in the Rfam database (http://rfam. xfam.org) (Nawrocki et al., 2015) were discarded. We mapped the remaining reads to the S. scrofa reference genome (v10.2) using miRDeep2(v2.0.0.8) software (Friedlander et al., 2012). The sequences of known mature miRNAs and their precursors were downloaded from miRBase (http://www.mirbase.org) (Kozomara and Griffiths-Jones, 2014), and the expression level of each miRNA was normalized using the transcripts per kilobase million (TPM) method.

Gene Expression Analyses
In our analyses, all PCG-miRNAs that had an RPKM or TPM greater than 0.1 in at least one sample were considered to be expressed. In addition, Pearson correlation coefficients were calculated to examine the similarities and correlations of mRNA expression in different samples. In our analyses, universally expressed genes are tissue-conserved expressed genes whose RPKM values are greater than 10 in every tissue. For miRNAs, two criteria were used to define universally expressed miRNAs across different tissues: 1) the TPM value in every tissue was more than 1, and 2) the coefficient of variation across all tissues was less than 0.5.
The tissue-associated genes for any given tissue were identified according to a previous study with the Z-score cutoff ≥ 1.5 and RPKM ≥ 1 (Li et al., 2014). We identified tissue-specific genes as having RPKM or TPM ≥ 10, with expression levels in a given tissue being greater than 10-fold higher than the mean expression value of any other tissues.

Co-Expression Network Analysis
We used RNA libraries of nine different tissue types, with all samples collected at postnatal day 240, for network construction. Based on the mRNA expression matrix, we constructed a weighted co-expression network using BioLayout Express (3D) with a Pearson correlation threshold cutoff ≥ 0.90 and a Markov clustering algorithm of 2.2.

miRNA-mRNA Interaction Analyses
miRNA-host gene co-expression pairs were identified based on the following pipeline: 1) the miRNA coordinate overlapped 100% of the protein-coding gene; 2) the host gene and the miRNA are transcribed from the same strand of DNA; 3) the miRNA does not have extra copies in other parts of the genome, since the transcription of each copy of the miRNA gene could be regulated by different mechanisms that would confound the results of our analyses; 4) the intragenic miRNAs and the host genes must be expressed in at least five tissues (TPM or RPKM ≥ 0.1); and 5) significant correlations were identified as p < 0.05 and r > 0.6 (Kaminski et al., 2013;Lin et al., 2016).
mRNA and miRNA pairs were subjected to Pearson correlation analysis and those pairs with r < −0.5 were chosen for further investigation. Then, we used the RNAhybrid (v2.1.2) (Kruger and Rehmsmeier, 2006) and TargetScan algorithms (Riffo-Campos et al., 2016) to detect whether the 3′-untranslated region (3′UTR) of the mRNA in each pair matched the seed region of the corresponding miRNA. The pairs that satisfied those two conditions were used in both KEGG pathway analysis and Gene Ontology (GO) enrichment analysis to further investigate the biological processes and functions associated with those negative correlations. These analyses were based on human annotation using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) web server (http://david.abcc.ncifcrf.gov/) with the EASE value set to 0.05 (Huang et al., 2007;Huang Da et al., 2009).

Differential Expression Analysis
We used R (v3.2.0) DESeq2 package (Love et al., 2014) for differential expression analyses of PCGs and miRNAs in skeletal muscle at three developmental (D) stages (D0, D30, and D240 postnatal days). Analysis was performed between muscle_D0 and muscle_D30 and between muscle_D30 and muscle_D240. In this study, significant differentially expressed genes (DEGs) met the criteria log2-FC ≥ 1 and FDR < 0.05, where FC is fold change and FDR is false discovery rate, and those acceptable DEGs were further examined with GO and KEGG analyses, as above.

Vector Construction, Cell Culture, and Dual Luciferase Reporter Assay
Two miRNA-mRNA interaction pairs (ACTN4/ssc-miR-133a-3p and Prox1/ssc-miR-338) were randomly selected to verify the expression profiles of miRNA and its target mRNA. Through the PCR method, the 3′UTR fragments (3′UTR-wt) flanking miRNA binding sites of these two genes were amplified and then cloned into pmirGLO Dual-Luciferase Vector through the Homologous Recombination Kit (Qingke, China). The mutant types of these two genes with the 3′UTR region (3′UTR-edt) were made by the Homologous Recombination Kit (Qingke, China) and confirmed by sequencing. Primer sequences are listed in Table 1.

Overview of mRNA and miRNA Profiling
To systematically investigate genome-wide expression profiles of S. scrofa PCGs and miRNAs, we first performed RNA-seq and small RNA-seq on nine tissues, as well as on skeletal muscle tissue from three developmental stages of Guizhou miniature pigs. A total of 824,887,548 reads were obtained for RNA-seq analysis, as we described in our previous study . We then measured the expression abundance of PCGs in each tissue ( Figure 1A), detecting 18,576 expressed PCGs (RPKM > 0.1) representing 85.97% of annotated PCGs in the porcine reference genome, and 46.14% of these PCGs were constitutively expressed through all selected tissues. The numbers of expressed PCGs ranged from 10,845 to 15,552 in the nine different tissues ( Table 2). The smallest numbers of expressed PCGs were in skeletal muscle and the largest numbers were in the reproductive tissues (ovary and testis). Also, the distribution of genes with high RPKM (>10) values was larger in the reproductive tissues than in the other tissues ( Figure 1B). These findings indicated that transcriptome complexity and activity in the reproductive system were higher than those in other tissues. We next examined the similarities and correlations between tissues based on global expression profiling. Clustering analysis suggested that the highest transcriptome similarity was shown between skeletal muscle and heart (Pearson correlation, R = 0.79), whereas the testis and liver showed the lowest expression correlation (Pearson correlation, R = 0.41) ( Figure 1C).
Small RNA analysis revealed 131,789,681 high-quality clean reads, accounting for 93.64% of the total reads. Analysis of the size distribution of all reads showed that the major class of reads peaked at 22-23 nt within most of the libraries. However, in liver and testis, the majority of clean reads were at 28-30 nt, followed by 22 nt, thus implying that Piwi-interacting RNAs were enriched in those two tissues. For all 11 libraries, 99,812,523 (75.75%) of the clean reads were mapped to the porcine reference genome. The composition of each RNA library ( Figure 1D) shows our findings that 319 known and 442 novel miRNAs (TPM > 0.1) were expressed in all the tissues we investigated. The number of expressed miRNAs in these libraries ranged from 475 to 594, including 206 to 292 novel miRNAs and 269 to 302 known miRNAs, respectively ( Table 3). The miRNAs with high RPKM  (>10) values had a larger distribution in testis than in other tissues ( Figure 1E).

Universally and Specifically Expressed mRNAs and miRNAs Across Tissues
Focusing on universally expressed mRNAs and miRNAs, we found that 209 mRNAs (Table S1), representing tissue-conserved expressed genes whose RPKM values were greater than 10 in all tissues, were abundantly and stably expressed. This dataset included some well-known housekeeping genes such as GAPDH, ACTB, RPS18, B2M, RPL4, RPL37, and RPL38. According to GO analysis, these genes were significantly enriched in the areas of translation, peptide biosynthetic, and amide biosynthetic ( Table S2), indicating that they have important roles in maintaining essential basal cellular functions. Additionally, we identified 43 universally expressed miRNAs ( Table S3). One of those miRNAs, miR-16, is most likely an important biomarker for several diseases including lung cancer, rheumatoid arthritis, and sepsis (Wang et al., 2012a;Sromek et al., 2017;Dunaeva et al., 2018) in humans; is abundantly expressed in all tissues; and has been used as a control in several systems, including animal models. These miRNAs may serve as candidate reference to normalize miRNA expression across tissues. In order to capture the functional differences in gene expression between tissues, we next analyzed tissue-specific mRNA and miRNA. As in a previous study , the genes whose abundance in one tissue was more than fourfold the mean expression value of that in other tissues were defined as tissue-specific genes. While for some genes the expression levels were very low, we defined the tissue-specific mRNAs and miRNAs as the genes whose abundance in one tissue was more than 10-fold the mean expression value in other tissues, with RPKM or TPM ≥10. Finally, we identified 1,819 tissuespecific PCGs, with the number ranging from 15 to 992 for a given tissue (Table S4). Testis had the largest number of tissue-specific genes (54.5%, 992/1819), compared with other tissues. In contrast, only 15 genes were specifically expressed in skeletal muscle. We found that the expression of GAPDH was significantly higher in skeletal muscle than in other tissues, a result consistent with a previous mouse study (Fortes et al., 2016). Meanwhile, we identified 96 tissue-specific miRNAs (including 48 novel and 48 known miRNAs). The number of miRNAs ranged from 2 to 40 for a given tissue. The largest numbers of tissue-specific miRNAs were found in testis and only two tissue-specific miRNAs were found in fat ( Table 4). The well-known myomiRs (miRNA-1, miR-133a/b, and miR-206) were specifically expressed in skeletal muscle. In testis, the top two tissue-specific miRNAs, miR-34c and miR-202-5p, were shown to possess an important influence in spermatogenesis regulation (Dabaja et al., 2015;Wang et al., 2018b). For fat, only one known miRNA (miR-224), which was reported to have an important role in adipogenesis development, was specifically expressed (Peng et al., 2013).

Tissue-Associated mRNAs Capture the Structure and Functional Features of Different Tissues
A weighted and undirected co-expression network analysis was performed to understand interactions between PCGs, and it generated 229 distinct clusters containing 13,894 nodes (Figure 2A). Cluster names were based on the tissues in which the genes were expressed the most. The three largest clusters (including 4,997 nodes) were groups of highly expressed genes in ovary and testis. These results indicated that the majority of PCGs showed a tissue-restricted expression pattern. Thus, we tried to capture the basic characteristics of gene expression for each tissue by using ''tissue-associated genes," genes that were highly expressed in one tissue relative to other tissues (see    Materials and Methods). The number of associated PCGs and miRNAs (Tables S5, S6) for a given tissue ranged from 287 to 5,606 and from 28 to 132, respectively, suggesting that the tissues with higher transcriptional activities, such as testis, had more associated genes (Figures 2B, C). GO analysis, based on tissue-associated PCGs, revealed physiological features for each organ ( Figure 2D). For example, GO terms for muscle tissue development were significantly enriched in skeletal muscle and heart, while the genes associated with spermatogenesis, lipid metabolism, and immune response were enriched in testis, adipose, and spleen, respectively. Also, we observed common GO terms shared in different tissues. For instance, the GO term related to the cell cycle and metabolic processes was obviously enriched in both ovary and testis, confirming that the reproductive tissues were highly proliferative. The GO terms for small-molecule catabolic and other metabolic processes were markedly enriched in liver and kidney. The genes associated with stimulus response and signal transduction were significantly enriched in lung and spleen. In summary, GO terms of tissue-associated genes agreed with the physiologies of the corresponding organs (Table S7).

Differentially Expressed mRNA and miRNA During Skeletal Muscle Development
To understand postnatal skeletal muscle development, we assessed the differentially expressed PCGs and miRNAs in skeletal muscle across three developmental stages at 0, 30, and 240 days after birth (D0, D30, and D240, respectively). Between D0 and D30, we detected 1,515 DEGs (Table S8), including 911 up-regulated and 604 down-regulated genes, respectively ( Figure 3A). GO analysis (Table S9) suggested that the up-regulated genes were involved mainly in vasculature development, the intracellular signaling cascade, blood vessel development, and enzyme linked receptor protein signaling pathway (Figure 4A), and the down-regulated genes were associated mainly with translation, ribonucleoprotein complex biogenesis, and RNA and ncRNA processing ( Figure 4B). Between D30 and D240, we identified 1,011 DEGs (Table S10) including 338 up-regulated and 673 down-regulated genes ( Figure 3B). According to GO analysis (Table S9), up-regulated genes were involved mainly in protein catabolic process, modification-dependent macromolecules, and modification-dependent protein catabolic process ( Figure 4C). The down-regulated genes were involved mainly in cell adhesion, biological adhesion, and skeletal systems development ( Figure 4D). These findings indicated that proliferative cell activity decreased, while cellular metabolic ability increased with age during postnatal skeletal muscle development and growth. Subsequently, we focused on a series of dynamic expression patterns exhibited during skeletal muscle development. Venn diagram (Figure 5) analysis for DEGs demonstrated that the greatest overlap (276 genes) occurred between both up-regulated genes in group D0 versus D30 and down-regulated genes in group D30 versus D240 ( Table S11). The largest cluster of overlapping genes were significantly associated with vasculature, the cardiovascular and circulatory systems, and blood vessel development, and obviously gathered in the focal adhesion, Notch signaling, and protein digestion and absorption pathways (Tables S12, S13). The second largest overlap cluster contained 94 genes and were present between the down-regulated genes in group D0 versus D30 and the up-regulated genes in group D30 versus D240 (Table S11). These overlapping genes functioned in the ATP metabolic process, purine ribonucleoside triphosphate metabolic process, and ribonucleoside triphosphate metabolic process and were significantly enriched in KEGG pathways for oxidative phosphorylation, proteasome, and Parkinson's disease (Tables S12, S13). However, there were only eight genes up-regulated and seven genes down-regulated throughout D0 to D240 (Table S11).
These miRNAs also played an important role in skeletal muscle development. In this study, we identified 70 and 85 differentially expressed miRNAs in groups D0 versus D30 and D30 versus D240, respectively (Tables S14, S15). Between D0 and D30, 56 miRNAs were up-regulated and 14 were down-regulated ( Figure 3C). Functional analysis suggested that the genes targeted by down-regulated miRNAs were significantly enriched in the insulin, Wnt, and Notch signaling pathways and in fatty acid biosynthesis, while the target genes for up-regulated miRNAs associated mainly with the MAPK and TGF-beta signaling, insulin secretion, and ECM-receptor interaction pathways. In groups D30 versus D240, we detected 21 up-regulated and 64 down-regulated miRNAs, respectively ( Figure 3D). The targets for the up-regulated miRNAs were significantly involved in the MAPK, focal adhesion, and insulin signaling pathways, while the targets for the down-regulated miRNAs were enriched mainly in the MAPK, TGF-beta, Notch, calcium, Wnt, insulin secretion, and focal adhesion pathways.

miRNA-mRNA Interaction Network Associated With Skeletal Muscle Development
miRNAs can affect gene expression by inhibiting protein translation or by causing mRNA degradation (Tang et al., 2015). When we evaluated the expression relationships between miRNA and mRNA using Pearson correlations, we detected 253,057 miRNA-mRNA interactions that were negatively correlated (r < −0.5) and 2,194 pairs (1,605 mRNAs and 263 miRNAs) with binding sites for miRNAs at mRNA 3′UTRs. GO enrichment analysis suggested that these miRNA-mRNA interactions associated mainly with protein catabolic processes, muscle cell differentiation, and muscle organ development. KEGG analysis revealed that the interactions were significantly enriched in pathways for the citrate cycle, axon guidance, purine and pyruvate metabolisms, as well as the gonadotropin-releasing hormone, MAPK, adipocytokine, and insulin signaling pathways.

Dual Luciferase Reporter Assay Validated the Interaction Between miRNAs and Their Target Genes
We randomly selected two miRNA-mRNA pairs (ACTN4/ssc-miR-133a-3p and Prox1/ssc-miR-338), which showed significant FIGURE 5 | Co-expression pattern analysis for genes during skeletal muscle development. August 2019 | Volume 10 | Article 756 Frontiers in Genetics | www.frontiersin.org negative expression correlations (r < −0.5) through all three skeletal muscle developmental stages, to verify whether the interaction between them were really exit. The dual luciferase reporter assay successfully validated the interaction between miRNA and their target genes. As shown in Figure 6, by binding to the 3′UTR region, all of the two miRNAs could markedly decrease the luciferase activity of the wild-type target genes (3′UTR-wt), while for the mutant type (3′UTR-mt), this repression was relieved. These results further confirmed the interaction between miRNA-mRNA pairs, which we discovered in the skeletal muscle development.

DISCUSSION
Transcriptome profiling is a good way to understand physiological functions and developmental regulations of organ tissues in plants and animals (Mcloughlin et al., 2014;Santos et al., 2014). The pig is not only an important livestock animal but also an important model organism in biomedical research. Therefore, a comprehensive atlas of gene expression for S. scrofa tissues and developmental stages is essential for both breeding and biomedical research . Here, we used RNA-seq to profile mRNAs and miRNAs found in nine different organ tissues and in three developmental stages of the Guizhou miniature pig, a breed widely used as a model organism in biomedical research. In all, the 18,576 PCGs we detected represent 85.97% of the annotated PCGs in the porcine reference genome. The numbers of PCGs ranged from 10,845 to 15,552 in different tissues and developmental stages, respectively, but only 46.14% of PCGs were constitutively expressed in all nine tissues. This indicates that PCG expressions are greatly temporally and spatially specific. Interestingly, the reproductive organs (ovary and testis) harbored the most complex transcriptomes, a finding consistent with those in humans and rats in which the testis and ovary also expressed the most genes (Yu et al., 2014).
We further analyzed the tissue-associated, tissue-specific, and universally expressed mRNA and miRNA across different tissues and showed that tissue-specific and -associated mRNAs and miRNAs are needed to maintain specific functions in a given tissue type. In this study, the numbers of tissue-specific and tissue-associated PCGs ranged from 15 to 992 and from 287 to 5,606, respectively. High variances reflected the differences in cell homogeneity and activity between different tissues. Tissues with higher transcriptional activities, such as testis and ovary, contained more associated and specific genes than other tissue types. The complex transcript of the Guizhou miniature pig testis was similar to that of other pig varieties. For instance, during the testis development of the Shaziling pig (a Chinese indigenous breed), 8,343 DEGs were identified and more than 50,000 miRNA-mRNA interaction sites were predicted (Ran et al., 2015). Additionally, many of these tissue-associated and tissuespecific genes are well correlated with physiological functions of each organ. For example, the testis-associated gene PAK2 (p21activated kinase 2) was reported to play a crucial regulatory role in porcine spermatogenesis apoptosis and when the PAK2 gene was knocked down by related siRNA, the mitotic activity for Sertoli cells was significantly repressed (Ran et al., 2018b). Interestingly, while PAK2 is one target gene of miR-26a, another miR-26a target gene, ULK2, was also a testis-associated gene in our study. ULK2, when knocked down, will inhibit swine Sertoli cell autophagy (Ran et al., 2018a). Meanwhile the testis-specific genes identified in this study such as SPEM1, TNP1, PRM1, DAZL (Hashemi et al., 2018), and CABYR (Shen et al., 2019) were also published as specific expression in human or mouse. As we know, tissue-specific functions are a result of specific expression and regulation of genes across an organism's lifespan. The transcriptome's degree of correlation suggests both similar and different biological functions between tissues. These data aid in the understanding of organ physiologies and molecular functions of genes in mammals. Skeletal muscle is an important organ for maintaining movement and energy metabolism in animals (Liu et al., 2018), and pig skeletal muscle is a protein resource for humans Yang et al., 2017). Thus, a systematic study of skeletal muscle development is essential to improving animal breeding as well as aiding biomedical research. Many studies have suggested that the PCGs, miRNAs, and the interactions between them are most important for cellular regulatory processes . Nielsen et al. (2010) analyzed the miRNA in pig longissimus dorsi by using deep sequencing, and they found that highly expressed miRNAs were involved in skeletal muscle development and regeneration (Nielsen et al., 2010). However, the understanding of skeletal muscle development based on a comprehensive profiling of mRNAs and miRNAs had been largely unclear. To fill that void, we carried out RNA-seq and small RNA-seq analysis on pig skeletal muscle at 0, 30, and 240 days after birth. In the D0 versus D30 and D30 versus D240 groups, 1,515 and 1,011 mRNAs, respectively, were differentially expressed. Functional analysis suggested the presence of significant differences in physiological characteristics at different developmental stages. Between D0 and D30, for example, genes functionally associated with translation, ribonucleoprotein complex biogenesis, and RNA and ncRNA processing were down-regulated, and in the D30 versus D240 groups, down-regulating genes were obviously involved in cell and biological adhesion and in skeletal system development.In the D0 versus D30 and D30 versus D240 groups, we detected 70 and 85 differentially expressed miRNAs, respectively. The miRNAs regulated biological processes by binding mRNAs at the 3′UTR. In the current study, 2194 negatively correlated (r < −0.5) miRNA-mRNA interaction pairs with binding sites for miRNAs at mRNA 3′UTRs were predicted. Of those pairs, 37 new miRNA-mRNA interaction pairs were associated with muscle cell differentiation and muscle organ development and were negatively correlated (r < −0.5) in the D0, D30, and D240 groups. Most of the predicted target mRNA in these pairs were reported to function in muscle. For instance, SRF and MBNL1 (serum response factor and muscleblind-like splicing regulator 1) genes are reported to regulate muscle atrophy in mice (Collard et al., 2014), AMOT (the angiomotin gene) may influence human aortic smooth muscle cell migration (Wang et al., 2018a), and QKI (the protein quaking gene) regulates smooth muscle cell differentiation . Notable miRNAs that we found include ssc-miR-744, which is reported to significantly up-regulate in muscles after ischemia-reperfusion injury ; ssc-miR-195, which induces postnatal quiescence of skeletal muscle stem cells (Sato et al., 2014); and ssc-miR-423-3p and ssc-miR-133a-3p, which each showed high correlations with mouse skeletal muscle C2C12 myoblast differentiation (Siengdee et al., 2015;Wang et al., 2018c). By targeting PCGs, miRNAs play important roles in regulating the complex processes of muscle development. Analysis of miRNA and mRNA expression profiles together was an effective way to minimize false-positive rates in miRNA-mRNA interaction pair predictions. In order to discover more miRNA-mRNA interaction pairs, we first analyzed the transcriptomes of both miRNA and mRNA in nine different tissues and then validated the miRNA-mRNA interaction pairs associated with muscle development in three different muscle development stages. For further verification, we also randomly selected two miRNA-mRNA interaction pairs to verify the expression profiles of each miRNA and its target mRNA by using a dual luciferase reporter assay. The results showed that, by binding to the 3′UTR region, miRNA could markedly decrease the luciferase activity of the target genes. Hence, the miRNA-mRNA interaction pairs predicted in this study most likely participate in the regulation of muscle development. The genes in these two pairs were ACTN4, a transcriptional regulator of myocyte enhancer factor that is associated with skeletal muscle differentiation , and Prox1, an essential gene for satellite cell differentiation and muscle fiber-type regulation (Kivela et al., 2016). They were paired, respectively, with important muscle-associated miRNAs ssc-miR-133a-3p (Wang et al., 2018c) and ssc-miR-338 (Mcdaneld et al., 2009). Although most of these miRNA-mRNA pairs have been reported to participate in muscle development, the actual interactions between them were still unclear and thus in need of further validation. The data from our study provide a rich resource for determining key interactions of miRNA-mRNA in muscle development.

DATA AVAILABILITY
Sequencing data for miRNA have been deposited to Sequence Read Archive at the National Center for Biotechnical Information under accession number PRJNA552780. The RNA-seq data for mRNA were deposited in the Gene Expression Omnibus under accession number GSE73763.

ETHICS STATEMENT
All animal handlings were approved by the Ethical Committee of the Faculty of Veterinary Medicine (EC2013/118) of Ghent University. All methods were performed in accordance with the relevant guidelines and regulations.

AUTHOR CONTRIBUTIONS
ZT and KL designed the experiment and wrote and revised the manuscript. MC and YaY wrote the manuscript and analyzed the data. YiY completed dual luciferase reporter assay. YT and SL were involved in sample collection and total RNA extraction. MZ installed the software used in this study.

FUNDING
This work was supported by the National Natural Science Foundation of China (31830090), the National Key Project (2016ZX08009-003-006), the Shenzhen Science, Technology and Innovation Commission (JCYJ20170307160516413), the Special Fund for Industrial Development of Dapeng New Area at Shenzhen (KY20180114), and the Agricultural Science and Technology Innovation Program (ASTIP-AGIS5). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.