The distinct transcriptomes of fast-twitch and slow-twitch muscles in Mongolian horses

Journal Pre-proof The distinct transcriptomes of fast-twitch and slow-twitch muscles in Mongolian horses

Tugeqin Bao, Haige Han, Bei Li, Yiping Zhao, Gerelchimeg Bou, Xinzhuang Zhang, Ming Du, Ruoyang Zhao, Togtokh Mongke, Laxima, Wenqi Ding, Zijie Jia, Manglai Dugarjaviin, Dongyi Bai PII:

S1744-117X(19)30178-9

DOI:

https://doi.org/10.1016/j.cbd.2019.100649

Reference:

CBD 100649

To appear in:

Comparative Biochemistry and Physiology - Part D: Genomics and Proteomics

Received date:

11 October 2019

Revised date:

27 November 2019

Accepted date:

27 November 2019

Please cite this article as: T. Bao, H. Han, B. Li, et al., The distinct transcriptomes of fast-twitch and slow-twitch muscles in Mongolian horses, Comparative Biochemistry and Physiology - Part D: Genomics and Proteomics(2019), https://doi.org/10.1016/ j.cbd.2019.100649

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier.

Journal Pre-proof The distinct transcriptomes of fast-twitch and slow-twitch muscles in Mongolian horses

Tugeqin Baoa, Haige Hana, Bei Lia, Yiping Zhaoa, Gerelchimeg Boua, Xinzhuang Zhanga, Ming Dua, Ruoyang Zhaoa, Togtokh Mongkea, Laximaa, Wenqi Dinga,

a

of

Zijie Jiaa, Manglai Dugarjaviina, #, Dongyi Baia, #, *

Inner Mongolia Key Laboratory of Equine Genetics, Breeding and Reproduction; Scientific

ro

Observing and Experimental Station of Equine Genetics, Breeding and Reproduction,

-p

Ministry of Agriculture and Rural Affairs; Equine Research Center, College of animal science,

re

Inner Mongolia Agricultural University; Hohhot 010018, China

lP

# These co-corresponding authors contributed equally to the study

na

Running title: Transcriptome Analysis of different type muscles in horse

Jo ur

ms. has 24 pages, 4 figures, 1 table, 6 suppl. files

* Corresponding Author:

Dongyi Bai

College of animal science, Inner Mongolia Agricultural University, Hohhot Zhao Wu Da Road 306 010018 Inner Mongolia, China Tel: +86 13848127409; Fax: +86 04716385690 [email protected]

Journal Pre-proof Abstract Skeletal muscle is the largest organ system in the mammalian body and plays a key role in locomotion of horses. Fast and slow muscle fibers have different abilities and functions to adapt to exercises. To investigate the RNA and miRNA expression profiles in the muscles with different muscle fiber compositions on Mongolian horses. We examined the muscle fiber type population and produced deep RNA sequencing for different parts of skeletal muscles. And chose two of them with the highest difference in fast and slow muscle fiber population

of

(splenius and gluteus medius) for comparing the gene expression profile of slow and fast muscle fiber types. We identified a total of 275 differentially expressed genes (DEGs), and 11

ro

differentially expressed miRNAs (DEmiRs). In addition, target gene prediction and

-p

alternative splicing analysis were also performed. Significant correlations were found between the differentially expressed gene, miRNAs, and alternative splicing events. The

re

result indicated that differentially expressed muscle-specific genes and target genes of

na

Mongolian horses.

lP

miRNAs might co-regulating the performance of slow and fast muscle fiber types in

Jo ur

Key words：Mongolian horse; skeletal muscle; fiber type; transcriptome; alternative splicing

1. Introduction

The Mongolian horse is one of the most ancient horse breeds in the world. This breed is characterized by its endurance, strength and disease resistance. Mongolian horses are of a small body, it has strong muscles and great stamina, they can finish the Mongolian traditional horse race with a long distance, up to 30 km (Matthew and Press 2010). Although the Mongolian breed is responsible for the main racing breed in traditional events, its locomotion ability has not been rigorously studied yet. Skeletal muscle constitutes the largest organ system in the mammalian body and is essential for movement and force generation. The skeletal muscle fiber of horses is a characteristic of investigating the breed’s aptitude for athletic activities (Rezende et al. 2016). The physiological, molecular and metabolic properties of skeletal muscle fibers are different,

Journal Pre-proof which can be mainly divided into slow and fast muscle fiber types by their contractile features (Bassel-Duby and Olson 2006). Previous studies suggest that muscle RNA landscapes differed genome-wide between slow and fast muscles (Gao et al. 2017; Raz et al. 2018). Slow muscle fibers contain relatively high capillaries, mitochondria, and myoglobin, they generate energy by oxidative metabolism for continuous contractions with less fatigue. Fast muscle fibers have a high myosin ATPase activity level and glycolytic capacity for fast movements (Schiaffino and Reggiani 2011).

of

As the main adaptations to endurance activities, metabolic patterns and contractile properties of slow and fast muscle fibers are quite different (Russell et al. 2013). The transition between

ro

muscle fibers are correlated with altered physiologic conditions. Endurance training has been

-p

shown to affect metabolic and structural changes in slow and fast muscle fibers, and lead to a switch to slow type (Seene et al. 2004; Handschin et al. 2007).

re

Unraveling the molecular mechanisms that shape the muscle-specific transcriptome will

lP

facilitate an understanding of the regulation of myofiber type switching in physiologic conditions with altered muscle function. In our study, we examined the skeletal muscle fiber

na

type population of Mongolian horses and chose two muscles of them with the highest divergence of slow and fast muscle fiber composition (splenius and gluteus medius) for

Jo ur

RNA-seq. To compare the transcriptome differences among slow and fast type muscles.

2. Materials and methods

2.1

Animals

Four different parts of muscle samples (including splenius, triceps brachii, longissimus dorsi, and gluteus medius) were respectively taken from clinically healthy untrained Mongolian horses (three males; mean age = 5 ± 1 year) by trained veterinary technicians using standard protocols and all those were approved by Animal Ethical Committee of Inner Mongolia Agricultural University. All horses were of a similar level of fitness and the same grazing condition. Every effort was made to minimize horses suffering. To represent every major part of the horse’s body, four muscles were chosen for sampling, including splenius from the neck,

Journal Pre-proof triceps brachii from the forelimb, longissimus dorsi from the back, gluteus medius from the hindlimb. Each muscle was isolated, and then 1 cm3 blocks were taken from the center of the superficial part. These blocks were divided into two parts. Half of the samples were put into liquid nitrogen immediately then stored at -80 ℃. The other half were paraformaldehyde fixed for 20h then gradient ethanol dehydrated to 100% ethanol and stored at -20 ℃.

2.2

Histochemical Analysis

Samples immersed in 100% ethanol were embedded into paraffin blocks. Several cross

of

sections of 10 μm thickness were obtained from each paraffin block on Microtomes (Leica

ro

RM 2245). The UltraSensitiveTMS-P Kit (Maixing, Fuzhou, China) was applied for

-p

immunohistochemistry, operated in accordance with the manufacturer's procedure. The primary antibody was: (1) Anti-Fast Myosin Skeletal Heavy chain antibody (rabbit, 1:500,

re

Bioss, Beijing, China), which specifically reacts with MHC-II; (2) Anti-Myosin-7 antibody

lP

(rabbit, IHC-P=1:500, Bioss, Beijing, China), which specifically reacts with MHC-I. The secondary antibody was: anti-mouse (dilution, 1:200, Maixing, Fuzhou, China). After the

na

staining was performed, signals were detected and measured through digital photographs using the microscopy imaging system (Axio Observer D1, ZEISS). Based on the

Jo ur

immunohistochemical staining images, the fibers were classified as Type I (slow) and II (fast) fibers (Fig.1A), and the population of each muscle fiber type was calculated.

2.3

RNA Library Construction and Sequencing

Total RNA was extracted from twelve samples (3 × 4 sites) using the Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure and Animal Tissue RNA Purification Kit, TRK1002 (LC Science, Houston, TX). The quantity and quality of total RNA were analysis of 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 libraries were sequenced on an Illumina Hiseq 4000 platform, and 300 bp (± 50 bp) paired-end reads were generated. Finally, the raw data was uploaded to the NCBI Sequence Read Archive and the accessions of processed

Journal Pre-proof BioProject

PRJNA573500.

All

raw

sequences

were

deposited

as

BioSamples

SAMN12812445, SAMN12812446, SAMN12812447, SAMN12812448, SAMN12812449, SAMN12812450, SAMN12812451, SAMN12812452, SAMN12812453, SAMN12812454, SAMN12812455 and SAMN12812456.

2.4

RNAseq analysis

The raw reads were first processed using Cutadapt (Martin 2011), and the adaptor sequences

of

and low-quality sequences were filtered from the raw data. The quality of these data was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Then clean

reads

were

mapped

to

the

reference

ro

these

genome

of

Equus

-p

caballus(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/863/925/GCF_002863925.1_EquCa b3.0/GCF_002863925.1_EquCab3.0_genomic.fna.gz). Bowtie2 (Langmead and Salzberg

re

2012) and TopHat (Trapnell et al. 2009) was used for mapping. Mapping results were

lP

evaluated by the RSeQC (Wang et al. 2012) and assembled using Cufflinks (Trapnell et al. 2010). Then, all transcripts generated by each Cufflinks were merged to a more

na

comprehensive transcriptome using Cuffmerge. After the transcriptome was generated, the gene expression levels were estimated by the Cuffdiff.

Jo ur

The expression level of the four group muscles was compared by the differential expression analysis software. FPKM was chosen for differential expression analysis, and identified by the DEseq2 (Love et al. 2014), the filter criteria of differential expression genes (DEGs): |log 2 (Fold Change)| > 1 & qvalue < 0.05. To compare the gene expression profile of fast and slow muscle fiber types, the results of the two groups with the highest difference in fast and slow muscle fiber population (splenius and gluteus medius) were used for the following studies. Functional analysis was conducted using gene ontology (GO) annotations by the clusterProfiler (Yu et al. 2012), and the tools of DAVID (Huang da et al. 2009) was used for KEGG (http://www.genome.jp/kegg/) pathway analysis of differentially expressed genes. Differential splicing, also known as alternative splicing, refers to the process from the precursor of mRNA to the mature mRNA, in which different splicing modes enable the same gene to produce multiple different mature RNA, resulting in different proteins. Alternative

Journal Pre-proof splicing was determined using the rMATS (Shen et al. 2014). This also is a computational tool to detect differential alternative splicing events from RNA-Seq data. The number of reads unique align to transcripts was defined as inclusion level by the statistical model of rMATS. And calculates the P-value with the likelihood-ratio test. Corrected FDR values were used to measure the difference between different groups. The

miRNA

data

were

also

mapped

to

the

reference

genome

of

Equus

caballus(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/863/925/GCF_002863925.1_EquCa

of

b3.0/GCF_002863925.1_EquCab3.0_genomic.fna.gz), miRDeep2 (Friedlander et al. 2012) was used for mapping, and identification of novel and known miRNAs. Differential

target

genes

of

DE

-p

the

ro

expression analysis of miRNA was identified by the edgeR (Robinson et al. 2010), to predict

(http://cbio.mskcc.org/microrna_data/m-iRanda-aug2010.tar.gz)

used

for

target

lP

2.5

was

miRanda

re

prediction.

miRs,

Quantitative RT-PCR validation

na

cDNA was generated from total RNA using the HiScript® II qRT SuperMix for qPCR (Vazyme). Then qPCR was performed with 11 primers (Supplementary File 1) and

Jo ur

SYBR®Premix Ex Taq™ II (TaKaRa) on CFX96 Real-Time PCR Detection System (Bio-Rad). Glyceraldehyde-3-phosphate dehydrogenase (Gapdh) was chosen as the housekeeping gene, which shows stable expression level between splenius and gluteus medius, has been used as a reference gene of skeletal muscles in horse and human. (Aleman and Nieto 2010; Eivers et al. 2012; Jemiolo and Trappe 2004). The statistical analysis method of the Pearson correlation coefficient was applied to investigating the correlation between the results of RNA-seq and qRT-PCR. The R script of correlation analysis: result <- cor (data, method = "pearson").

3. Result

3.1

Muscle fiber type populations

Journal Pre-proof We investigated the muscle fiber populations of skeletal muscles taken from three Mongolian Horses. In total, 12 (4 muscles × 3) samples were investigated including splenius from the neck, triceps brachii from the forelimb, longissimus dorsi from the back, gluteus medius from the hindlimb. As shown in Fig.1B, the type I (slow) muscle fibers populations of Mongolian Horses showed prevalence in splenius (60.83%) followed by triceps brachii (20.50%), longissimus dorsi (14.33%) and gluteus medius (8.87%).

Analysis sequencing and transcriptome annotations

of

3.2

In splenius, 91.17% of clean reads were mapped to the genome. In triceps brachii, 91.93% of

ro

clean reads were mapped to the genome. In longissimus dorsi, 92.30% of clean reads were

-p

mapped to the genome. In gluteus medius, 91.87% of clean reads were mapped to the

Gene expression and identification of DEGs

lP

3.3

re

genome.

In this experiment, four different muscles were compared, discriminated against by the presence of

na

slow and fast myofibers. While data analysis by DESeq2, some genes were differentially expressed in multiple comparison groups, especially the group with the greatest difference in muscle fiber type

Jo ur

distribution (Supplementary File 2). The results showed that there are two groups with the greatest difference in phenotype (Fig.1B) and higher gene expression difference (splenius and gluteus medius). Muscles were discriminated against by the presence of slow and fast myofibers. A total of 386 differentially expressed genes (DEGs) were identified by comparing the gene expression between splenius and gluteus medius (Fig.2 A), from which 260 genes were up-regulated in splenius and 126 genes were up-regulated in gluteus medius.

3.4

GO enrichment of DEGs

Close to three quarters (73%) of DEGs were annotated by the GO enrichment and mainly enriched in biological pathways. In biological processes, the GO terms show a high divergence between the splenius and gluteus medius, the DEGs which higher expressed in

Journal Pre-proof splenius were significantly enriched in skeletal muscle contraction (GO:0003009), cardiac muscle hypertrophy in response to stress (GO:0014898) and muscle cell development (GO:0055001), especially significantly in transition between fast and slow fiber (GO:0014883). Unlike the up-regulated DEGs in splenius, gluteus medius showed up-regulated DEGs in glucan biosynthetic process (GO:0009250) and glycogen biosynthetic process (GO:0005978), which including ATP generation from ADP (GO:0006757) (Fig.2 B).

KEGG pathway analysis

of

3.5

To gain further insight into the molecular functions that differ between splenius and gluteus

ro

medius, the DEGs were evaluated using the functional annotation for the KEGG pathways.

-p

Twelve pathways were significantly enrichment between these two muscles (Table 1). Among these pathways, glycolysis /gluconeogenesis (ecb00010) pathway is a key metabolic

re

pathway in fast-twitch fibers, which differed between two muscles as expected, and this

lP

pathway was greater in gluteus medius with a high proportion of fast muscle fibers. Two circulatory system pathways were higher enrichment in splenius, adrenergic signaling in

na

cardiomyocytes (ecb04261) and cardiac muscle contraction (ecb04260), there are some genes

pathways.

3.6

Jo ur

(Myl2, Myh7, Tnni1, Tnnc1, Tpm3) related to muscle contraction were found among these two

DE microRNAs and miRNA target prediction

MicroRNAs are small non-coding RNAs that can regulate their target genes at the post-transcriptional level, a small RNA library was also constructed for splenius and gluteus medius. After sequence alignment, 691 mature miRNAs and 218 novel miRNAs were found. From those, 11 microRNAs were DE microRNAs among these two muscles (Fig.3 A). There myomiRs (myosin miRNAs), miR-499-3p, miR-499-5p and miR-206 show a higher level in splenius. The miRNAs high in gluteus medius, miR-196a, miR-196b, and miR-10b were previously reported as cancer biomarkers (Ma et al. 2007; Ma 2010; Lu et al. 2016; Mazumder et al. 2019). Through analysis of DE miRNAs target genes, it is predicted that there were 1506 target genes, including previously identified 37 DEGs.

Journal Pre-proof 3.7

Alternative splicing analysis

Alternative splicing was detected using rMATS, five types of alternative splice events were figured out: Exon Skipping (ES), Intron Retention (RI), Mutually Exclusive Exons (MXE), Alternative 5’ splice site (A5SS), Alternative 3’ splice site (A3SS). As shown in Table (Supplementary File 3), ES was the main alternative splice type with the least RI type. In total, 281 gene differential alternative splicing variants were found between splenius and gluteus medius. Among these genes, the inclusion level of some muscle cell development-related

of

genes (Ttn, Mybpc1, Neb, Tnnt1, Tnnt3) were altered between these two muscles in ES. Ttn and Neb exhibited a large number of ES events and show a higher inclusion level in splenius

ro

(Supplementary File 4). Confirming rMATS analyses, the Sox6 was found an exon skipping

-p

in exon 10, which exhibited in both two muscles and show a higher inclusion level in gluteus

The Sox6-mymiRs model

lP

3.8

re

medius (Supplementary File 5).

A co-regulatory model was constructed (Fig.4) in regulating the proportion of fast and slow

na

fiber types, based on our results and findings have been reported. Showing the interactions of Sox6 with Myh7b and myomiRs (miR-206 and miR-499), and exon skipping event happened

3.9

Jo ur

in exon 10 of Sox6.

Quantitative real-time PCR validation

Validation of RNA-seq results was performed for 11 DEGs using quantitative real-time PCR. Those genes were involved in skeletal muscle contraction, transition between fast and slow fiber, and glycolytic metabolism. The result (Supplementary File 6) showed that the similar expression trend of DEGs was found for the RNA-seq and qRT-PCR. The statistical analysis revealed a high correlation between the results of RNA-seq and qRT-PCR. These data indicate that our RNA-seq analysis was reliable for the present study.

4. Discussion

Journal Pre-proof In this study, we investigated the distribution of slow and fast skeletal muscle fibers obtained from four different parts of the Mongolian horse. In addition to the neck, the other muscles contain a much higher percentage of fast type fibers (Fig.1B). This composition is similar to the results of the previous study that investigated the skeletal muscle fiber composition of six Thoroughbred horses (Kawai et al. 2009). The characteristic was that the gluteus medius from hindlimb was significantly higher in fast type fiber than the triceps brachii from forelimb. The predominance of fast type fiber is crucial for the animals to perform powerful energy for a

of

short period of time (Rezende et al. 2016). This finding supports the hypothesis that the hindlimbs muscle perform a more propulsive power than forelimbs during activities (Payne et

ro

al. 2005; Dutto et al. 2006). In general, it is suggested that Type I (slow type) fiber is

-p

indispensable for larger mammals to sustain their heavy body weight (Kawai et al. 2009). The splenius from the neck of the Mongolian horse showed a much higher percentage of slow type

re

fibers, it can be enduring support for its particular large head. Thus, there was a large

lP

variability in the muscle fiber type population of these muscles. To determine and compare the function of slow and fast muscle fiber types, we investigated

na

the transcriptome of splenius and gluteus medius, the function enrichment presents a high divergence between these two muscles. As the metabolism pathways of fast muscle fibers, the

Jo ur

GO results of gluteus medius mainly enriched in pathways related to glycolysis metabolism. The function analysis also showed that DEGs involved in the Calcium signaling pathway and Cardiac muscle contraction were higher expressed at splenius, including some genes Myl2, Myh7, Tnni1, Tnnc1, Tpm3, were mainly expressed in cardiac and slow-twitch skeletal muscle. Several muscle-specific genes and miRNAs were identified among DEGs and DE miRs, miRNAs also play a key role in skeletal muscle differentiation and metabolism (Dumortier et al. 2013; Luo et al. 2013). These could suggest that high divergence between fast and slow type muscles were regulated by the miRNA-mRNA interactive network. Glycogen metabolism rapidly produce energy for quick and forceful muscle contraction during intensive high activity, but depletion of glycogen results in muscle fatigue (Berman and North 2010). The α-actin-3 protein, coded by the Actn3 gene, is critical in muscle contraction, which can regulate the glycogenolysis by influence glycogen phosphorylase in

Journal Pre-proof the Z line (Sjoblom et al. 2008). Because of its potential implication in endurance and sprint performance of humans (Eynon et al. 2013; Pereira et al. 2013), the structure of this gene has been analyzed in many horse breeds. Some single nucleotide polymorphism sites of Actn3 may have an impact on exercise performance. (Mata et al. 2012; Thomas et al. 2014; Ropka-Molik et al. 2019). A transcriptome analysis of skeletal muscle in Arabian horses showed that the expression level of Actn3 decreased with the increasing of training intensity (Ropka-Molik et al. 2017). In Actn3 knockout mice, its deletion exacerbates the shift in

of

muscle metabolism towards more oxidative metabolism (Berman and North 2010). In our study, the Actn3 gene showed higher expression level in the gluteus medius and longissimus

ro

dorsi compared with the splenius, the fast-twitch fiber population of these two muscles was

-p

significantly higher as well. All of these suggested that there may have a connection between myofiber type switching influenced by Actn3 and adaptation to exercises in horses.

re

The calcineurin-NFAT signaling pathway, which is known as a regulator in programming

lP

slow-oxidative skeletal muscle (Naya et al. 2000), dephosphorylation of NFATs by calcineurin promotes their nuclear translocation, and conjunction with myocyte enhancer

na

factor 2 (MEF2), resulting in a slow-twitch muscle gene program (Chin et al. 1998). Two calcineurin related DEGs were found, Sln, and Myoz2, Sln was up-regulated in the fast-twitch

Jo ur

muscle and Myoz2 was up-regulated in the slow-twitch muscle. Sarcolipin (Sln) is a novel regulator of SR Ca2+-ATPase (SERCA) in muscle (Mall et al. 2006; Babu et al. 2007), Sln/SERCA interaction could increase intracellular free Ca2+ and activates calcineurin signaling pathway to maintain the oxidative fiber phenotype (Maurya et al. 2015). Unlike the Actn3, which is also highly expressed in fast-twitch muscle, Sln deficiency can lead to the conversion to slow-to-fast fiber type in the skeleton muscle (Fajardo et al. 2017). The research performed on Myoz2 knockout mice confirmed that the absence of Myoz2 activates the calcineurin signaling and deficiency of Myoz2 increases the composition of slow-twitch muscle fibers (Naya et al. 2000; De Windt et al. 2001). MicroRNAs can regulate expression by promoting degradation or repressing the translation of target transcripts. We found a total of 11 DE miRNAs when comparing splenius and gluteus medius, some muscle-specific enriched miRNAs were identified. From those, miR-499 and

Journal Pre-proof miR-206 were previously reported to involve in muscle performance, metabolism (Hitachi and Tsuchida 2013; Wang 2013) and muscle disorders (McCarthy 2014; Horak et al. 2016). These two muscle-specific miRNAs are also called myomiR, can regulate myoblast differentiation and formation of skeletal muscle fibers through the repression of target genes (Dey et al. 2011; Winbanks et al. 2013). The miR-206 show a higher level in muscles with predominantly slow-twitch fibers than in fast (Allen et al. 2009), which was predicted to regulate the expression of the slow myosin transcriptional repressors (Sox6, Purβ, and Sp3)

of

(Hagiwara et al. 2007; Ji et al. 2007). Unloading of miR-206 caused a shift towards fast-twitch fibers (Allen et al. 2009). In our study, miR-499 and Myh7b show a relatively high

ro

level in splenius. The miR-499 are encoded by introns of the myosin gene Myh7b (Rossi et al.

-p

2010). The specification and maintenance of the slow and fast twitch fiber are regulated by miR-499 through the repressing of the Sox6 mRNA (van Rooij et al. 2009; Duran et al. 2015).

re

MiR-499-5p and miR-206 were higher in splenius compared with soleus and showed a

lP

reverse expression trend to Sox6 (Fig.3 B). Prior studies have shown that horse miR-499 and miR-206 were up-regulated after endurance exercise (Mach et al. 2016), suggesting that slow

na

type muscle-specific myomiRs may play a possible role in continuous exercise by promoting the formation of slow-twitch fibers, which had better resistance to fatigue and high aerobic

Jo ur

energy production capacity.

In different stages of development or tissue, alternative splicing is not invariable, and alternative splicing variants will be generated in specific tissues or conditions. Our results show that few DE alternatively spliced transcripts contribute to muscle cell development and show different inclusion levels between splenius and gluteus medius. The correlation between gene expression level and alternative splicing level of Sox6 was studied. These two values both showed a higher level in splenius. There was a high correlation between reads count numbers of Sox6 and junction reads of exon skipping (Fig.3 C). This suggests that Sox6 expression in slow and fast muscles could be regulated by different levels of alternative splicing. Skeletal muscles are composed of heterogeneous. According to our results and published finds, fast and slow type muscle formation and function require different transcriptional

Journal Pre-proof regulators, microRNAs, and alternative splicing events. MiR-206 and miR-499 are positively associated with splenius which shows a higher composition of slow muscle fibers as they repress transcriptional repressors of slow type muscle formational gene Sox6. We found alternative splicing in Sox6, and show a high divergence between slow and fast fiber type muscles. The different levels of alternative splicing events may also affect the expression level of the Sox6 between slow and fast type muscles. Our results indicate that mRNA

of

processing mechanisms contribute to muscle-specific expression profiles in skeletal muscles.

5. Conclusion

ro

The DEG-myomiR co-regulatory network profiled in skeletal muscle reflects the specification

-p

and maintenance of the fast and slow fiber type muscles of the Mongolian horse. Actn3, Myoz2, and Sln showed a different expression level between gluteus medius and splenius,

re

which can influence the transition between fast and slow fiber by regulating glycogen

lP

metabolism and calcineurin-NFAT signaling pathway. The Sox6-myomiR model plays a

na

possible role in regulating the performance of fast and slow muscle fiber types.

Acknowledgments. Our work was supported by the National Natural Science Foundation of

Jo ur

China [grant numbers 31960657, 31860636, 31902188, 31360538,31472070]; the National Key Research and Development Program of China [grant number 2017YFE0108700]; the Natural Science Foundation Major Project of Inner Mongolia [grant numbers 2017ZD06, 2019MS03064]; the project of special funds for landmark achievements of College of Animal Science, Inner Mongolia Agricultural University [grant number BZCG201901]; the project of talent cultivation project of double-first-class discipline innovation team construction of Inner Mongolia Agricultural University [grant number NDSC2018-05]; the project of special funds for subsidies after scientific and technological innovation of herbivorous livestock of College of Animal Science, Inner Mongolia Agricultural University [grant number QN201909]; the project of special funds for transforming scientific and technological achievements breeding (cultivating) for new varieties of animals and plants of Inner Mongolia Agricultural

Journal Pre-proof University [grant number YZGC2017001]; the startup project for High-level Talents of Inner Mongolia Agricultural University [grant number NDGCC2016-01].

Author contributions. B.T. designed the experiments, performed histological experiments, analyzed the data and wrote the manuscript, L.B. applied the Transcriptional Profiling, B.D.Y. assisted with experimental design, Z.Y.P, T.M, Z.R.Y. and the others helped with sample collection. G.B offered help in lab training. B.D.Y. and M.D. were the supervisor of this

ro

of

project and secured funding, conceived, revised and submitted the manuscript.

-p

Conflicts of interest. The authors declared that they have no conflicts of interest to this work.

lP

re

References

na

Aleman, M., and J. E. Nieto, 2010. Gene expression of proteolytic systems and growth regulators of skeletal muscle in horses with myopathy associated with pituitary pars

Jo ur

intermedia dysfunction. Am J Vet Res. 71, 664-70. Allen, D. L., E. R. Bandstra, B. C. Harrison, S. Thorng, L. S. Stodieck, P. J. Kostenuik, S. Morony, D. L. Lacey, T. G. Hammond, L. L. Leinwand, W. S. Argraves, T. A. Bateman, and J. L. Barth, 2009. Effects of spaceflight on murine skeletal muscle gene expression. J Appl Physiol (1985). 106, 582-95. Babu, G. J., P. Bhupathy, V. Timofeyev, N. N. Petrashevskaya, P. J. Reiser, N. Chiamvimonvat, and M. Periasamy, 2007. Ablation of sarcolipin enhances sarcoplasmic reticulum calcium transport and atrial contractility. Proc Natl Acad Sci U S A. 104, 17867-72.

Journal Pre-proof Bassel-Duby, R., and E. N. Olson, 2006. Signaling pathways in skeletal muscle remodeling. Annu Rev Biochem. 75, 19-37. Berman, Y., and K. N. North, 2010. A gene for speed: the emerging role of alpha-actinin-3 in muscle metabolism. Physiology (Bethesda). 25, 250-9. Chin, E. R., E. N. Olson, J. A. Richardson, Q. Yang, C. Humphries, J. M. Shelton, H. Wu, W.

of

Zhu, R. Bassel-Duby, and R. S. Williams, 1998. A calcineurin-dependent

ro

transcriptional pathway controls skeletal muscle fiber type. Genes Dev. 12, 2499-509. De Windt, L. J., H. W. Lim, O. F. Bueno, Q. Liang, U. Delling, J. C. Braz, B. J. Glascock, T. F.

-p

Kimball, F. del Monte, R. J. Hajjar, and J. D. Molkentin, 2001. Targeted inhibition of

lP

3322-7.

re

calcineurin attenuates cardiac hypertrophy in vivo. Proc Natl Acad Sci U S A. 98,

na

Dey, B. K., J. Gagan, and A. Dutta, 2011. miR-206 and -486 induce myoblast differentiation by downregulating Pax7. Mol Cell Biol. 31, 203-14.

Jo ur

Dumortier, O., C. Hinault, and E. Van Obberghen, 2013. MicroRNAs and metabolism crosstalk in energy homeostasis. Cell Metab. 18, 312-24. Duran, B. O., G. J. Fernandez, E. A. Mareco, L. N. Moraes, R. A. Salomao, T. Gutierrez de Paula, V. B. Santos, R. F. Carvalho, and M. Dal-Pai-Silva, 2015. Differential microRNA Expression in Fast- and Slow-Twitch Skeletal Muscle of Piaractus mesopotamicus during Growth. PLoS One. 10, e0141967. Dutto, D. J., D. F. Hoyt, H. M. Clayton, E. A. Cogger, and S. J. Wickler, 2006. Joint work and power for both the forelimb and hindlimb during trotting in the horse. J Exp Biol. 209, 3990-9.

Journal Pre-proof Eivers, S. S., B. A. McGivney, J. Gu, D. E. MacHugh, L. M. Katz, and E. W. Hill, 2012. PGC-1alpha encoded by the PPARGC1A gene regulates oxidative energy metabolism in equine skeletal muscle during exercise. Anim Genet. 43, 153-62. Eynon, N., E. D. Hanson, A. Lucia, P. J. Houweling, F. Garton, K. N. North, and D. J. Bishop, 2013. Genes for elite power and sprint performance: ACTN3 leads the way. Sports

of

Med. 43, 803-17.

ro

Fajardo, V. A., B. A. Rietze, P. J. Chambers, C. Bellissimo, E. Bombardier, J. Quadrilatero, and A. R. Tupling, 2017. Effects of sarcolipin deletion on skeletal muscle adaptive

-p

responses to functional overload and unload. Am J Physiol Cell Physiol. 313,

re

C154-c61.

lP

Friedlander, M. R., S. D. Mackowiak, N. Li, W. Chen, and N. Rajewsky, 2012. miRDeep2

na

accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37-52.

Jo ur

Gao, K., Z. Wang, X. Zhou, H. Wang, D. Kong, C. Jiang, X. Wang, Z. Jiang, and X. Qiu, 2017. Comparative transcriptome analysis of fast twitch muscle and slow twitch muscle in Takifugu rubripes. Comp Biochem Physiol Part D Genomics Proteomics. 24, 79-88. Hagiwara, N., M. Yeh, and A. Liu, 2007. Sox6 is required for normal fiber type differentiation of fetal skeletal muscle in mice. Dev Dyn. 236, 2062-76. Handschin, C., S. Chin, P. Li, F. Liu, E. Maratos-Flier, N. K. Lebrasseur, Z. Yan, and B. M. Spiegelman, 2007. Skeletal muscle fiber-type switching, exercise intolerance, and myopathy in PGC-1alpha muscle-specific knock-out animals. J Biol Chem. 282, 30014-21.

Journal Pre-proof Hitachi, K., and K. Tsuchida, 2013. Role of microRNAs in skeletal muscle hypertrophy. Front Physiol. 4, 408. Horak, M., J. Novak, and J. Bienertova-Vasku, 2016. Muscle-specific microRNAs in skeletal muscle development. Dev Biol. 410, 1-13. Huang da, W., B. T. Sherman, and R. A. Lempicki, 2009. Bioinformatics enrichment tools:

of

paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids

ro

Res. 37, 1-13.

Jemiolo, B., and S. Trappe, 2004. Single muscle fiber gene expression in human skeletal

-p

muscle: validation of internal control with exercise. Biochem Biophys Res Commun.

re

320, 1043-50.

lP

Ji, J., G. L. Tsika, H. Rindt, K. L. Schreiber, J. J. McCarthy, R. J. Kelm, Jr., and R. Tsika, 2007.

na

Puralpha and Purbeta collaborate with Sp3 to negatively regulate beta-myosin heavy chain gene expression during skeletal muscle inactivity. Mol Cell Biol. 27, 1531-43.

Jo ur

Kawai, M., Y. Minami, Y. Sayama, A. Kuwano, A. Hiraga, and H. Miyata, 2009. Muscle fiber population and biochemical properties of whole body muscles in Thoroughbred horses. Anat Rec (Hoboken). 292, 1663-9. Langmead, B., and S. L. Salzberg, 2012. Fast gapped-read alignment with Bowtie 2. Nat Methods. 9, 357-9. Love, M. I., W. Huber, and S. Anders, 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550.

Journal Pre-proof Lu, Y. C., J. T. Chang, E. C. Chan, Y. K. Chao, T. S. Yeh, J. S. Chen, and A. J. Cheng, 2016. miR-196, an Emerging Cancer Biomarker for Digestive Tract Cancers. J Cancer. 7, 650-5. Luo, W., Q. Nie, and X. Zhang, 2013. MicroRNAs involved in skeletal muscle differentiation. J Genet Genomics. 40, 107-16.

of

Ma, L., 2010. Role of miR-10b in breast cancer metastasis. Breast Cancer Res. 12, 210.

ro

Ma, L., J. Teruya-Feldstein, and R. A. Weinberg, 2007. Tumour invasion and metastasis initiated by microRNA-10b in breast cancer. Nature. 449, 682-8.

-p

Mach, N., S. Plancade, A. Pacholewska, J. Lecardonnel, J. Riviere, M. Moroldo, A. Vaiman, C.

re

Morgenthaler, M. Beinat, A. Nevot, C. Robert, and E. Barrey, 2016. Integrated mRNA

lP

and miRNA expression profiling in blood reveals candidate biomarkers associated

na

with endurance exercise in the horse. Sci Rep. 6, 22932. Mall, S., R. Broadbridge, S. L. Harrison, M. G. Gore, A. G. Lee, and J. M. East, 2006. The

Jo ur

presence of sarcolipin results in increased heat production by Ca(2+)-ATPase. J Biol Chem. 281, 36597-602.

Martin, Marcel, 2011. Cutadapt removes adapter sequences from high-throughput sequencing reads. Embnet Journal. 17. Mata, X., A. Vaiman, A. Ducasse, M. Diribarne, L. Schibler, and G. Guerin, 2012. Genomic structure, polymorphism and expression of the horse alpha-actinin-3 gene. Gene. 491, 20-4. Matthew, Davis, and St. Martin'S Press, 2010. When Things Get Dark: A Mongolian Winter`s Tale.

Journal Pre-proof Maurya, S. K., N. C. Bal, D. H. Sopariwala, M. Pant, L. A. Rowland, S. A. Shaikh, and M. Periasamy, 2015. Sarcolipin Is a Key Determinant of the Basal Metabolic Rate, and Its Overexpression Enhances Energy Expenditure and Resistance against Diet-induced Obesity. J Biol Chem. 290, 10840-9. Mazumder, S., S. Datta, J. G. Ray, K. Chaudhuri, and R. Chatterjee, 2019. Liquid biopsy:

of

miRNA as a potential biomarker in oral cancer. Cancer Epidemiol. 58, 137-45.

ro

McCarthy, J. J., 2014. microRNA and skeletal muscle function: novel potential roles in exercise, diseases, and aging. Front Physiol. 5, 290.

-p

Naya, F. J., B. Mercer, J. Shelton, J. A. Richardson, R. S. Williams, and E. N. Olson, 2000.

lP

Chem. 275, 4545-8.

re

Stimulation of slow skeletal muscle fiber gene expression by calcineurin in vivo. J Biol

na

Payne, R. C., P. Veenman, and A. M. Wilson, 2005. The role of the extrinsic thoracic limb muscles in equine locomotion. J Anat. 206, 193-204.

Jo ur

Pereira, A., A. M. Costa, M. Izquierdo, A. J. Silva, E. Bastos, and M. C. Marques, 2013. ACE I/D and ACTN3 R/X polymorphisms as potential factors in modulating exercise-related phenotypes in older women in response to a muscle power training stimuli. Age (Dordr). 35, 1949-59. Raz, V., M. Riaz, Z. Tatum, S. M. Kielbasa, and P. A. C. t Hoen, 2018. The distinct transcriptomes of slow and fast adult muscles are delineated by noncoding RNAs. Faseb j. 32, 1579-90. Rezende, Adalgiza Souza Carneiro De, Mayara Gonçalves Fonseca, Lilian De Rezende Jordão, Maria Luiza Mendes De Almeida, Antônio De Queiroz Neto, Guilherme De

Journal Pre-proof Camargo Ferraz, and José Luis L. Rivero, 2016. Skeletal Muscle Fiber Composition of Untrained Mangalarga Marchador Fillies. Journal of Equine Veterinary Science. 36, 101-04. Robinson, M. D., D. J. McCarthy, and G. K. Smyth, 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26,

of

139-40.

ro

Ropka-Molik, K., M. Stefaniuk-Szmukier, Z. Ukowski K, K. Piorkowska, and M. Bugno-Poniewierska, 2017. Exercise-induced modification of the skeletal muscle

-p

transcriptome in Arabian horses. Physiol Genomics. 49, 318-26.

re

Ropka-Molik, K., M. Stefaniuk-Szmukier, A. D. Musial, K. Piorkowska, and T. Szmatola, 2019.

lP

Sequence analysis and expression profiling of the equine ACTN3 gene during

na

exercise in Arabian horses. Gene. 685, 149-55. Rossi, A. C., C. Mammucari, C. Argentini, C. Reggiani, and S. Schiaffino, 2010. Two

Jo ur

novel/ancient myosins in mammalian skeletal muscles: MYH14/7b and MYH15 are expressed in extraocular muscles and muscle spindles. J Physiol. 588, 353-64. Russell, A. P., S. Lamon, H. Boon, S. Wada, I. Guller, E. L. Brown, A. V. Chibalin, J. R. Zierath, R. J. Snow, N. Stepto, G. D. Wadley, and T. Akimoto, 2013. Regulation of miRNAs in human skeletal muscle following acute endurance exercise and short-term endurance training. J Physiol. 591, 4637-53. Schiaffino, S., and C. Reggiani, 2011. Fiber types in mammalian skeletal muscles. Physiol Rev. 91, 1447-531.

Journal Pre-proof Seene, T., P. Kaasik, K. Alev, A. Pehme, and E. M. Riso, 2004. Composition and turnover of contractile proteins in volume-overtrained skeletal muscle. Int J Sports Med. 25, 438-45. Shen, S., J. W. Park, Z. X. Lu, L. Lin, M. D. Henry, Y. N. Wu, Q. Zhou, and Y. Xing, 2014. rMATS: robust and flexible detection of differential alternative splicing from replicate

of

RNA-Seq data. Proc Natl Acad Sci U S A. 111, E5593-601.

ro

Sjoblom, B., A. Salmazo, and K. Djinovic-Carugo, 2008. Alpha-actinin structure and regulation. Cell Mol Life Sci. 65, 2688-701.

-p

Thomas, K. C., N. A. Hamilton, K. N. North, and P. J. Houweling, 2014. Sequence analysis of

re

the equine ACTN3 gene in Australian horse breeds. Gene. 538, 88-93.

lP

Trapnell, C., L. Pachter, and S. L. Salzberg, 2009. TopHat: discovering splice junctions with

na

RNA-Seq. Bioinformatics. 25, 1105-11.

Trapnell, C., B. A. Williams, G. Pertea, A. Mortazavi, G. Kwan, M. J. van Baren, S. L. Salzberg,

Jo ur

B. J. Wold, and L. Pachter, 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 28, 511-5.

van Rooij, E., D. Quiat, B. A. Johnson, L. B. Sutherland, X. Qi, J. A. Richardson, R. J. Kelm, Jr., and E. N. Olson, 2009. A family of microRNAs encoded by myosin genes governs myosin expression and muscle performance. Dev Cell. 17, 662-73. Wang, L., S. Wang, and W. Li, 2012. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 28, 2184-5.

Journal Pre-proof Wang, X. H., 2013. MicroRNA in myogenesis and muscle atrophy. Curr Opin Clin Nutr Metab Care. 16, 258-66. Winbanks, C. E., C. Beyer, A. Hagg, H. Qian, P. V. Sepulveda, and P. Gregorevic, 2013. miR-206 represses hypertrophy of myogenic cells but not muscle fibers via inhibition of HDAC4. PLoS One. 8, e73589.

of

Yu, G., L. G. Wang, Y. Han, and Q. Y. He, 2012. clusterProfiler: an R package for comparing

Jo ur

na

lP

re

-p

ro

biological themes among gene clusters. Omics. 16, 284-7.

Journal Pre-proof

Legends Figure 1. Population of muscle fiber types in Mongolian horses. (A) Serial sections of the splenius from a Mongolian horse. The sections were stained with monoclonal antibodies against specific myosin heavy chain (MHC) isoforms. Muscle fibers were classified into Slow twitch phenotype (expressing MHC-I) and Fast twitch phenotype (expressing MHC-II). Bar=100μm. (B) Comparison of the mean value of muscle fiber type

of

population (%) among four different part skeletal muscles of the Mongolian Horse. **: Significant differences in percentage of Slow twitch phenotype when compared with SP. SP:

-p

ro

splenius, TB: triceps brachii, LD: longissimus dorsi, and GM: gluteus medius.

Figure 2. Divergence of RNA expression profiles between splenius and gluteus medius.

re

(A) Volcano plot of adjusted P-value vs. log2(Fold change). The DEGs (padj < 0.05) are

lP

depicted in grey. Eleven selected genes whose expression was validated by qRT-PCR are depicted in other colors. Greater expression in splenius is shown as the log2 (Fold change) > 1

na

and higher expression in gluteus medius as the log2 (Fold change) < -1. (B) GO Enrichment Analysis of DEG between splenius (SP) and gluteus medius (GM). DEGs annotated by the

Jo ur

GO enrichment were ordered by log2 (Fold change).

Figure 3. Differential expression of miRNAs and alternative splicing events of target gene Sox6. (A) Expression levels of DE miRs between splenius (SP) and gluteus medius (GM). (B) Differential expression analysis of Sox6 and myomiRs. In the left picture, horizontal histogram shows reads count numbers of different genes and microRNAs. The right picture shows the difference in expression level between SP and GM. Sox6 showed a reverse expression trend to miR-499 and miR-206. (C) Correlation analysis of gene expression and junction reads of exon skipping in target gene Sox6. R: Pearson correlation coefficient between reads count and junction reads of Sox6. p: p-value of the significance level.

Journal Pre-proof Figure 4. Model of DEGs and myomiRs-mediated regulation of skeletal muscles in the Mongolian horse. MiR-499 and miR-206 involved in muscle performance by regulating the Sox6 (repressor of slow-twitch genes), and miR-499 encoded with the intron of Myh7b, the

Jo ur

na

lP

re

-p

ro

of

exon skipping event happened in exon 10 positively associated with the expression of Sox6.

Journal Pre-proof

Jo ur

na

lP

re

-p

ro

of

Figure 1

Journal Pre-proof

Jo ur

na

lP

re

-p

ro

of

Figure 2

Journal Pre-proof

Jo ur

na

lP

re

-p

ro

of

Figure 3

Journal Pre-proof

Jo ur

na

lP

re

-p

ro

of

Figure 4

Journal Pre-proof Table 1

Significantly enrichment KEGG pathways of DEGs

ID

Description

GeneRatio

p.adjust

ecb00010

Glycolysis / Gluconeogenesis

9/97

8.87E-06

ACSS1/PFKL/LDHB/BPGM/ LDHA/PKM/PGK1/ALDOA/ PFKM

ecb04922

Glucagon signaling pathway

9/97

0.000264

PDE3B/PFKL/PRKACA/LDHB/ PRKAB2/LDHA/PKM/PPP3R1/ PFKM

ecb01230

Biosynthesis of amino acids

7/97

0.001323

PFKL/PKM/ALDH18A1/PGK1/ ASNS/ALDOA/PFKM

ecb04066

HIF-1 signaling pathway

8/97

0.001323

PFKL/LDHB/LDHA/BCL2/PGK1 /PFKFB3/ALDOA/PFKM

9/97

0.001323

Cardiac muscle contraction

7/97

ecb05410

Hypertrophic cardiomyopathy (HCM)

7/97

ecb05414

Dilated cardiomyopathy (DCM)

ecb00051

Fructose and mannose metabolism

ecb00620

Pyruvate metabolism

ecb05230

Central carbon metabolism in cancer

Carbon metabolism

0.001323

MYL2/MYH7/ATP2A2/TPM3/ TNNC1/COX6A2/ATP1A4 MYL2/DES/MYH7/ATP2A2/ TPM3/TNNC1/PRKAB2

0.002232

MYL2/DES/MYH7/PRKACA/ ATP2A2/TPM3/TNNC1

4/97

0.010238

PFKL/PFKFB3/ALDOA/PFKM

4/97

0.01929

ACSS1/LDHB/LDHA/PKM

5/97

0.01929

PFKL/LDHB/LDHA/PKM/PFKM

6/97

0.035142

ACSS1/PFKL/PKM/PGK1/ ALDOA/PFKM

re

-p

0.001926

7/97

lP

na

Jo ur

ecb01200

ro

cardiomyocytes ecb04260

MYL2/MYH7/PRKACA/ATP2A2/ TPM3/TNNC1/ATP2B3/ATP1A4 /BCL2

of

Adrenergic signaling in ecb04261

GeneSymbol

Journal Pre-proof

Jo ur

na

lP

re

-p

ro

of

Graphical abstract

Journal Pre-proof Highlights

of ro -p re lP

4.

na

3.

We examined the skeletal muscle fiber type population muscle of Mongolian horses. We generated deep mRNA sequencing data in fast and slow type muscles of Mongolian horses We found differentially expressed myomiRs (miR-499 and miR-206) and their target gene Sox6. Confirming rMATS analyses, the Sox6 was found an exon skipping in exon 10, which exhibited in both two muscles and show a higher inclusion level in gluteus medius, there was a high correlation between reads count of Sox6 and junction reads of exon skipping.

Jo ur

1. 2.