Examining the Genetic Background of Porcine Muscle Growth and Development Based on Transcriptome and miRNAome Data

Recently, selection in pigs has been focused on improving the lean meat content in carcasses; this focus has been most evident in breeds constituting a paternal component in breeding. Such sire-breeds are used to improve the meat quantity of cross-breed pig lines. However, even in one breed, a significant variation in the meatiness level can be observed. In the present study, the comprehensive analysis of genes and microRNA expression profiles in porcine muscle tissue was applied to identify the genetic background of meat content. The comparison was performed between whole gene expression and miRNA profiles of muscle tissue collected from two sire-line pig breeds (Pietrain, Hampshire). The RNA-seq approach allowed the identification of 627 and 416 differentially expressed genes (DEGs) between pig groups differing in terms of loin weight between Pietrain and Hampshire breeds, respectively. The comparison of miRNA profiles showed differential expression of 57 microRNAs for Hampshire and 34 miRNAs for Pietrain pigs. Next, 43 genes and 18 miRNAs were selected as differentially expressed in both breeds and potentially related to muscle development. According to Gene Ontology analysis, identified DEGs and microRNAs were involved in the regulation of the cell cycle, fatty acid biosynthesis and regulation of the actin cytoskeleton. The most deregulated pathways dependent on muscle mass were the Hippo signalling pathway connected with the TGF-β signalling pathway and controlling organ size via the regulation of ubiquitin-mediated proteolysis, cell proliferation and apoptosis. The identified target genes were also involved in pathways such as the FoxO signalling pathway, signalling pathways regulating pluripotency of stem cells and the PI3K-Akt signalling pathway. The obtained results indicate molecular mechanisms controlling porcine muscle growth and development. Identified genes (SOX2, SIRT1, KLF4, PAX6 and genes belonging to the transforming growth factor beta superfamily) could be considered candidate genes for determining muscle mass in pigs.


Introduction
In pigs, as one of the major domesticated meat animals, processes regulating muscle growth and development have been widely investigated. However, despite numerous studies focused on the regulation of skeletal muscle growth and performed on various breeds, the obtained findings did not provide a clear view of the process. The next-generation sequencing (NGS) technology provides new possibilities for gene expression measurements and can be helpful to identify the genetic background of phenotypic traits [1].
In pigs, the analysis of global gene expression profiles using NGS technology in semimembranosus and longissimus dorsi muscles indicated genes potentially related to muscle growth [2,3]. Furthermore, the global gene expression profile was also analysed in terms of porcine meat quality [4]. Authors compared muscle transcriptomes of three pig breeds (Pietrain, Polish Landrace, and Pulawska pig) characterized by different production traits and selected genes potentially being under long-term selection focusing on improving the meat content in carcasses. Other studies confirmed the significant role of miRNAs in skeletal muscle growth and development as well as muscle disorders (atrophies and myopathies) [5][6][7]. The comparison of miRNA profiles between prenatal and postnatal ontogenesis periods allowed us to detect miRNAs and their muscle-specific targets associated with muscle development [8]. Tang et al. [9], performing a comprehensive analysis of miRNA-mRNA profiles in muscle tissue of local pig breeds, indicated the set of genes and miRNAs critical to prenatal skeletal muscle development.
A broad range of miRNA effects on gene expression and the translation process indicate the necessity of a comprehensive analysis of global gene and miRNA expression profiles. Such an approach will allow the detection of mechanisms responsible for gene expression regulation in the processes of muscle growth and development. The aim of the present research was to comprehensively identify genes and miRNAs potentially related to muscle growth in pig species. The RNA-seq method was applied to compare whole gene expression and miRNA profiles of muscle tissue collected from two sire-line pig breeds, which are characterized by high muscularity but differ in terms of this feature within each breed.

Transcriptome Analysis-Differentially Expressed Genes (DEGs)
According to the RNA sequencing data, an average of 17.5 M reads per sample were obtained. On average, 98.3% of all sequences passed quality control (QC), of which 75% (12.6 M) mapped to the reference genome. Additionally, 74% of reads were mapped to exons, 18% to intronic regions and 8% to intergenic regions.

Figure 1.
The differentially expressed genes were up-regulated in muscles of pigs with higher muscle mass. Detected genes were coloured and showed no more than 20 interactions (genes marked in grey); line shape indicates the predicted mode of action. Lines coloured indicate the predicted mode of action (red-interactions that were experimentally determined; blue-interactions from curated databases; black-co-expression; green-text mining associations and interactions based on relevant publications mentioning a transfer from other organisms, yellow-transcriptional regulation). Network were presented as a scalable vector graphic and based on Sus scrofa reference. The differentially expressed genes were down-regulated in muscles of pigs with higher muscle mass. Detected genes were coloured and showed no more than 20 interactions (genes marked in grey); line shape indicates the predicted mode of action. Lines coloured indicate the predicted mode of action (blue-interactions from curated databases; black-co-expression). Network were presented as a scalable vector graphic and based on Sus scrofa reference. In both breeds, in the pigs with higher muscularity, the most up-regulated UCP3 gene coding for mitochondrial uncoupling protein 3 was primarily expressed in skeletal muscle and regulating the tissue respiratory processes (fold change 2.73 and 1.68 for Hampshire and Pietrain, respectively). On the other hand, genes showing the greatest significant expression decrease in the pigs with higher loin weight (compared to the animals with lower muscle mass) were MED13 (encoding Mediator Complex Subunit 13; FC = −2.40 and −1.39 for Hp and Pi pigs) and CYP51 (encoding cytochrome P450 family 51 subfamily A member 1; FC = −2.26 and −1.41 for Hp and Pi pigs) genes.
The RNA-seq experiment was submitted to the NCBI Gene Expression Omnibus (GEO) functional genomics data and assigned the GEO accession number GSE107207.

miRNAome Analysis-Differentially Expressed miRNA
As a result of the miRNAome sequencing, between 2.5 M and 7.6 M raw sequences in individual samples were identified. After the filtering, between 2.4 M (94.6%) and 5.8 M (77%) sequences were further mapped to the reference genome. miRNAome sequencing allowed the detection of 196 known and 149 potentially new microRNAs in both breeds. In the Pietrain samples, 160 known and 59 potentially novel miRNAs were identified, while for the Hampshire pigs, 192 known and 138 potentially new miRNAs were detected. The comparison of miRNA profiles between pig groups with different loin weights showed differential expression of 57 microRNAs for Hampshire pigs (24 up-regulated and 34 down-regulated; Figures S1 and S2) and 34 miRNAs for Pietrain pigs (28 up-regulated and 6 down-regulated; Figures S3 and S4). The most significant (p < 0.01) miRNAs detected in Pietrain and Hampshire pigs are shown in Figures 3 and 4. Furthermore, 18 common microRNAs with differential expression were identified for both tested breeds (miR-499a-5p; miR-206; miR-146a-5p; miR-133a-3p; miR-378b; miR-128-3p; miR-378b-3p; miR-10b; miR-451a; miR-125b; miR-30a-5p; let-7f-5p; let-7i-5p; let-7g-5p; miR-7-5p; miR-26a-5p; miR-185-5p; miR-126-3p). This finding may indicate the existence of muscle growth with the regulating pathways common for both investigated high-muscularity pig breeds.  Eighteen DE microRNAs detected in both breeds were further subjected to Gene Ontology and KEGG Pathway analyses (mirPath v.3.0, DIANA-Lab, Athens, Greece). According to an experimentally validated target gene database (Tarbase v7.0, (DIANA-Lab, Athens, Greece)), the most over-represented GO processes were identified, such as cell cycle, fatty acid biosynthesis, regulation of actin cytoskeleton as well as lysine degradation and ubiquitin mediated proteolysis ( Table 2). All identified GO biological processes related to 18 DE microRNAs are shown in Figure 5. The GO terms in each breed and taking into account the direction of expression change is presented in Figures S1-S4.    The miRNA experiment was submitted to the NCBI GEO database and assigned the GEO accession number GSE109215.

Integrated Analysis of Transcriptome and miRNAome
The pathway analysis based on 43 DEGs allowed us to identify three pathways in which two genes coding for collagens (COL6A1 and COL6A2) were involved: the PI3K-Akt signalling (ssc04151), Focal adhesion (ssc04510) and ECM-receptor interaction (ssc04512) pathways. Furthermore, two genes (CYP51-Cytochrome P450 51A; PIGX-Phosphatidylinositol Glycan Anchor Biosynthesis Class X) belonging to the metabolic pathway (ssc01100) were identified. The KEGG pathway analysis based on miRNA results indicated that the most significant pathway was the Hippo signalling pathway, which is also connected with the TGF-β signalling pathway and controls organ size via the regulation of ubiquitin-mediated proteolysis, cell proliferation and apoptosis ( Figure 6). The most interesting target genes involved in identified pathways (FoxO signalling pathway, TGF-β signalling pathway, signalling pathways regulating pluripotency of stem cells and PI3K-Akt signalling pathway) were SOX2, SIRT1, KLF4, PAX6 (Figure 7) and other genes belonging to the transforming growth factor beta superfamily (Table 3). The Hippo signalling pathway identified as the most regulated by detected differentially expressed miRNAs depends on muscle mass. Green frame represents genes engaged in the pathway. Targeted genes regulated by identified miRNAs were marked orange (genes included in more than one pathway) and yellow (genes included only in one pathway) (KEGG database; hsa04390). Black arrows stand for a molecular interaction or relation, while dotted arrows stand for an indirect link or unknown reaction. Figure 7. The signalling pathways regulating pluripotency of stem cells regulated by detected differentially expressed miRNAs depend on muscle mass. Green frame represents genes engaged in the pathway. Targeted genes regulated by identified miRNAs were marked orange (genes included in more than one pathway) and yellow (genes included only in one pathway) (KEGG database; hsa04550). Black arrows stand for a molecular interaction or relation, while dotted arrows stand for an indirect link or unknown reaction. hsa-miR-125b-5p; hsa-miR-30a-5p; hsa-miR-378a-3p; hsa-miR-10b-5p; hsa-let-7i-5p; hsa-let-7g-5p; hsa-let-7f-5p; hsa-miR-7-5p; hsa-miR-26a-5p; hsa-miR-146a-5p; hsa-miR-185-5p; hsa-miR-126-3p; hsa-miR-499a-5p; hsa-miR-451a Moreover, the analysis of miRNA-target gene interactions was performed with the use of the DIANA-mirExTra v2.0 [10] (DIANA-Lab, Athens, Greece) tool to identify functional microRNAs potentially responsible for changes in gene expression identified in the present study as DEGs. mRNA and microRNA differential expression analysis results were simultaneously analysed, and important miRNA regulators were identified based on functional analysis of their target genes. The obtained results confirmed the significant association of 5 detected microRNAs and 9 DEGs (Table 4). Two genes were identified as the most frequently regulated by detected miRNAs, namely, CYP51A1, which is related to cholesterol metabolism (hsa-miR-155-5p; hsa-miR-30c-5p; hsa-miR-199b-5p), as well as the RNF19A gene, which is associated with ubiquitination (hsa-miR-155-5p; hsa-miR-133a-3p; hsa-miR-126-5p) ( Table 4).

Discussion
Since the beginning of domestication, selection in pigs has been focused on improving the lean meat content in carcasses. Currently, this direction of selection is most evident in breeds constituting a paternal component in breeding (Duroc, Pietrain and Hampshire pigs), which are used to improve the meat quantity of cross-breed pig lines. Furthermore, even in one breed, a significant variation in the meatiness level can be observed. In the present study, a comprehensive analysis of genes and microRNA profiles in porcine muscle tissue was performed to identify the genetic background of meat content in carcasses. Moreover, NGS analyses were applied to compare global genes and miRNA expression profiles between groups differing in weight of loin muscle in each breed. Such an approach allowed us to avoid a potential breed-specific effect on transcript abundance and to detect differential expression of genes related to the trait of interest.
The global gene expression profile in a cell or tissue is a result of the transcription rate as well as other post-transcriptional modifications, including mRNA silencing. One of the main gene's expression regulators is microRNA molecules, which predominantly promote transcript degradation [11]. Commonly, miRNAs regulate transcript stability via binding to the 3'UTR of their targeted genes with the RISC multiprotein complex (RNA-induced silencing complex) [12]. microRNAs can control transcript silencing by involvement in transcriptional inhibition through microRNA-mediated chromatin reorganization, mRNA destabilization by decay or cleavage and regulation of proteins belonging to transcription complexes [13,14].
To date, efforts have been made to identify regulatory mechanisms related to muscle growth and development during ontogenesis in pigs. Using a pathway-focused oligo microarray, Li et al. [15] compared the expression profile of genes involved in muscle growth and fatty acid biosynthesis in muscle tissue collected from two pig breeds at six growth stages. Among the genes significantly deregulated during growth periods, the authors detected the UCP3 gene (uncoupling protein 3), whose expression was regulated by 30 other identified genes. The UCP3 gene encodes mitochondrial anion carrier protein, and together with UCP2, it is essential for body metabolism regulation [16]. Both genes are mainly expressed in muscle tissue and control mitochondrial fatty acid transport and glucose metabolism [17]. In our research, we confirmed the differential expression of the UCP3 gene between pigs with varied muscling in both breeds. Furthermore, a miRNA-mRNA target analysis showed an interaction between the UCP3 gene and identified in Hampshire pigs miR-30c-5p, which was established to regulate myoblast differentiation [18]. MicroRNAs belonging to the miR-30 family also affect the activity of miR-206-the main modulator of skeletal muscle development [19,20]. Muscle-specific miR-206 (myomiR) was detected in the present study as differentially expressed in Pietrain and Hampshire pigs differing in muscle weight. This observation confirmed that miR-206 plays an important role in muscle growth and development, regardless of breed tested.
Sheng et al. [21] performed whole miRNAome profiling in pig breeds with high lean meat percentage (Large White) and low body weight (Min Pigs) at several postnatal stages and showed that the one of the most abundant myomiRs was miR-206. Interesting results were reported by Cai et al. [22], who confirmed that the microRNA profile in skeletal muscle was significantly altered by castration. The authors identified seven differentially expressed miRNAs between intact and castrated pigs: miR-206, let-7c and let-7a (down-regulated) and miR-1, miR-133a, miR-26a, and miR-133b (up-regulated). Cai et al. [22] indicated the key impact of castration on miRNA expression, which can explain the main mechanism underlying variations in muscle growth in intact and castrated pigs. Our study allowed us to detect three of seven previously described miRNAs (miR-206; miR-26a; miR-133a). Moreover, we detected the expression of three microRNAs (let7i; let-7g; let-7f) belonging to the let-7 family, which are also known as important developmental regulators.
Previous studies performed by Huang et al. [23] and Yan et al. [24] showed that skeletal muscle hypertrophy can be controlled by the IGF-1-Akt pathway and myostatin signalling pathway. It results from the regulation of IGF-1 or IGFR expression by miR-133a/b and miR-206. Furthermore, during exercise-induced muscle response, IGF-1 signalling can be modified by miR-126, which also impacts MyoD and Myf5 genes strongly related to myogenesis [25]. In humans, several muscle-specific miRNAs, including miR-133a/b and miR-206, are associated with muscular dystrophies (Duchenne and Becker) and are called dystromirs. Their serum expression profile can be used as a biomarker of such muscle disorder [26]. Apart from miR-206 and miR-133a, our results showed significant differences in miR-126-3p expression and significant changes in the expression of its target gene, RNF19A (Ring Finger Protein 19A). It has been established that the RNF19A gene mainly plays a role in neurons causing amyotrophic lateral sclerosis or Parkinson's disease. On the other hand, it was proposed as one of interesting genes related to muscle response during stabilized weight loss [27].
Another gene with modified transcript levels dependent on loin mass was the SETD2 gene, whose expression is regulated by miR-30c-5p. The SETD2 gene, which encodes a histone methylotransferase, is involved in the regulation of transcription processes via histone methylation. Recent studies have shown that the SETD2 gene is involved in epigenetic mechanisms regulating myoblast proliferation and differentiation [28]. Furthermore, using the CRISPR/CAS9 system (Clustered Regularly-Interspaced Short Palindromic Repeats), the authors indicated that silencing the SETD2 gene resulted in a decrease in the expression of major regulators of the cell cycle, repression of myogenin-MyoG-transcription and up-regulation of the cyclin-dependent kinase inhibitor p21. The arrest of the cell cycle caused by the silencing of the SETD2 gene demonstrated that this gene is critical during myoblast proliferation and differentiation [28]. In our study, a significant up-regulation of the SETD2 gene in pigs characterized by higher muscle mass in carcasses was observed. These results support previous findings and indicate a strong relationship between the SETD2 gene and myogenesis in pigs.
The comparison of the global miRNA profile in muscle tissue in the pigs differing in terms of muscularity allowed for the detection of some pathways regulated by identified differentially expressed miRNAs. The most significantly enriched pathways were the Hippo signalling pathway, FoxO signalling pathway and signalling pathways regulating pluripotency of stem cells. Numerous studies have shown that the Hippo signalling pathway plays a key role in muscle cell proliferation, differentiation and apoptosis and, as a result, controls myogenesis processes [24,[29][30][31]. The Hippo signal transduction network can regulate muscle mass and organ size via the activation of cascades of serine/threonine kinases [32,33]. Moreover, Gnimassou et al. [33] suggested that this pathway can contribute to changes and adaptation of muscle mass after exercise. microRNAs identified in the present study can be upstream regulators of the Hippo pathway and, as described previously, are considered to be myomiRs closely related to muscle growth. Our results showed that several miRNAs (miR-30a-5p; miR-206; miR-26a-5p; miR-499a-5p; miR-146a-5p) also regulated the FoxO, TGF-β, PI3K-Akt signalling pathways and signalling pathways regulating pluripotency of stem cells. The identified miRNAs regulate genes coding for proteins belonging to the 14-3-3 family (YWHAE; YWHAB; YWHAQ), which are considered important versatile regulators of the cell cycle, metabolism, cell signalling and apoptosis [34,35]. It was established that 14-3-3 proteins act through regulation of FoxO transcription factors-especially FoXO3, whose significant role in controlling modifications of muscle mass (hypertrophy or atrophy) was confirmed [36]. Furthermore, the identified microRNAs control the transcript abundance of genes regulating different aspects of development: TGF-β proteins (TGFBR1; BMP2; BMP4BMP5; BMPR1B; BMPR2) and genes recognized as pluripotency markers (SOX2; KLF4). The deregulation of the detected DE miRNAs and genes can provide new information about processes related to the growth and development of muscle tissue in pigs.
The PCOLCE gene (Procollagen C-Endopeptidase Enhancer), which is involved in muscle growth disorders in humans, was detected in the set of DEGs in Pietrain and Hampshire breeds. This gene encodes an enzyme responsible for the transformation of procollagens to collagens and is associated with the deregulation of extracellular matrix (ECM) functions and, as a result, with the appearance of late-onset muscle fibrosis and oculopharyngeal muscular [37]. The extracellular matrix plays a key role during signal transduction, and any disorders of this function can lead to myopathies. However, the exact mechanism of the association between impaired ECM matrix function and muscle growth remains unknown. The differential expression of the PCOLCE gene in muscle tissue diverting in mass suggests that this gene can contribute to porcine muscle growth via modification of the ECM role. In both breeds, the differential expression of two genes coding for collagens were detected, namely, COL6A1 and COL6A2. Both COL6A2 and its paralog COL6A1 contribute to the organization of matrix components by binding several extracellular matrix proteins and are related to the occurrence of muscle dystrophy and myopathy (bethlem mopathy 1 and ullrich congenital muscular dystrophy 1) [38,39]. These findings also support the hypothesis that ECM matrix remodelling is associated with the variation of muscle tissue growth and development.
Furthermore, the TBX2 gene coding for the T-box transcription factor was identified in our study as up-regulated in pigs with higher muscle mass. In humans, it was confirmed that TBX2 is responsible for developing tissue sarcomas in which this gene is overexpressed. Zhu et al. [40] indicated that TBX2 regulates the skeletal muscle cell cycle via an interaction with the myogenic transcription factor MyoD and the cell cycle regulators p21 and p14. The authors suggested that in sarcomas, TBX2 promoted cell proliferation, and inhibition of its expression led to the repression of tumour growth [40,41].

Animals
RNA-seq and miRNAome analyses were performed on longissimus lumborum samples collected from 13 pigs representing two sire-line breeds-Pietrain and Hampshire. Animals were maintained at the same housing and feeding conditions according to the procedure previously described [42]. When animals reached a weight of 105 kg (±2.5), they were slaughtered, and after 24 h of chilling, the half carcasses were dissected. Based on dissection data, pigs were selected from a larger population in terms of weight of loin muscle to obtain groups with possibly extreme phenotypes in each breed (Table 7). Pigs included in the present study were unrelated (minimally three crosses back). The research was performed on biological material derived from pigs maintained and slaughtered in the Test Pig Station (National Research Institute of Animal Production), where pigs are slaughtered, dissected and after carcass evaluation, meat is standard intended for consumption. Therefore, research performed does not require the approval of Animal Experimentation committee.
Immediately after slaughter, the tissue samples were collected into the RNAlater solution (Ambion, Thermo Fisher Scientific, Waltham, MA, USA) and stored at −20 • C. Total RNA was isolated from muscle tissue using Direct-zol RNA Mini Prep (Zymo Research, Irvine, CA, USA). The quality and quantity of the obtained RNA were checked using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific; Waltham, MA, USA) and a TapeStation 2200 system (RNAScreen Tape, Agilent Technologies, Santa Clara, CA, USA). Samples with RIN (RNA integrity number) values above 7.5 were used to construct libraries. The pipeline used for raw sequence data quality control (QC), estimation of gene transcript levels and differential expression analysis were described in detail previously [43]. The obtained raw reads were aligned to the Sus scrofa reference genome (assembly Sscrofa10.2 (GCF_000003025.5)). Differentially expressed genes (DEGs; fold change ≥ |1.2|; adjusted p-value < 0.05) were determined separately for each breed using Deseq2 software ( version 1.12.4, [44] between pig groups differing in weight of loin muscle. The GO and pathway analyses were performed using Panther and David software (version 13.1 and v6.8; respectively) [45,46] with Sus scrofa reference. The p-value in enrichment tests was determined based on the Mann-Whitney U Test (Wilcoxon Rank-Sum Test) (Panther software, v 1.12.4 [45], while significantly enriched pathways were detected based on Fisher's exact test (David software, v6.8) [46].

miRNAome Sequencing and Data Analysis
MicroRNA libraries were prepared using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, MA, United States) according to the standard protocol. Briefly, the first step was the 3 adaptor ligation, followed by hybridization with the Reverse Transcription Primer and ligation with the 5 adaptor. The RNA-adaptor ligation products were subjected to reverse transcription. Then, PCR amplification with 12 different indexed primers was performed to allow further multiplexing of the samples. The amplified samples were purified and size-selected on a Novex 6% TBE PAGE gel (Invitrogen, Thermo Fisher Scientific, Waltham, MA USA). After the overnight elution from the gel, the libraries were precipitated and purified with ethanol. Next, they were subjected to a concentration measurement with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and a size assessment with a 2200 TapeStation instrument (Agilent). Finally, the obtained libraries were pooled and sequenced by applying 36 cycles on HiScan SQ (Illumina) according to the manufacturer's protocol.
The obtained reads were converted to FastQ files, subjected to demultiplexing using bcl2fastq software (version, Illumina), and quality controlled using FastQC software [47]. The obtained sequences were analysed using UEA sRNA Workbench v3.2 [48]. The first step was 3 adaptor trimming and filtering according to the following parameters: minimum abundance of at least 6 supporting reads, 17-25 nt of length, low complexity as well as tRNA and rRNA sequences removed. MicroRNAs were identified by applying the miRCat tool v1.0with default parameters for animals, except for minimum abundance (6 reads), minimum length (17 nt), and maximum length (25 nt). The Sus scrofa genome (assembly Sscrofa 10.2, Ensemble database) and miRBase v21.0 (Griffiths-Jones lab at the Faculty of Biolo Manchester; USA) [49,50] were used as references. Detected potentially novel miRNAs were additionally controlled for the presence of other non-coding RNAs using the RNAcentral database (The RNAcentral Consortium) [51]. In the last step, the identification of microRNA length and sequence variants-isomiRs-was performed. The analysis was carried out using isomiR-SEA software (Version 1.60) [52] with the default settings.
miRNAs identified as differentially expressed in the examined samples were analysed using the mirPath v.3.0 DIANA Tools web application (DIANA-Lab, Athens, Greece) [53] to identify their target genes and biological processes enriched by them. DIANA-TarBase v7.0 (DIANA-Lab, Athens, Greece) [54] with experimentally validated target genes was used as a reference database of target genes, and the KEGG Pathway Database was used as a reference database of biological pathways. The analysis was performed on the basis of human homologues deposited in miRBase (21.0) due to the lack of data for pig microRNAs.

Validation of NGS Data by qPCR
The validation of RNA-seq results was performed for 6 genes. The exact transcript abundance was measured using real-time PCR analysis. The primers for all genes were designed using Primer 4.0 software based on the Ensemble reference genome. cDNA was synthesized from 200 ng of total RNA using a High-Capacity RNA-to-cDNA™ Kit (Applied Biosystems, Thermo Fisher, Waltham, MA, USA) according to the protocol. qPCR was performed in 45 cycles on a QuantStudio 7Flex system (Applied Biosystems) with three replicates for each sample (AmpliQ 5x HOT EvaGreen qPCR Mix; Novazym; Poland). OAZ1 and RPL27 were used as endogenous controls. The details of the validated genes are presented in Supplementary File 1- Table. The expression levels were calculated using the delta-delta Ct method [55]. The correlation coefficients between RQ (relative quantity) values and normalized read counts (RNA-seq) were calculated using Pearson's correlation analysis (SAS software, v. 8.02, Wadowice, Poland)).
RT-qPCR method was also employed to validate the expression levels of 9 selected microRNAs in all samples (Table S2). Reverse transcription reactions were performed using a TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific) according to the manufacturer's protocol. Next, quantitative PCR reactions were run using TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific) and commercially available TaqMan microRNA Advanced Assays (Thermo Fisher Scientific), following the protocol and including Non-Template Control (NTC) for each microRNA assay. The reactions were run in triplicate on a QuantStudio 7 Flex Real-Time PCR System (Thermo Fisher Scientific). The efficiency of amplification reactions for each of the miRNA assays was calculated based on the standard curve method.
The NormFinder software [56] was employed to choose an endogenous control that is a microRNA with the most stable expression profile (miR-103a-3p; miR-20a-5p). Finally, relative expression levels of the examined microRNAs were calculated using the ∆∆Ct method, including reaction efficiency E [55].

Conclusions
The comprehensive analysis of the transcriptome and miRNAome profile in porcine muscle tissue differing in terms of muscle mass allowed us to identify differentially expressed genes and miRNAs potentially related to muscle growth and development. We detected microRNAs characteristic of muscle tissue, such as microRNAs belonging to the miR-30 family, miR-206, miR-26a and miR-133a, which have been established as the main modulators of skeletal muscle development. The identification of DE genes and microRNAs enabled us to pinpoint interesting molecular pathways and biological processes essential in controlling muscle growth and development in the pig. The close relationship of some identified in the present study genes and miRNAs with muscle disorders in humans can confirm their significant role in the myogenesis process in pigs. The presented molecular networks show new possibilities for searching candidate genes related to the process of developing muscle tissue and as a result with important production traits in pigs.