Integrated analysis of miRNA-mediated ceRNA networks in ovine skeletal muscle development

Background: The embryo stage is a key period for sheep skeletal muscle growth and development. Proliferation, differentiation, and hypertrophy of �bers affect muscle growth potential directly. Analyzing transcriptome data is of great signi�cance for revealing important time nodes of fetus muscle development and screening related regulation factors. Muscle development is a complex biological process, including a intricate network of multiple factor interactions. Among them, non-coding RNA, especially miRNA-mediated regulation, plays a �ne regulatory role. The purpose of this study was to investigate the important genes and transcripts involved in the genetic mechanism of embryos skeletal muscle development in late pregnancy. Results: Herein we did a small RNA sequencing(RNA-Seq) of embryo at 85 days (D85N), 105 days (D105N) and 135 days(D135N), then performed bioinformatic analysis in order to identify the miRNA-mediated co-expression networks. Our �ndings identi�ed 505 DE-miRNAs. Integrating the current miRNA data and the previously obtained lncRNA data, multiple networks were constructed, including miRNA-mRNA, miRNA-target gene(TG)-pathway, lncRNA-miRNA-mRNA, and miRNA-TG-transcription factor (TF) network. The results showed that the miRNA-mRNA network and lncRNA-miRNA-mRNA network identi�ed three important lncRNAs (MSTRG.3533, MSTRG.4324, and MSTRG.1470) and three miRNAs(miR-493-3p, miR-3959-3p and miR-410-5p). The four genes ( TEAD1 , ZBTB34 , GSK3B, and POGLUT1 ) and three transcription factors (C / EBPbeta, TFIID, and PR B) play a key regulatory role in the miRNA-TG-TF network. Notably, a similar trend of gene expression was reported by RT-qPCR for RNA-seq data. Conclusions: This study identi�ed three miRNAs, three lncRNAs, four genes, and three transcription factors, and revealed their crucial role in fetal �brogenesis and lipid metabolism. It also shows that D105N is a


Background
Skeletal muscle is one of the most important heterogeneous and complicated tissues in mammals.
During the prenatal and postnatal stages, skeletal muscle shows durative growth, intense plasticity and enormous variability in physiological and functional features [1].The growth rate of prenatal skeletal muscle is signi cantly quicker than that of postpartum, moreover, the muscle development in gestation has a potential effect on postpartum growth performance and muscle maintenance in adult animals.We divided the gestation period into the early and late-stage according to the development characteristics, which included myogenesis and brogenesis [2,3].Previous studies have investigated that myogenesis of early-stage fetus is more intense, in which at least three waves occurred in sheep [4,5], and brogenesis dominances in late-stage accompanied by a slight proliferation of small-diameter ber.Wilson and Mater demonstrated that sheep muscle development is followed through the second half of gestation [6].
Histological characteristics of the ovine fetus showed that myotube fusion, myo ber transformations, and modulation are the main factors to promote late-stage growth [7].As such, the analysis of the temporal and spatial characteristics of ber transformations and a comprehensive molecular understanding of the mechanisms that occur during this stage is of utmost importance to interpret the mechanism of muscle development.
The research of miRNA and the interpretation of its regulatory functions are an important supplement to the content of the genetic central dogma [8].The study of competitive endogenous RNA (ceRNA) network has identi ed new regulatory mechanisms for post-transcriptional regulation, which has facilitated the exploration of the regulation mechanism of the genetic information [9,10].The mechanism of the ceRNA network is that all types of RNA transcripts (lncRNA, circular RNAs, genes and pseudogenes) can communicate with each other by competing with the binding of shared miRNA binding sites (MRE) [11].
In this study, we constructed a putative ceRNA network by integrating lncRNA, miRNA, and mRNA expression based on high-throughput RNA sequencing.Obviously, although the speci c functions of most miRNAs have yet to be elucidated, it is certain that miRNA-mediated regulation is an indispensable and important part of the gene expression regulation mechanism [12,13].More and more evidence have shown that miRNA has a very wide range of gene expression regulation functions, which is closely related to the regulation of gene expression in various life phenomena.
Here we integrate the small RNA sequencing results and the team's previous whole transcriptome sequencing data [14].By constructing an integrated ceRNA network to identify key genes, miRNAs, lncRNAs, transcription factors, and signaling pathways, and study the molecular interaction mechanism which regulates embryonic muscle ber formation, hypertrophy, and maintenance in the late embryo stage.

Small-RNA sequencing
Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure.The total RNA quantity and purity were analyzed of Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0.Approximately 1 ug of total RNA was used to prepare small RNA libraries according to the protocol of TruSeq Small RNA Sample Prep Kits(Illumina, San Diego, USA).And then we performed the single-end sequencing (36 bp or 50 bp ) on an Illumina Hiseq 2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol.The miRNA raw data was uploaded to the GEO (Gene Expression Omnibus) database GSE127287 dataset, it also contains data of lncRNA and mRNA obtained previously.

Differential expression and functional enrichment analysis
Differential expression of miRNAs based on normalized deep-sequencing counts was analyzed by selectively using the Fisher exact test, Chi-squared n × n test, and ANOVA based on the design of the experiment.The signi cance threshold was set to be 0.01 and 0.05 in each test [15,16].In addition, functional enrichment analysis of the differentially expressed miRNAs(DE-miRNAs) identi ed Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, performed using the online resource DAVID v6.8 [17] (https://david.ncifcrf.gov/).Construction of miRNA-mRNA networks miRNA plays a biological role mainly by negatively regulating the expression of downstream genes [18,19].Therefore, we have predicted downstream target genes of differential expressed miRNAs based on the miRanda, TargetScan and miRNet database.miRNet is a comprehensive forecasting database [20,21].The database collects a total of 11 forecasting database results, and the forecasting results are more accurate.Finally, the target genes predicted by the up-regulated (down-regulated) miRNA and the downregulated (up-regulated) mRNA were subjected to intersection analysis to select reliable target genes and construct a miRNA-TG network [22].

Construction of the miRNA-TG-pathway network
The differentially expressed miRNA were uploaded to the BINGO plugin of Cytoscape(v3.6.1).Then, we download the gene annotation le of the corresponding species (http://geneontology.org/) and import it into BINGO with a P threshold for drawing [23].Combined with the results of the miRNA-mRNA network, the same differentially expressed genes(DE-genes) were extracted and uploaded to the ClueGo + CluePedia plugin to enrich the function [24,25].The signi cance threshold was set to 0.05, and the minimum limit of gene clustering was 20.

Construction of integral lncRNA-miRNA-mRNA interaction networks
The Miranda software was used to identify sites on differentially expressed lncRNAs(DE-lncRNAs) to which miRNAs bind to build a lncRNA-miRNA interactome, and the TargetScan explored the potential interactome of DE-miRNAs and its target genes [26,27].Python was used to build relationships between the three types of RNA, and next, a lncRNA-miRNA-mRNA ceRNA network was visualized by Cytoscape, following which, modules representing densely connected nodes were extracted [28,29].

Construction of miRNA-gene-TF network in WGCNA modules
Combined with the WGCNA results of lncRNA and mRNA [14], we selected signi cant modules for integrated analysis in depth.The data set was exported, and functional enrichment analysis was performed on the genes in selected modules using the Cytoscape plugin ClueGO.The DE-genes interacting with the DE-miRNAs were extracted from these signi cant modules, and the miRNA-TG network was merged.miRNA target genes can be involved in cellular activities as transcription factors upstream of the pathway.Therefore, we predict the transcription factors of important genes through the UCSC database(https://genome.ucsc.edu/)and the PROMO database(http://alggen.lsi.upc.es/cgibin/promo_v3/promo/promoinit.cgi?dirDB=TF_8.3).Finally, the gene-transcription factor interaction pairs obtained were tested by the GeneCards database.The tanscription factor genes information refers to TRRUST (http://www.grnpedia.org/trrust/)database and FANTOM 5 (http://fantom.gsc.riken.jp/5/)database.Important genes and transcription factors were extracted and a miRNA-gene-TF network was constructed.

Validation of RNA-Seq data
Five miRNAs were selected to validate the accuracy of RNA sequencing via quantitative real-time PCR(qRT-PCR).Primers were designed for selected transcripts from the transcriptome database(Additional le 1), and real-time PCR was performed using the SYBR Green I master mix (Roche, GmBH, Basel, Switzerland) on the CFX-Connect TM Real-Time System (BioRad).The relative expression of the transcripts was calculated using the ΔΔCt method.

Characterization of miRNA expression pro les
We calculated the length distribution of the total number of valid data, and certain RNA sequences, including rRNA, tRNA, snRNA, and snoRNA, were searched against and ltered from the raw reads by aligning the sequences to RFam and Repbase databases.The length distribution showed that over half of the clean reads (52.6%) were 22 nt in length, which was a good cue that most of the clean reads obtained were miRNA sequences.All clean reads were aligned to the sheep genome and mammalian miRbase using the ACGT-miR101 program.Finally, 4,752 miRNAs were detected in total, and 2,275 of them have never been reported.

Identi cation of DE-miRNAs and functional analysis
The DE-genes were screened with p ≤ 0.05 as a threshold.In the expression pro le, a total of 505 miRNAs were differentially expressed.According to the expression characteristics of pro les, there were ve categories: incremental type (107), decreasing type (151), high-low-high type (91), and low-high-low type (45), and irregularity type(111).We performed a pairwise comparison of three groups on the basis of lter criteria(D85N vs D105N, D105N vs D135N, and D85N vs D135N).There were 63 up-regulated miRNAs and 43 down-regulated miRNAs in D85N vs D105N, 106 up-regulated miRNAs and 123 down-regulated miRNAs from D105N to D135N, and 172 up-regulated miRNAs and 112 down-regulated miRNAs between D85N and D135N.We then performed overlapping statistics on these data, and 16 of the DE-miRNAs were the same in three groups (Fig. 1).The detailed information of these 16 miRNAs was extracted from the miRNA pro le (Table 1).According to the sequence, it is found that 2 of them are the same miRNA.Duplicate data was deleted and 14 miRNAs were used for subsequent network construction.signal transduction processes such as MAPK, Wnt, regulation of actin cytoskeleton, Ras, focal adhesion and the ErbB signaling pathway.Based on the top 20 GO enrichment analysis of miRNAs (Fig. 2a), including biological process, cellular component, molecular function.In biological process, they were riched in protein ubiquitination involved in ubiquitin-dependent protein catabolic process, protein dephosphorylation, protein glycosylation, microtubule cytoskeleton organization, regulation of gastrulation, peptidyl-tyrosine dephosphorylation and actin lament organization, which suggested that spatio-temporal embryonic development was closely related to intracellular protein breakdown and maintenance, energy and material metabolism, signal transduction, and et al.We then further analyzed the 16 miRNAs selected above to better understand their functions.In the top 20 KEGG pathway(Fig.2b), the most highly enriched pathway was nucleotide excision repair.Furthermore, many pathways related to metabolism were also enriched, such as Wnt signaling pathway, insulin signaling pathway, focal adhesion, and ErbB signaling pathway.

miRNA-mRNA networks
With 505 differential expressed miRNAs, more than 100,000 target genes were screened in the database.Therefore, we focused on the 14 important miRNAs selected above for target prediction and screened 1,137 target genes in total.Then, integrating the differential expressed mRNAs in pro le, 379 target genes were retained, and 944 miRNA-mRNA interaction pairs related to skeletal muscle ber development were constructed (Fig. 3).Each miRNA regulates multiple mRNAs and vice versa.A total of 253 target genes were obtained from 7 down-regulated miRNAs and a total of 126 target genes were obtained from 7 upregulated miRNAs.Among them, oar-miR-150 regulates the most target genes, with a total of 100.The average number of target genes regulated by up-regulated and down-regulated miRNAs is more than 35, except for bta-miR-450a.

miRNA-TG-pathway network
On the basis of miRNA-mRNA analysis, we performed an analysis of the miRNA-TG-pathway network in order to nd the biological processes and signaling pathways that the miRNA target genes are mainly involved in.From the results of BINGO, it can be seen that the DE-genes are mainly concentrated in 198 biological processes (Additional le 2), which mainly cover 5 aspects, such as multicellular organismal process, intracellular signal transduction, reproductive process, animal organ development and metabolic process, of which muscle tissue development, actin lament-based process, muscle structure development, and muscle cell differentiation processes have attracted our attention (Fig. 4a).In the diagram drawn by ClueGo and CluePedia (Fig. 4b), 5 pathways were enriched, cytoskeleton organization, positive regulation of synaptic transmission, regulation of small GTPase mediated signal transduction, peptidyl-serine phosphorylation, regulation of neuron differentiation, among which More than 47% of genes are involved in the cytoskeleton tissue process.Finally, the DE-genes involved in these ve pathways were extracted, and miRNAs from these DE-genes were extracted to construct a miRNA-TGpathway regulatory network.In the miRNA-TG-pathway network (Fig. 5), the enriched pathways are Insulin signaling pathway, Wnt signaling pathway, Cell cycle, Focal adhesion, Adherens junction, Ubiquitin mediated proteolysis, ECM-receptor interaction, ErbB signaling pathway, Tight junction, p53 signaling pathway.Among them, Focal adhesion pathway enriched the most differentially expressed target genes.Seven miRNAs enriched in this network were used for in-depth analysis.

miRNA-TG-TF network in WGCNA modules
In this section, we combine the previous WGCNA results [14] to perform in-depth analysis of the screened miRNAs and target genes above.The WGCNA analysis of lncRNA and mRNA was clustered into 25 modules.It can be seen from the that of all the modules, the turquoise module has the highest frequency in the previously study, so it has the most enrichment signi cance.Most differentially expressed RNAs in the ceRNA network fall into this turquoise module.The functional enrichment analysis of DE-genes in the turquoise module shows that these genes participate in many signaling pathways, including ECMreceptor interaction, Adhesion, mTOR signaling pathway, Type diabetes.It is consistent with the above results.The hub genes were extracted from it and predict their TF using USCS.Then, integrating the interaction pairs, a miRNA-TF-gene network was constructed(Fig.7a).The top 30 nodes with the highest degrees were identi ed based on its topological characteristics and were constructed a sub-network (Fig. 7b).It can be seen that 4 genes(TEAD1, ZBTB34, GSK3B and POGLUT1) and 2 miRNAs(oar-miR-493-3p and oar-miR-3959-3p) play important roles in sub-network.GSK3B and POGLUT1 as DE-genes that act as TFs are shown to be targeted by oar-miR-3959-3p and seen to regulate the expression of many TFs (i.e., C / EBPbeta, Sp1, NF1, TFIID and PR B).TEAD1 and ZBTB34 act as TFs are shown to be targeted by oar-miR-493-3p, they can also regulate TFs.In addition, C / EBPbeta, Sp1, NF1, TFIID and PR B can interact with at least three DE-genes.Therefore, according to their functions in skeletal muscle development, C / EBPbeta, TFIID and PR B were selected as important muscle development-related transcription factors.

Validation of RNA-Seq data results
The validation results for the ve miRNAs selected to substantiate the accuracy of sequencing are displayed in Fig. 8.The results indicate that there is a similar expression pattern of miRNAs generated from RNA-Seq and qRT-PCR data.

Discussion
From the data analysis results, in the D85N vs D135N group, there are a great number of DE-miRNAs and DE-genes.The gene expression patterns of D85N and D135N are very different, which supports previous studies [30].There are several genes in muscle tissue that directly or indirectly regulate muscle development, and they are continuously expressed in different stages.We speculate that embryos still have a small range of cell proliferation, differentiation, and migration during the later stages of pregnancy.In D135N, more DE-genes may be related to fat metabolism, muscle disease, carbohydrate metabolism, and protein digestion and absorption, such as DMD, ITGB6, NEDD4L, etc, but not to myo brillar proliferation.This is similar to the results of porcine embryonic skeletal muscle multi-omics studies.Transcription factors are mainly involved in cell proliferation and apoptosis, energy metabolism, adipocyte differentiation, and skeletal muscle contraction [31].Compared with the three groups (D85N vs D135N, D85N vs D105N, and D105N vs D135N), it was found that some DE-genes which exert similar functions were different between D85N and D105N.Therefore, it can be speculated that D105N may be the ber development turning point from myotube differentiation to ber hypertrophy.
Through the construction of multiple networks, genes and transcription factors considered to be key factors were found.TEAD1 is an important member of the TEA domain family [32].It can speci cally regulate the expression of cardiac troponin T, β-myosin heavy chain, smooth muscle α-actin and skeletal muscle α-actin in mammalian and avian skeletal muscle [33,34].In addition, studies have shown that TEAD1 plays a vital role in myoblast growth, skeletal muscle development, muscle ber hypertrophy, muscle regeneration, and myocardial development [35,36].Therefore, during embryonic skeletal muscle development, TEAD1 may be involved in the regulation of muscle ber differentiation, which is similar to the research of some poultry, pig, and mouse [37,38].POGLUT1 is associated with muscle disorders.Muscle satellite cells are activated after muscle damage.Studies have shown that mutations in POGLUT1 inhibit the repair and regeneration of skeletal muscle injury [39,40].ZBTB34 is a new member of the BTB / POZ protein family and gene function is not clear.It can be used as an important transcriptional regulatory protein in the cell growth process to participate in the regulation of the transcriptional activity of certain downstream genes [41].GSK3 is an important factor involved in glucose metabolism, lipid metabolism, insulin and Wnt pathway [42].Activated GSK3 can accelerate the inhibition of adipocyte differentiation and lead to gain body weight [43,44].The transcription factor C / EBPbeta is associated with in ammatory responses, and GSK3 also has a similar function [45].promote the growth and hypertrophy of postpartum skeletal muscle [49].DE-miRNAs and lncRNAs affect insulin by regulating target genes and ultimately participate in the regulation of skeletal muscle cells into adipocytes [50].It shows that these non-coding RNAs play important regulatory roles in the induction of adipogenesis in skeletal muscle and adipose tissue and myoblasts.
In this study, the functional analysis of multiple miRNA-mediated networks in skeletal muscle at different periods, and the embryonic age demarcation of sheep embryonic muscle ber differentiation and hypertrophy was estimated.This study shows that most of the pivotal genes associated with brous hypertrophy and nutritional diseases play a role in fat synthesis, energy metabolism, and immune defense pathways.The intricate energy metabolism pathways indicate that the fetus is slightly affected by the nutritional supply of the parent during the late stages of development(after brogenesis).In addition, some genes and non-coding RNAs are involved in fetal reproductive system development, gonad development, sex differentiation, and other processes.It is speculated that these genes may be related to mammalian sex determination.Regulation of neuron differentiation process may be related to the physiological needs of the embryo during the late stages of development, such as skeletal muscle contraction, and the regulation of fast and slow muscles is innervated.The speci c regulation and control methods need further research.

Conclusion
In conclusion, we constructed a co-expression network with whole transcriptome analysis and revealed that D105N is a pivotal embryo age in the transition from myotube differentiation to ber hypertrophy.Three miRNAs (miR-493-3p, miR-3959-3p, and miR-410-5p), three lncRNAs (MSTRG.3533,MSTRG.4324 and MSTRG.1470),four genes(TEAD1, ZBTB34, GSK3B, and POGLUT1), and three transcription factors (C/EBPbeta, TFIID, and PR B ) were identi ed and emerged as critical in brogenesis and lipid metabolism in the fetus.It provided new ideas to further elucidate the molecular mechanisms of muscle development.

Supplementary Files
This is a list of supplementary les associated with this preprint.Click to download.

Figure 6 a
Figure 6

Figure 7 a
Figure 7

Table 1 The
upTo better understand the functions of DE-miRNAs, we conducted GO term and KEGG pathway analyses.A total of 643 unique biological processes and 73 enriched pathways were identi ed (p ≤ 0.05), including [47]efore, we speculate that in this study, transcription factors such as C / EBPbeta may be closely related to obesity or fat metabolism.PR B and TFIID are basic transcription factors that can recognize both speci c site promoter sequences and chromatin modi cations associated with activated promoters and can function in tissue-speci c and development-speci c pathways[46].miRNA(miR-493-3p,miR-3959-3p, and miR-410-5p) and lncRNA (MSTRG.3533,MSTRG.4324andMSTRG.1470)identied in the ceRNA network were riched in energy metabolism, muscle contraction and oxidation phosphorylation pathway.All pathway analysis results indicate that metabolic and oxidative phosphorylation pathways are signi cantly related to muscle development[47].Heeley et al.showed that high levels of oxidative phosphorylation existed during the rapid development of mammalian skeletal muscle and the formation of myo brils[48].So we predicted that metabolism not only provides energy throughout embryonic skeletal muscle development but may also ne-tune embryonic muscle developmental patterns.Wnt signaling pathway is a potential downstream target of myostatin, which can