Muya Chen1,2†Yi Long Yao1,2†
Yalan Yang1,2
Min Zhu1,2Yijie Tang1,2Siyuan Liu1,2Kui Li3
Zhonglin Tang1,2,3*- 1Research Centre for Animal Genome, Agricultural Genome Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China
- 2Genome Analysis Laboratory of the Ministry of Agriculture, Agricultural Genome Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China
- 3Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing, China
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 tissue-specific 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.
Materials and Methods
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 reverse-transcription 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 (Tang et al., 2017; 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.
Table 1 Primer sequences.
The HEK293 cells were cultured at 37°C with Dulbecco’s modified Eagle’s medium (Sigma), 10% FBS (Gibco), 1% penicillin/streptomycin (Gibco), and 5% CO2. The miR-and -133b mimics (double-stranded RNA oligonucleotides) and negative control duplexes were synthesized by GenePharma. The pmirGLO-3′UTR-wt, pmirGLO-3′UTR-mt and miRNA (mimic/negative control) were co-transfected into HEK293 cells. The co-transfection assays were performed in 12-well plates with Lipofectamine 2000 reagent (Invitrogen) according to the manufacturer’s instructions and harvested after 24 h. Finally, the dual-luciferase assay system (Promega) was used to examine the activity of renilla and firefly luciferase.
Results
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 (Liang et al., 2017). 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).
Figure 1 The transcriptome profile of multiple tissues in pigs. (A) Genome-wide mean expression value profile. (B) Genome-wide RPKM distribution of mRNA. (C) Hierarchical clustering generated using Pearson correlation coefficients of log2-transformed RPKM (mRNA) values. (D) Composition of the RNA library in each tissue. (E) Genome-wide RPKM distribution of miRNA.
Table 2 Numbers of expressed genes in different tissues.
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).
Table 3 Numbers of miRNA identified in different tissues.
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 (Li et al., 2017), 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 tissue-specific 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).
Table 4 Tissue-specific miRNA in different tissues.
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).
Figure 2 The expression characteristic of tissue-associate genes. (A) Co-expression network of the protein coding genes. (B) Heatmap of tissue-associated mRNA in different pig tissues. (C) Heatmap of tissue-associated miRNA in different pig tissues. (D) GO biological process analysis of tissue-associated mRNA.
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.
Figure 3 Differentially expressed genes during skeletal muscle development. (A) D0 vs. D30 mRNA. (B) D30 vs. D240 mRNA. (C) D0 vs. D30 miRNA. (D) D30 vs. D240 miRNA.
Figure 4 GO analysis of differentially expressed genes during skeletal muscle development. (A) D0 vs. D30 up-regulated genes. (B) D0 vs. D30 down-regulated genes. (C) D30 vs. D240 up-regulated genes. (D) D30 vs. D240 down-regulated genes.
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).
Figure 5 Co-expression pattern analysis for genes during skeletal muscle development.
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. In addition, there were 79 miRNA–mRNA interaction pairs that demonstrated functions in muscle cell differentiation and muscle organ development (Table S16). Of them, 37 miRNA–mRNA pairs showed significant negative expression correlations (r < −0.5) through all three skeletal muscle development stages (Table 5). Many mRNA genes have been reported as regulators of either muscle differentiation or development, including MyoD1 (Blum et al., 2012), CSRP2 (Herrmann et al., 2006), MBNL1 (Chen et al., 2016), FHOD1 (Staus et al., 2011), MET (Park et al., 2015), SOD1 (Sakellariou et al., 2018), RCAN1 (Emrani et al., 2015), SGCA (Fougerousse et al., 1998), SRF (Ding et al., 2017), MEF2A (Yuan et al., 2014), MTM1(Bachmann et al., 2017), LBX1 (Chao et al., 2011), MEF2D (Runfola et al., 2015), IGFBP5 (Zhang et al., 2017), PDGFA (Tallquist et al., 2000), AMOT (Wang et al., 2018a), PDPK1 (Mora et al., 2003), SIRT2 (Arora and Dey, 2014), SIX4 (Chakroun et al., 2015), NOS1 (Villmow et al., 2015), MYL1 (Burguiere et al., 2011), ACTA1 (Hu et al., 2015), PROX1 (Kivela et al., 2016), and QKI (Wu et al., 2017). These interaction pairs included 7 known and 27 novel miRNAs, of which the miRNAs ssc-miR-744 (Yang et al., 2015), ssc-miR-497 (Sato et al., 2014), ssc-miR-338 (Mcdaneld et al., 2009), ssc-miR-423-3p (Siengdee et al., 2015), and ssc-miR-133a-3p (Wang et al., 2018c) were reported to have important roles in myogenesis.
Table 5 The expression correlation of miRNA and their target mRNA during three muscle development stages.
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 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.
Figure 6 Dual luciferase reporter assay. (A) The secondary structure for miRNA and target gene 3′UTR. (B) Schematic of the wild-type and mutation-type 3′UTR vector with miRNA binding sites. Red is the mutation sites, and blue is the wild-type binding sites. (C) luciferase assays were performed.
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 (Yang et al., 2016). 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 tissue-specific genes are well correlated with physiological functions of each organ. For example, the testis-associated gene PAK2 (p21-activated 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 (Tang et al., 2017; 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 (Hou et al., 2016). 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 (Wu et al., 2017). Notable miRNAs that we found include ssc-miR-744, which is reported to significantly up-regulate in muscles after ischemia–reperfusion injury (Yang et al., 2015); 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 (An et al., 2014), 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.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00756/full#supplementary-material
References
An, H. T., Kim, J., Yoo, S., Ko, J. (2014). Small leucine zipper protein (sLZIP) negatively regulates skeletal muscle differentiation via interaction with alpha-actinin-4. J. Biol. Chem. 289, 4969–4979. doi: 10.1074/jbc.M113.515395
Anders, S., Pyl, P. T., Huber, W. (2015). HTSeq—A Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169. doi: 10.1093/bioinformatics/btu638
Arora, A., Dey, C. S. (2014). SIRT2 negatively regulates insulin resistance in C2C12 skeletal muscle cells. Biochim. Biophys. Acta 1842, 1372–1378. doi: 10.1016/j.bbadis.2014.04.027
Bachmann, C., Jungbluth, H., Muntoni, F., Manzur, A. Y., Zorzato, F., Treves, S. (2017). Cellular, biochemical and molecular changes in muscles from patients with X-linked myotubular myopathy due to MTM1 mutations. Hum. Mol. Genet. 26, 320–332. doi: 10.1093/hmg/ddw388
Bai, L., Liang, R., Yang, Y., Hou, X., Wang, Z., Zhu, S., et al. (2015). MicroRNA-21 regulates PI3K/Akt/mTOR signaling by targeting TGFbetaI during skeletal muscle development in pigs. PLoS One 10, e0119396. doi: 10.1371/journal.pone.0119396
Blum, R., Vethantham, V., Bowman, C., Rudnicki, M., Dynlacht, B. D. (2012). Genome-wide identification of enhancers in skeletal muscle: the role of MyoD1. Genes Dev. 26, 2763–2779. doi: 10.1101/gad.200113.112
Burguiere, A. C., Nord, H., Von Hofsten, J. (2011). Alkali-like myosin light chain-1 (myl1) is an early marker for differentiating fast muscle cells in zebrafish. Dev. Dyn. 240, 1856–1863. doi: 10.1002/dvdy.22677
Chakroun, I., Yang, D., Girgis, J., Gunasekharan, A., Phenix, H., Kaern, M., et al. (2015). Genome-wide association between Six4, MyoD, and the histone demethylase Utx during myogenesis. FASEB J. 29, 4738–4755. doi: 10.1096/fj.15-277053
Chao, Z., Wu, J., Zheng, R., Li, F. E., Xiong, Y. Z., Deng, C. Y. (2011). Molecular characterization and expression patterns of Lbx1 in porcine skeletal muscle. Mol. Biol. Rep. 38, 3983–3991. doi: 10.1007/s11033-010-0516-1
Chen, G., Masuda, A., Konishi, H., Ohkawara, B., Ito, M., Kinoshita, M., et al. (2016). Phenylbutazone induces expression of MBNL1 and suppresses formation of MBNL1-CUG RNA foci in a mouse model of myotonic dystrophy. Sci. Rep. 6, 25317. doi: 10.1038/srep25317
Collard, L., Herledan, G., Pincini, A., Guerci, A., Randrianarison-Huetz, V., Sotiropoulos, A. (2014). Nuclear actin and myocardin-related transcription factors control disuse muscle atrophy through regulation of Srf activity. J. Cell. Sci. 127, 5157–5163. doi: 10.1242/jcs.155911
Dabaja, A. A., Mielnik, A., Robinson, B. D., Wosnitzer, M. S., Schlegel, P. N., Paduch, D. A. (2015). Possible germ cell-Sertoli cell interactions are critical for establishing appropriate expression levels for the Sertoli cell-specific microRNA, miR-202-5p, in human testis. Basic Clin. Androl. 25, 2. doi: 10.1186/s12610-015-0018-z
Ding, X., Zhou, S., Li, M., Cao, C., Wu, P., Sun, L., et al. (2017). Upregulation of SRF is associated with hypoxic pulmonary hypertension by promoting viability of smooth muscle cells via increasing expression of Bcl-2. J. Cell. Biochem. 118, 2731–2738. doi: 10.1002/jcb.25922
Dunaeva, M., Blom, J., Thurlings, R., Pruijn, G. J. M. (2018). Circulating serum miR-223-3p and miR-16-5p as possible biomarkers of early rheumatoid arthritis. Clin. Exp. Immunol. 193, 376–385. doi: 10.1111/cei.13156
Emrani, R., Rebillard, A., Lefeuvre, L., Gratas-Delamarche, A., Davies, K. J., Cillard, J. (2015). The calcineurin antagonist RCAN1-4 is induced by exhaustive exercise in rat skeletal muscle. Free Radic. Biol. Med. 87, 290–299. doi: 10.1016/j.freeradbiomed.2015.06.023
Fortes, M. A., Marzuca-Nassr, G. N., Vitzel, K. F., Da Justa Pinheiro, C. H., Newsholme, P., Curi, R. (2016). Housekeeping proteins: how useful are they in skeletal muscle diabetes studies and muscle hypertrophy models? Anal. Biochem. 504, 38–40. doi: 10.1016/j.ab.2016.03.023
Fougerousse, F., Durand, M., Suel, L., Pourquie, O., Delezoide, A. L., Romero, N. B., et al. (1998). Expression of genes (CAPN3, SGCA, SGCB, and TTN) involved in progressive muscular dystrophies during early human development. Genomics 48, 145–156. doi: 10.1006/geno.1997.5160
Friedlander, M. R., Mackowiak, S. D., Li, N., Chen, W., Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Giuffra, E., Kijas, J. M., Amarger, V., Carlborg, O., Jeon, J. T., Andersson, L. (2000). The origin of the domestic pig: independent domestication and subsequent introgression. Genetics 154, 1785–1791.
Groenen, M. A., Archibald, A. L., Uenishi, H., Tuggle, C. K., Takeuchi, Y., Rothschild, M. F., et al. (2012). Analyses of pig genomes provide insight into porcine demography and evolution. Nature 491, 393–398. doi: 10.1038/nature11622
Hashemi, M. S., Mozdarani, H., Ghaedi, K., Nasr-Esfahani, M. H. (2018). Among seven testis-specific molecular markers, SPEM1 appears to have a significant clinical value for prediction of sperm retrieval in azoospermic men. Andrology 6, 890–895. doi: 10.1111/andr.12528
Herrmann, J., Borkham-Kamphorst, E., Haas, U., Van De Leur, E., Fraga, M. F., Esteller, M., et al. (2006). The expression of CSRP2 encoding the LIM domain protein CRP2 is mediated by TGF-beta in smooth muscle and hepatic stellate cells. Biochem. Biophys. Res. Commun. 345, 1526–1535. doi: 10.1016/j.bbrc.2006.05.076
Hou, X., Tang, Z., Liu, H., Wang, N., Ju, H., Li, K. (2012). Discovery of microRNAs associated with myogenesis by deep sequencing of serial developmental skeletal muscles in pigs. PLoS One 7, e52123. doi: 10.1371/journal.pone.0052123
Hou, X., Yang, Y., Zhu, S., Hua, C., Zhou, R., Mu, Y., et al. (2016). Comparison of skeletal muscle miRNA and mRNA profiles among three pig breeds. Mol. Genet. Genomics 291, 559–573. doi: 10.1007/s00438-015-1126-3
Hu, Q., Tong, H., Zhao, D., Cao, Y., Zhang, W., Chang, S., et al. (2015). Generation of an efficient artificial promoter of bovine skeletal muscle alpha-actin gene (ACTA1) through addition of cis-acting element. Cell. Mol. Biol. Lett. 20, 160–176. doi: 10.1515/cmble-2015-0009
Huang Da, W., Sherman, B. T., Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Huang, D. W., Sherman, B. T., Tan, Q., Kir, J., Liu, D., Bryant, D., et al. (2007). DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 35, W169–W175. doi: 10.1093/nar/gkm415
Iwakiri, J., Terai, G., Hamada, M. (2017). Computational prediction of lncRNA–mRNA interactionsby integrating tissue specificity in human transcriptome. Biol. Direct 12, 15. doi: 10.1186/s13062-017-0183-4
Kaminski, M. J., Kaminska, M., Skorupa, I., Kazimierczyk, R., Musial, W. J., Kaminski, K. A. (2013). In-silico identification of cardiovascular disease-related SNPs affecting predicted microRNA target sites. Pol. Arch. Med. Wewn. 123, 355–363. doi: 10.20452/pamw.1819
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14, R36. doi: 10.1186/gb-2013-14-4-r36
Kim, P., Park, A., Han, G., Sun, H., Jia, P., Zhao, Z. (2018). TissGDB: tissue-specific gene database in cancer. Nucleic Acids Res. 46, D1031–D1038. doi: 10.1093/nar/gkx850
Kivela, R., Salmela, I., Nguyen, Y. H., Petrova, T. V., Koistinen, H. A., Wiener, Z., et al. (2016). The transcription factor Prox1 is essential for satellite cell differentiation and muscle fibre-type regulation. Nat. Commun. 7, 13124. doi: 10.1038/ncomms13124
Koh, W., Pan, W., Gawad, C., Fan, H. C., Kerchner, G. A., Wyss-Coray, T., et al. (2014). Noninvasive in vivo monitoring of tissue-specific global gene expression in humans. Proc. Natl. Acad. Sci. U. S. A. 111, 7361–7366. doi: 10.1073/pnas.1405528111
Kozomara, A., Griffiths-Jones, S. (2014). miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi: 10.1093/nar/gkt1181
Kruger, J., Rehmsmeier, M. (2006). RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 34, W451–W454. doi: 10.1093/nar/gkl243
Lage, K., Hansen, N. T., Karlberg, E. O., Eklund, A. C., Roque, F. S., Donahoe, P. K., et al. (2008). A large-scale analysis of tissue-specific pathology and gene expression of human disease genes and complexes. Proc. Natl. Acad. Sci. U. S. A. 105, 20870–20875. doi: 10.1073/pnas.0810772105
Li, B., Qing, T., Zhu, J., Wen, Z., Yu, Y., Fukumura, R., et al. (2017). A comprehensive mouse transcriptomic bodymap across 17 tissues by RNA-seq. Sci. Rep. 7, 4200. doi: 10.1038/s41598-017-04520-z
Li, J. J., Huang, H., Bickel, P. J., Brenner, S. E. (2014). Comparison of D. melanogaster and C. elegans developmental stages, tissues, and cells by modENCODE RNA-seq data. Genome Res. 24, 1086–1101. doi: 10.1101/gr.170100.113
Liang, G., Yang, Y., Niu, G., Tang, Z., Li, K. (2017). Genome-wide profiling of Sus scrofa circular RNAs across nine organs and three developmental stages. DNA Res. 24, 523–535. doi: 10.1093/dnares/dsx022
Lin, H. Y., Cheng, C. H., Chen, D. T., Chen, Y. A., Park, J. Y. (2016). Coexpression and expression quantitative trait loci analyses of the angiogenesis gene–gene interaction network in prostate cancer. Transl. Cancer Res. 5, S951–S963. doi: 10.21037/tcr.2016.10.55
Liu, H., Xi, Y., Liu, G., Zhao, Y., Li, J., Lei, M. (2018). Comparative transcriptomic analysis of skeletal muscle tissue during prenatal stages in Tongcheng and Yorkshire pig using RNA-seq. Funct. Integr. Genomics 18, 195–209. doi: 10.1007/s10142-017-0584-6
Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi: 10.1186/s13059-014-0550-8
Mcdaneld, T. G., Smith, T. P., Doumit, M. E., Miles, J. R., Coutinho, L. L., Sonstegard, T. S., et al. (2009). MicroRNA transcriptome profiles during swine skeletal muscle development. BMC Genomics 10, 77. doi: 10.1186/1471-2164-10-77
Mcloughlin, K. E., Nalpas, N. C., Rue-Albrecht, K., Browne, J. A., Magee, D. A., Killick, K. E., et al. (2014). RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with Mycobacterium bovis. Front. Immunol. 5, 396. doi: 10.3389/fimmu.2014.00396
Mora, A., Davies, A. M., Bertrand, L., Sharif, I., Budas, G. R., Jovanovic, S., et al. (2003). Deficiency of PDK1 in cardiac muscle results in heart failure and increased sensitivity to hypoxia. EMBO J. 22, 4666–4676. doi: 10.1093/emboj/cdg469
Nawrocki, E. P., Burge, S. W., Bateman, A., Daub, J., Eberhardt, R. Y., Eddy, S. R., et al. (2015). Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 43, D130–D137. doi: 10.1093/nar/gku1063
Neville, M. J., Collins, J. M., Gloyn, A. L., Mccarthy, M. I., Karpe, F. (2011). Comprehensive human adipose tissue mRNA and microRNA endogenous control selection for quantitative real-time-PCR normalization. Obesity (Silver Spring) 19, 888–892. doi: 10.1038/oby.2010.257
Nielsen, M., Hansen, J. H., Hedegaard, J., Nielsen, R. O., Panitz, F., Bendixen, C., et al. (2010). MicroRNA identity and abundance in porcine skeletal muscles determined by deep sequencing. Anim. Genet. 41, 159–168. doi: 10.1111/j.1365-2052.2009.01981.x
Park, M., Lee, B. S., Jeon, S. H., Nam, H. J., Lee, G., Kim, C. H., et al. (2015). A novel isoform of met receptor tyrosine kinase blocks hepatocyte growth factor/Met signaling and stimulates skeletal muscle cell differentiation. J. Biol. Chem. 290, 1804–1817. doi: 10.1074/jbc.M114.596957
Peng, Y., Xiang, H., Chen, C., Zheng, R., Chai, J., Peng, J., et al. (2013). MiR-224 impairs adipocyte early differentiation and regulates fatty acid metabolism. Int. J. Biochem. Cell. Biol. 45, 1585–1593. doi: 10.1016/j.biocel.2013.04.029
Ran, M., Chen, B., Wu, M., Liu, X., He, C., Yang, A., et al. (2015). Integrated analysis of miRNA and mRNA expression profiles in development of porcine testes. RSC Adv. 5, 63439–63449. doi: 10.1039/C5RA07488F
Ran, M., Li, Z., Cao, R., Weng, B., Peng, F., He, C., et al. (2018a). miR-26a suppresses autophagy in swine Sertoli cells by targeting ULK2. Reprod. Domest. Anim. 53, 864–871. doi: 10.1111/rda.13177
Ran, M., Weng, B., Cao, R., Li, Z., Peng, F., Luo, H., et al. (2018b). miR-26a inhibits proliferation and promotes apoptosis in porcine immature Sertoli cells by targeting the PAK2 gene. Reprod. Domest. Anim. 53, 1375–1385. doi: 10.1111/rda.13254
Riffo-Campos, A. L., Riquelme, I., Brebi-Mieville, P. (2016). Tools for sequence-based miRNA target prediction: what to choose? Int. J. Mol. Sci. 17, E1987. doi: 10.3390/ijms17121987
Roux, J., Gonzalez-Porta, M., Robinson-Rechavi, M. (2012). Comparative analysis of human and mouse expression data illuminates tissue-specific evolutionary patterns of miRNAs. Nucleic Acids Res. 40, 5890–5900. doi: 10.1093/nar/gks279
Runfola, V., Sebastian, S., Dilworth, F. J., Gabellini, D. (2015). Rbfox proteins regulate tissue-specific alternative splicing of Mef2D required for muscle differentiation. J. Cell. Sci. 128, 631–637. doi: 10.1242/jcs.161059
Sakellariou, G. K., Mcdonagh, B., Porter, H., Giakoumaki, I. I., Earl, K. E., et al. (2018). Comparison of whole body SOD1 knockout with muscle-specific SOD1 knockout mice reveals a role for nerve redox signaling in regulation of degenerative pathways in skeletal muscle. Antioxid. Redox Signal. 28, 275–295. doi: 10.1089/ars.2017.7249
Santos, C. A., Blanck, D. V., De Freitas, P. D. (2014). RNA-seq as a powerful tool for penaeid shrimp genetic progress. Front. Genet. 5, 298. doi: 10.3389/fgene.2014.00298
Sato, T., Yamamoto, T., Sehara-Fujisawa, A. (2014). miR-195/497 induce postnatal quiescence of skeletal muscle stem cells. Nat. Commun. 5, 4597. doi: 10.1038/ncomms5597
Shen, S., Li, D., Liang, J., Wang, J. (2019). Testis-specific calcium-binding protein CBP86-IV (CABYR) binds with phosphoglycerate kinase 2 in vitro and in vivo experiment. Andrologia 10, e13287. doi: 10.1111/and.13287
Siengdee, P., Trakooljul, N., Murani, E., Schwerin, M., Wimmers, K., Ponsuksili, S. (2015). MicroRNAs regulate cellular ATP levels by targeting mitochondrial energy metabolism genes during C2C12 myoblast differentiation. PLoS One 10, e0127850. doi: 10.1371/journal.pone.0127850
Sromek, M., Glogowski, M., Chechlinska, M., Kulinczak, M., Szafron, L., Zakrzewska, K., et al. (2017). Changes in plasma miR-9, miR-16, miR-205 and miR-486 levels after non-small cell lung cancer resection. Cell. Oncol. (Dordr) 40, 529–536. doi: 10.1007/s13402-017-0334-8
Staus, D. P., Blaker, A. L., Medlin, M. D., Taylor, J. M., Mack, C. P. (2011). Formin homology domain-containing protein 1 regulates smooth muscle cell phenotype. Arterioscler. Thromb. Vasc. Biol. 31, 360–367. doi: 10.1161/ATVBAHA.110.212993
Szabo, L., Morey, R., Palpant, N. J., Wang, P. L., Afari, N., Jiang, C., et al. (2015). Statistically based splicing detection reveals neural enrichment and tissue-specific induction of circular RNA during human fetal development. Genome Biol. 16, 126. doi: 10.1186/s13059-015-0690-5
Tallquist, M. D., Weismann, K. E., Hellstrom, M., Soriano, P. (2000). Early myotome specification regulates PDGFA expression and axial skeleton development. Development 127, 5059–5070.
Tang, Z., Wu, Y., Yang, Y., Yang, Y. T., Wang, Z., Yuan, J., et al. (2017). Comprehensive analysis of long non-coding RNAs highlights their spatio-temporal expression patterns and evolutional conservation in Sus scrofa. Sci. Rep. 7, 43166. doi: 10.1038/srep43166
Tang, Z., Yang, Y., Wang, Z., Zhao, S., Mu, Y., Li, K. (2015). Integrated analysis of miRNA and mRNA paired expression profiling of prenatal skeletal muscle development in three genotype pigs. Sci. Rep. 5, 15544. doi: 10.1038/srep15544
Villmow, M., Klockner, U., Heymes, C., Gekle, M., Rueckschloss, U. (2015). NOS1 induces NADPH oxidases and impairs contraction kinetics in aged murine ventricular myocytes. Basic Res. Cardiol. 110, 506. doi: 10.1007/s00395-015-0506-5
Wang, H., Zhang, P., Chen, W., Feng, D., Jia, Y., Xie, L. X. (2012a). Evidence for serum miR-15a and miR-16 levels as biomarkers that distinguish sepsis from systemic inflammatory response syndrome in human subjects. Clin. Chem. Lab. Med. 50, 1423–1428. doi: 10.1515/cclm-2011-0826
Wang, X., Du, C., He, X., Deng, X., He, Y., Zhou, X. (2018a). MiR-4463 inhibits the migration of human aortic smooth muscle cells by AMOT. Biosci. Rep. 38, BSR20180150. doi: 10.1042/BSR20180150
Wang, X., Zhang, P., Li, L., Che, D., Li, T., Li, H., et al. (2018b). miRNA editing landscape reveals miR-34c regulated spermatogenesis through structure and target change in pig and mouse. Biochem. Biophys. Res. Commun. 502, 486–492. doi: 10.1016/j.bbrc.2018.05.197
Wang, X. G., Yu, J. F., Zhang, Y., Gong, D. Q., Gu, Z. L. (2012b). Identification and characterization of microRNA from chicken adipose tissue and skeletal muscle. Poult. Sci. 91, 139–149. doi: 10.3382/ps.2011-01656
Wang, Y., Ma, J., Qiu, W., Zhang, J., Feng, S., Zhou, X., et al. (2018c). Guanidinoacetic acid regulates myogenic differentiation and muscle growth through miR-133a-3p and miR-1a-3p co-mediated Akt/mTOR/S6K signaling pathway. Int. J. Mol. Sci. 19, E2837. doi: 10.3390/ijms19092837
Wu, Y., Li, Z., Yang, M., Dai, B., Hu, F., Yang, F., et al. (2017). MicroRNA-214 regulates smooth muscle cell differentiation from stem cells by targeting RNA-binding protein QKI. Oncotarget 8, 19866–19878. doi: 10.18632/oncotarget.15189
Yang, J. C., Wu, S. C., Rau, C. S., Chen, Y. C., Lu, T. H., Wu, Y. C., et al. (2015). TLR4/NF-kappaB-responsive microRNAs and their potential target genes: a mouse model of skeletal muscle ischemia–reperfusion injury. Biomed. Res. Int. 2015, 410721. doi: 10.1155/2015/410721
Yang, Y., Liang, G., Niu, G., Zhang, Y., Zhou, R., Wang, Y., et al. (2017). Comparative analysis of DNA methylome and transcriptome of skeletal muscle in lean-, obese-, and mini-type pigs. Sci. Rep. 7, 39883. doi: 10.1038/srep39883
Yang, Y., Zhou, R., Mu, Y., Hou, X., Tang, Z., Li, K. (2016). Genome-wide analysis of DNA methylation in obese, lean, and miniature pig breeds. Sci. Rep. 6, 30160. doi: 10.1038/srep30160
Yu, Y., Fuscoe, J. C., Zhao, C., Guo, C., Jia, M., Qing, T., et al. (2014). A rat RNA-seq transcriptomic BodyMap across 11 organs and 4 developmental stages. Nat. Commun. 5, 3230. doi: 10.1038/ncomms4230
Yuan, H., Niu, Y., Liu, X., Fu, L. (2014). Exercise increases the binding of MEF2A to the Cpt1b promoter in mouse skeletal muscle. Acta Physiol. (Oxf) 212, 283–292. doi: 10.1111/apha.12395
Zeng, J., Liu, S., Zhao, Y., Tan, X., Aljohi, H. A., Liu, W., et al. (2016). Identification and analysis of house-keeping and tissue-specific genes based on RNA-seq data sets across 15 mouse tissues. Gene 576, 560–570. doi: 10.1016/j.gene.2015.11.003
Zhang, W. R., Zhang, H. N., Wang, Y. M., Dai, Y., Liu, X. F., Li, X., et al. (2017). miR-143 regulates proliferation and differentiation of bovine skeletal muscle satellite cells by targeting IGFBP5. In Vitro Cell. Dev. Biol. Anim. 53, 265–271. doi: 10.1007/s11626-016-0109-y
Zhang, X., Azhar, G., Williams, E. D., Rogers, S. C., Wei, J. Y. (2015a). MicroRNA clusters in the adult mouse heart: age-associated changes. Biomed. Res. Int. 2015, 732397. doi: 10.1155/2015/732397
Keywords: mRNAs, miRNAs, multiple organ, skeletal muscle, pig
Citation: Chen M, Yao YL, Yang Y, Zhu M, Tang Y, Liu S, Li K and Tang Z (2019) Comprehensive Profiles of mRNAs and miRNAs Reveal Molecular Characteristics of Multiple Organ Physiologies and Development in Pigs. Front. Genet. 10:756. doi: 10.3389/fgene.2019.00756
Received: 12 March 2019; Accepted: 17 July 2019;
Published: 28 August 2019.
Edited by:
Robert J. Schaefer, University of Minnesota Twin Cities, United StatesReviewed by:
Yun Li, Ocean University of China, ChinaYanzhi Jiang, Sichuan Agricultural University, China
Bin Chen, Hunan Agricultural University, China
Copyright © 2019 Chen, Yao, Yang, Zhu, Tang, Liu, Li and Tang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhonglin Tang, tangzhonglin@caas.cn
†These authors have contributed equally to this work.