LncRNAs and CircRNAs Control Myober-Specic Expression Proles by Acting as CeRNAs in Mongolian Horses

Background: The heterogeneity and plasticity of muscle bers are essential for the athletic performance of horses, mainly at the adaption of exercises and the effect on muscle diseases. Skeletal muscle bers can be generally distinguished by their characteristics of contraction as slow and fast type myobers. The diversity of contractile properties and metabolism enable skeletal muscles to respond to the variable functional requirements. We investigated the muscle ber composition and metabolic enzyme activities of splenius muscle and gluteus medius muscle from Mongolian horses. The deep RNA-seq analysis of detecting differentially expressed mRNAs, lncRNAs, circRNAs and their correlation analysis from two muscles were performed. Results: Splenius muscle and gluteus medius muscle from Mongolian horses showed a high divergence of myober compositions and metabolic enzyme activities. Corresponding to their phenotypic characteristics, 94 differentially expressed long noncoding RNAs and 91 differentially expressed circle RNAs were found between two muscles. The analysis results indicate multiple binding sites were detected in lncRNAs and circRNAs with myober-specic expressed miRNAs. Among which we found signicant correlations between the above noncoding RNAs, miRNAs, their target genes, myober-specic developmental transcript factors, and sarcomere genes. Conclusions: We suggest that the ceRNA mechanism of myober-specic expressed noncoding RNAs by acting as miRNA sponges could be ne tuners in regulating skeletal muscle ber composition and transition in horses, which will operate new protective measures of muscle disease and locomotor adaption for racehorses. type ratio was investigated using immunohistochemical staining. The paran blocks with muscle samples were sectioned serially (10 μm thickness) in Microtomes (Leica RM 2245). For myober typing, the reagents and procedures of UltraSensitiveTMS-P Kit (Maixing, Fuzhou, China) were applied for immunohistochemistry. Two antibodies were used to differentiate slow type and fast type muscle in enhancing SOX6 by sponging eca-miR-499, and also affect the muscle developmental transcription factors and sarcomere genes. These noncoding RNAs with their regulatory role will contribute to the study of new mechanisms for muscle plasticity while delivering novel designs of therapeutic conditions for myopathies and locomotor adaptions of skeletal muscle in racehorses.


Background
Mammalian skeletal muscles are heterogeneous, composed of different muscle bers, can be distinguished by their contractile properties as slow-twitch and fast-twitch ber types. Slow muscle bers have a relatively high resistance to fatigue correlated with their thick Z line and rich mitochondrial content, fast type bers have greater contractile speed correlated with their thin Z line and richly developed sarcoplasmic reticulum [1]. Transitions between muscle ber types are correlated with altered physiologic conditions. Endurance training has been shown to lead to a shift from fast-twitch ber types to slow-twitch ber types [2]. Rivero also reported a considerable increase in the proportion of slow-twitch bers has occurred in Arabian horses trained for endurance race [3]. These properties of plasticity and heterogeneity enable horse skeletal muscles to respond to various activities. In addition to affecting the athletic performance and exercise training adaptability of horses, the heterogeneity of skeletal muscle bers has been found associated with some myopathy in horses [4] and mouse models [5][6][7][8][9]. Increasingly, studies suggest that lncRNAs and circRNAs play crucial roles in muscle differentiation and proliferation via acting as competing endogenous RNAs (ceRNAs) or their genomic locations [10][11][12][13][14]. Noncoding RNAs contain binding sites with miRNAs, act as miRNA sponges in the cell, which relieves the inhibitory effect of miRNA on its target genes, increases the expression level of target genes. This competitive mode could be an additional regulator of the horse skeletal muscle regulation network.
The Mongolian horse is the main breed of traditional long-distance races in Inner Mongolia Autonomous Region and Mongolian, has excellent endurance and strong genetic diversity. After learning their superior sports performance and disease resistance, our study investigated the transcriptomic landscape of skeletal muscles with different muscle ber proportions in Mongolian horses using deep RNAseq and extracted the expression pro le of myo ber-speci c mRNAs and ncRNAs from the transcriptome. Association of expression level and binding sites of ncRNAs suggests that long noncoding RNAs and circle RNAs mediate the muscle type-speci c RNA pro les by acting as miRNA sponges.

Materials And Method
Animals and Sample collection In this study, splenius muscle (SPM) samples, gluteus medius muscle (GMM) samples, and additional two muscles (triceps brachii and longissimus dorsi) for correlation analysis were taken from three 5-year-old male Mongolian Horses. All three untrained Mongolian stallions have the same grazing environment and feeding level. The whole protocol of sample collection was approved by the Animal Ethical Committee of Inner Mongolia Agricultural University. All the muscle samples were separated into two pieces, rapidly frozen in liquid nitrogen and stored at -80℃, paraformaldehyde-xed and dehydrated in ascending series of ethanol for -20℃ storage.

Immunohistochemistry
The myo ber type ratio was investigated using immunohistochemical staining. The para n blocks with muscle samples were sectioned serially (10 μm thickness) in Microtomes (Leica RM 2245). For myo ber typing, the reagents and procedures of UltraSensitiveTMS-P Kit (Maixing, Fuzhou, China) were applied for immunohistochemistry. Two antibodies were used to differentiate slow type and fast type muscle bers: Anti-Myosin-7 antibody (rabbit, IHC-1:500) and Anti-Fast Myosin Skeletal Heavy chain antibody (rabbit, IHC-1:500) (all from Bioss, Beijing, China). Four images were captured from each section using the microscopy imaging system (Axio Observer D1, ZEISS). The Image-Pro Plus 6.0 was used for area measurement of muscle bers. The myo ber type ratios were calculated and presented by GraphPad Prism version 6.01 for Windows, GraphPad Software, San Diego, California USA, www.graphpad.com.

Biochemical Analysis
The enzyme activities of phosphofructokinase (PFK) and succinic dehydrogenase (SDH) were estimated to investigate the glycolytic and oxidative capacities of two muscles, respectively. The enzyme assays were performed using the reagents and protocol of PFK Activity Assay Kit (No: BC0530) and SDH Activity Assay Kit (No: BC0950) (all Kits from Beijing Solarbio Science & Technology Co., Ltd, China).
According to the procedure, the degradation rate of NADH which measured at 340 nm is used to re ect the activity of PFK, and the reduction rate of 2,6-dichlorophenol indophenol is determined by the change of absorbance at 600 nm, which represents the activity of the SDH enzyme. The muscle sample (0.1g) frozen by liquid nitrogen was chosen for enzyme extraction. Absorbance is measured using a spectrophotometer (BioPhotometer 6131, Eppendorf). The activities of the two enzymes were statistic by GraphPad Prism version 6.01. Both enzyme activities were represented in U/g fresh weight.
Transcriptome library construction and sequencing Total RNA from SPM, GMM, and additional two muscles were extracted by using the Trizol reagent (Invitrogen, CA, USA), and puri ed with Kit, TRK1002 (LC Science, Houston, TX). RNA quantity and quality were assessed by using NanoDrop 2000, Agilent 2100 Bioanalyzer, and Agilent RNA 6000 Nano Kit (Agilent, CA, USA). Then sequencing libraries were generated using the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). The twelve libraries were sequenced on Illumina Hiseq 4000 platform, and 300 bp (±50 bp) paired-end reads were generated. All the raw data have been deposited in NCBI at BioProject PRJNA573500.
Differentially expressed lncRNAs were identi ed by the ballgown (http://bioconductor.org/packages/release/bioc/html/ballgown.html) R package, with the criteria: |log 2 (Fold change) | >1 and p value < 0.05, and the DESeq2 [24] was used for differential expression analysis of circRNAs. The jvenn interactive Venn diagram viewer was used for computing the intersection [25].

Prediction of target genes
To get insight into molecular functions of lncRNAs and circRNA that differ between two muscles. The target genes of differentially expressed lncRNAs (DELs) and differentially expressed circRNAs (DECs) have been studied. According to the cis-acting regulation of lncRNA, lncRNAs have been demonstrated to modulate the expression of nearby target genes [26]. Protein-coding genes from 100kb upstream and downstream of the DELs were searched. MiRanda [27] and TargetScan [28] were used to predict the lncRNA targets of miRNAs by identifying the binding sites between the miRNA seed and lncRNAs. MiRanda and TargetScan were also used for analyzing the interaction of circRNA-miRNA. MiRNAs are differentially expressed between splenius muscle and gluteus medius muscle in Mongolian horses, which also have been reported in slow and fast type muscles. [29][30][31][32] Interaction between ncRNAs, miRNAs and myo ber-speci c genes The Cytoscape [33] was utilized to construct a network for DELs, DECs and differential expressed miRNAs which had positive expression patterns in SPM and GMM. Correlation analysis of the above ncRNAs, miRNAs, their target genes, myo ber-speci c developmental transcript factors, and sarcomere genes were investigated by R using the Spearman method.
Quantitative RT-PCR validation Eight lncRNAs and 8 circRNAs were selected as candidate genes for RT-PCR validation, all PCR primer sequences are given in (Additional le 1), the convergent primer and discrete primer were specially designed for each circRNA. SYBR ® Premix Ex Taq™ II (TaKaRa) was used to perform the reactions on CFX96 Real-Time PCR Detection System (Bio-Rad). All reactions were performed in triplicate for each sample. Fold change value was calculated by using 2 -ΔΔCt method.

Result
High divergence of phenotypic and biochemical pro les between the SPM and GMM The SPM and GMM differed in composition of MHC (myosin skeletal heavy chain). MHC-I marks slow-twitch myo bers and MHC-II marks fast-twitch myo bers. We investigated the myo ber type of these two muscles ( Fig. 1a-b). We found a signi cantly greater composition of fast-twitch myo bers in gluteus medius muscle (GMM) and a much higher composition of slow-twitch myo bers in splenius muscle (SPM) (Fig. 1c).
The glycolytic enzyme (PFK) and oxidative enzyme (SDH) activities of muscle samples are presented in (Fig. 1d). As shown in gures, two muscles with distinct myo ber compositions had different metabolic enzyme activities, the PFK activities of GMM higher than SPM, the SPM showed a relatively high SDH activity, which corresponding to their myo ber-speci c metabolic properties.
Identi cation of lncRNAs and circRNAs in skeletal muscles with phenotypic differences A total of 598,528,372 raw reads were generated for the SPM and GMM from 3 Mongolian horses by deep mRNA-seq (> 46 million paired reads per sample). On average, 97.19% of the raw reads passed quality control and most reads (89.3-92.9%) were mapped to the reference genome. We obtained 18,794 novel transcripts by comparing the merged gene transfer format le (GTF) and reference.gtf (GCF_002863925.1_EquCab3.0_genomic.gtf). To decrease the false positive rate and strengthen the output quality, identi cation of lncRNAs and circRNAs were all using collections of three forecasting software. The intersection of CPAT, CPC2, and CNCI gives 11731 novel lncRNAs (Additional le 2), and 605 known lncRNAs were identi ed from known transcripts, among which 33 DELs were up-regulated, and 61 were down-regulated in SPM compared to GMM (Fig. 2a). With the split-alignment-based approaches, the intersection of nd_circ, CIRCexplorer2, and CIRI2 gives 853 novel circRNAs (Additional le 2), including 91 DECs, among which 56 highly expressed in SPM, and 35 highly expressed in GMM (Fig. 2b). Validation of lncRNAs and circRNAs were performed for 8 genes respectively using qRT-PCR. The result (Fig. 2c) showed that the same ncRNA expressional trend was found between RNA-seq and qRT-PCR in SPM and GMM. As in the case of circ-EPS15L, the length of PCR products ampli ed using convergent primer and discrete primer were correct, and the divergent primer ampli ed circ-EPS15L in cDNA but not in gDNA (Fig. 2d). These results indicate that the expression pro les from the RNA-seq of lncRNAs and circRNAs were reliable for follow-up studies.

Prediction of target genes
To investigate the function of DELs in cis actions, we detected target genes by searching coding genes 100kb upstream and downstream of the lncRNAs (Additional le 3), among which MYH1, MYH2, MYH4, and MYH8 are MyHC isoforms mainly distributed in fast bers and developing muscle bers [1], which were located near MSTRG.11132.1 and LOC102147674 highly expressed in fast type GMM. MYOZ2 and CLCN1 are involved in muscle contraction, which was located near MSTRG.2829.1 and MSTRG.4532.1 highly expressed in slow type SPM, respectively.
Accumulating experimental evidence indicates that non-coding RNAs served as miRNA sponges or competing endogenous RNAs (ceRNAs) which have an impact on gene expression by binding miRNA, including lncRNAs and circRNAs (Tay et al. 2014). By investigating the binding sites between ncRNAs and miRNAs, ltered by the negative correlation of expression level, 29 interaction relationships (Additional le 4) were found between DELs and miRNAs, and only two interaction relationships were found between DECs and miRNAs.
Construction of DECs, DELs, and differential expressed miRNAs The differential expressed miRNAs between SPM and GMM had multiple binding relationships with DELs or DECs (Fig. 3). The most binding sites of long noncoding RNAs were found in GMM higher expressed lncRNAs. And some DELs with longer sequence lengths had multiple binding sites with different miRNAs, including LOC102150924 (72675 bp) and LOC106781663 (44732 bp). These indicate that a longer lncRNA has further opportunities to sponge the miRNAs, and the ceRNAs mechanism of lncRNA had more advantages in GMM with the higher composition of fast myo bers.

Discussion
In this study, we generated deep mRNA-sequencing data in splenius muscle (SPM) and gluteus medius muscle (GMM) from Mongolian horses and investigated expression pro les of lncRNAs and circRNAs. We found long noncoding RNAs and circle RNAs show a high divergence of expression level between two muscles differed in myo ber type and predominant metabolic enzyme activity. And relativity Analysis of these noncoding RNAs and speci c genes expressed in different muscle ber types were investigated. The results revealed that the distinct myo ber-speci c expression pro les between slow and fast type muscles were predicted to be induced by the differentially expressed noncoding RNAs in horses.
Long noncoding RNA can regulate gene expression in cis or trans actions, and emerging studies suggested that lncRNAs were assigned with functions in myogenic differentiation and muscle regeneration [34][35][36]. A linc-MYH located 54 kb upstream of the MYH2 gene is involved as a regulator in fast muscle ber speci cation through trans-regulation in mice, which is induced by the homeodomain transcription factor Six1. Fast myosin heavy chain genes (also known as MyHC-II) were found expressed in adult skeletal muscles, the gene position and their expression pattern ow the order MYH3 (MyHC-embryonic), MYH2 (MyHC-IIa), MYH1 (MyHC-IIx), MYH4 (MyHC-IIb), MYH8 (MyHC-developing) and MYH13 (MyHC-extraocular) [37,38]. Our results showed that two lncRNAs higher expressed at GMM occupy the fast MYH gene region, a novel lncRNA MSTRG.11132.1 is located between MYH1 and MYH4, and lncRNA LOC102147674 is assigned at the upstream of MYH2 on horse chromosome 11 (Fig. 4a). And we found the expression level of these two lncRNAs positively correlated with MYH1, SIX1, and SIX4, but irrelevant with adjacent MYH2 and MYH3 (Fig. 4b). This suggests that MSTRG.11132.1 and LOC102147674 could interact with six factors in muscles, but could not act in trans to induce the expression of all fast MyHC genes.
In addition to the above cis way, ceRNA mechanisms have been mostly reported in muscle development and proliferation [39][40][41]. Our results revealed that the great majority of binding sites with myo ber-speci c miRNAs were found in lncRNAs with higher expression in GMM. Among these, the LOC102147674 mentioned above was identi ed an exact match to the seed sequence of eca-miR-499 (Fig. 4c). MiR-499 hosted in MYH7B, which involved in post-transcriptional regulation of SOX6, a slow-speci c repressor. The suppression of SOX6 in slow muscle ber type-speci c genes has been reported by many ndings [42,43]. The overexpression of miR-499 leads to a switch from fast to slow ber type of soleus and a transition from MyHC-IIb to MyHC-IIx in II/fast muscles [44], these conversions are on the order of IIb → IIx → IIa → I. We showed a strong negative correlation of LOC102147674 with eca-miR-499 highly expressed in the SPM and positively correlated with the target gene SOX6 (Fig. 4d). A sarcomere is the basic unit of myo bril, the complex components can be divided into different isoforms expressed at slow and fast myo bers. A signi cant correlation was observed between LOC102147674 with a set of 10 fast-speci c sarcomere genes (Fig. 4d). This indicates that LOC102147674 could act as a competitive microRNA sponge of eca-miR-499, thus mediate the target gene SOX6 and ber type-speci c gene programs.
Muscle transcription factors (MTFs) are involved in maintaining and switching the different phenotypes during muscle development. In our study, the transcription factors of the MyoD family, MYOG and MYF5 were found prevalent in the SPM (Additional le 5), in agreement with their role in slow-type muscles. The homeodomain transcription factors Six family are major myogenic transcription factors, our results showed SIX1 and SIX4 expression were also greater in the SPM (Additional le 5). This could suggest that these two Six factors cooperate with MYOG and MYF5 activity in slow-type muscles of horses. The transcriptional coactivators of six factors, protein phosphatases EYA1 and EYA4, were enriched in the GMM (Additional le 5). EYA1 was previously reported can reprogram muscle from slow-twitch muscle into fast-twitch muscle by coactivating with SIX1 [45]. Interestingly, no fast-twitch speci c MyoD and Six transcription factors were found expressed in SPM or GMM of Mongolian horses, this should be investigated in future studies. The FOXO forkhead type transcription factor has a central role in control muscle atrophy. FOXO1 negatively regulates slow-twitch ber genes especially PGC1-α was reported [46]. The FOXO1 and FOXO3 showed higher expression levels in GMM, and FOXO4 higher in SPM (Additional le 5) is consistent with the previous study [47]. Furthermore, corresponding correlations of expression were found between fast-speci c LOC102147674 and transcription factors (Fig. 4d), in particular, negative correlation to slow muscle differentiation factors MYOG and MYF5 with their coactivators SIX1 and SIX4.
CircRNAs also play a central role in regulating skeletal muscle development mostly through sponging miRNA [48][49][50][51][52], but the circRNAs and their regulator role of horse's skeletal muscle remain unexplored. Our results of annotation revealed that some host genes of SPM higher expressed circRNAs involved muscle-related functions, including MYBPC1 in muscle contraction, MYPN in sarcomere organization, GLRX3 in Z disc, and FXN in oxidative phosphorylation and mitochondrion organization are associated with slow oxidative ber metabolism. After the screening of correlation, circ-EIF4G3 and circ-FKBP9 were found to have binding sites with myo ber-speci c expressed miRNAs. The expression pro le of circ-EIF4G3 was identical with the LOC102147674, the result also showed a high score from the binding site with eca-miR-499 (Fig. 4c), a signi cant positive correlation with SOX6, fast-speci c sarcomere genes, and muscle transcription factors (Additional le 5). The circ-FKBP9 showed a higher expression level in SPM, and had a binding site with eca-miR-196b was reported greater in the fast type muscles, but the role in myo ber and myogenesis of this miRNA has not been substantiated by research. Our study suggests that the circRNA pro le also distinguished between slow and fast type muscles, circ-EIF4G3 together with LOC102147674 could in uence myo berspeci c mRNA landscape by acting as a ceRNA for eca-miR-499 in horse skeletal muscles (Fig. 5).
The diversity and plasticity of muscle bers making skeletal muscles can able to cope with altered physiological conditions. Fiber type composition and rearrangements play a key role in the adaptation of various exercise activities for horses [53,54], all that involved muscle development and remodeling. By comparing the transcriptomes of skeletal muscles between SPM and GMM in Mongolian horses, the distinct noncoding RNAs were found base-paired binding sites with myo ber-speci c miRNAs, among which multiple lncRNAs and circRNAs compete with the well-researched slow myo ber gene program transcriptional repressor Sox6. These noncoding RNAs have the potential capacity to modulate the distribution of myo bers and their plasticity.

Conclusion
The distinct transcriptomes of splenius muscles and gluteus medius muscles with considerable variation in myo bers and metabolisms were investigated in Mongolian Horses. The results of deep sequencing showed expressional differences of long noncoding RNAs and circle RNAs between slow and fast type muscles. Multiple binding sites were found in these noncoding RNAs with myo ber-speci c expressed miRNAs. Association of expression levels matched the roles of LOC102147674 with fast-speci c locus and circ-EIF4G3 in enhancing SOX6 by sponging eca-miR-499, and also affect the muscle developmental transcription factors and sarcomere genes. These noncoding RNAs with their regulatory role will contribute to the study of new mechanisms for muscle plasticity while delivering novel designs of therapeutic conditions for myopathies and locomotor adaptions of skeletal muscle in racehorses.

Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval
All animal experiments in this study were approved by the Animal Ethical Committee of Inner Mongolia Agricultural University (Hohhot, China).

Consent for publication
Not applicable.
Competing interests Figure 1 Splenius and gluteus medius muscles from Mongolian horses are differentiated by myo ber composition and oxidative/glycolytic enzyme activities. (a) Section of the splenius muscle. The section was stained with monoclonal antibodies against speci c myosin heavy chain type I (Slow twitch phenotype). Bar=100μm. (b) Section of the gluteus medius muscle. The section was stained with monoclonal antibodies against speci c myosin heavy chain type II (Fast twitch phenotype). Bar=100μm. (c) Comparison of the mean value of myo ber composition (%) between Splenius and gluteus medius muscles. (d) The SDH (oxidative enzyme) and PFK (glycolytic enzyme) activities (U/g fresh weight) of two skeletal muscles Figure 2 Divergence of noncoding RNA expression pro les between splenius and gluteus medius muscles from Mongolian horses. (a) Differentially expressed lncRNAs between two muscles with their muscle types. (b) Differentially expressed circRNAs between two muscles with their muscle types. (c) Comparison of fold changes between RNAseq and qRT-PCR for 8 lncRNAs and 8 circRNAs, both results show the same expression trend. (d) CircRNA validation using PCR with divergent primer. The expected ampli cation size (147 bp) of divergent primer from an example of circ-EPS15L was detected in cDNA but not in genomic DNA (gDNA).