Next Article in Journal
Effect of 20(S)-Hydroxycholesterol on Multilineage Differentiation of Mesenchymal Stem Cells Isolated from Compact Bones in Chicken
Previous Article in Journal
Segregation Analysis of Rare NRP1 and NRP2 Variants in Families with Lymphedema
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Genes
Volume 11
Issue 11
10.3390/genes11111359
genes-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Gene Expression Profile in Similar Tissues Using Transcriptome Sequencing Data of Whole-Body Horse Skeletal Muscle

by
Ho-Yeon Lee
1,2,
Jae-Yoon Kim
1,2,
Kyoung Hyoun Kim
1,2,
Seongmun Jeong
1,
Youngbum Cho
1,2 and
Namshin Kim
1,2,*
1
Genome Editing Research Center, Korea Research Institute of Bioscience and Biotechnology, Daejeon 34141, Korea
2
Department of Bioinformatics, KRIBB School of Bioscience, University of Science and Technology (UST), Daejeon 34141, Korea
*
Author to whom correspondence should be addressed.
Genes 2020, 11(11), 1359; https://doi.org/10.3390/genes11111359
Submission received: 8 October 2020 / Revised: 11 November 2020 / Accepted: 14 November 2020 / Published: 17 November 2020
(This article belongs to the Section Animal Genetics and Genomics)
Browse Figures
Versions Notes

Abstract

:
Horses have been studied for exercise function rather than food production, unlike most livestock. Therefore, the role and characteristics of tissue landscapes are critically understudied, except for certain muscles used in exercise-related studies. In the present study, we compared RNA-Seq data from 18 Jeju horse skeletal muscles to identify differentially expressed genes (DEGs) between tissues that have similar functions and to characterize these differences. We identified DEGs between different muscles using pairwise differential expression (DE) analyses of tissue transcriptome expression data and classified the samples using the expression values of those genes. Each tissue was largely classified into two groups and their subgroups by k-means clustering, and the DEGs identified in comparison between each group were analyzed by functional/pathway level using gene set enrichment analysis and gene level, confirming the expression of significant genes. As a result of the analysis, the differences in metabolic properties like glycolysis, oxidative phosphorylation, and exercise adaptation of the groups were detected. The results demonstrated that the biochemical and anatomical features of a wide range of muscle tissues in horses could be determined through transcriptome expression analysis, and provided proof-of-concept data demonstrating that RNA-Seq analysis can be used to classify and study in-depth differences between tissues with similar properties.

1. Introduction

Genetic studies of livestock have primarily been aimed at increasing production. Most livestock animals raised today are for meat, and improvements have been made to control fat content and muscle properties to induce rapid muscle growth or maximize the favored types of tissues [1,2,3,4,5,6].
Horses have been identified as a suitable model for studying gene expression in skeletal muscle related to exercise ability [7,8,9,10,11], which is rare for livestock. The present study is similar to a prior report of muscle activity and associated genes in humans [12,13,14,15]. Skeletal muscles have different physiological characteristics depending on their role, which is related to the composition of muscle fibers according to the main purpose of the muscle in question.
The primary component that determines the physiological function of muscles is myosin heavy chain (MyHC), which is myofibril’s thick filament. MyHC is classified into three types, type I, type IIa, and type IIb/IIx, depending on their contraction rate and metabolic phenotype. Type I is a slow oxidative (SO) fiber that has a slow contraction speed and obtains ATP from oxidative phosphorylation (OXPHOS). Type IIb/IIx is a fast glycolytic muscle (FG) that contracts rapidly and take ATP from glycolysis, and type IIa is a fast-oxidative-glycolytic (FOG) fiber with intermediate properties of type I and type IIb/IIx.Type I is a slow oxidative (SO) fiber that has a slow contraction speed and obtains energy from OXPHOS.
Higher content of type I fibers is associated with slower contraction and highly oxidative metabolism. These slow-twitch muscles are reddish due to the high mitochondria, myoglobin, and capillary distribution. The force that can be created by muscle contraction is not strong, but their fatigue resistance makes them used for posture maintenance or continuous and repetitive activities. On the contrary, a higher content of type II (b) fibers is associated with faster contraction and a a glycolytic phenotype [16,17]. Fast-twitch muscles such as these produce ATP through anaerobic glycolysis, and they can exert a strong force in a short time due to their fast-contraction speed but become fatigued easily. Because they use less oxygen, they contain fewer mitochondria, blood vessels, and myoglobin. Consequently, they are less reddish than slow-twitch muscle and look beige or pale-colored. MyHC fibers vary in composition depending on the muscle’s use, and their ratio can be changed by muscle conditioning. Additionally, depending on the characteristics of muscles altered by MyHC content, gene expression associated with the corresponding functions will differ.
In general, to evaluate the distribution, or ratio, of each fiber, measurement of protein expression, or comparison of cross-section thickness to identify the fiber types via microscopy is used [7,16]. However, next-generation sequencing (NGS) analysis of RNA-Seq data can reveal the dominant characteristics of each tissue by comparing global gene expression between tissues, even though the actual fiber composition cannot be directly identified. In addition, it is also possible to compare one tissue repeatedly with tissues from different sites [18].
Most studies of equine muscle tissue are focused on the analysis and improvement of exercise function. The primary breeds used in these studies have been Thoroughbred and Arabian horses, which are used primarily for horse racing or riding, and only a few tissues, such as the blood or gluteus medius, have been evaluated [7,8,9,10,11]. In the present study, we analyzed transcriptome sequencing data from 18 different skeletal muscle tissues taken from the Jeju horse. The Jeju horse is a type of pony that lives on Jeju Island in Korea. It is ~110–120 cm in height and resilient to disease [19].
However, we could not study the unique characteristics of Jeju horse itself because a sufficient number of samples were not secured. Furthermore, because of the specificity of finely classified muscle data from the whole body, we could not obtain similar data from previous studies. Thus, we aimed to identify the transcription expression between skeletal muscle tissues by utilizing the particularity of the subdivided tissue data.
RNA-Seq data of 18 different muscle regions of the Jeju horse were subjected to differential expression (DE) analysis (Figure 1) [20,21], which was classified by k-mean clustering. We used RNA-Seq data to determine if the detailed classification of skeletal muscles with similar properties ware possible. Additionally, we sought to determine the physiological characteristics of each muscle compared with similar and proximal tissues by identifying specific gene expression patterns between tissues.

2. Materials and Methods

2.1. Ethics Statements

All animal experiments were performed in accordance with the guidelines of the Institutional Animal Care and Use Committee and were approved by the Animal Genomics and Bioinformatics Division, National Institute of Animal Science (No.2014-080). All efforts were made to minimize animal suffering.

2.2. Sample Collection

The specimens were obtained from a male horse, born on 12 June 2012, managed by the National Institute of Animal Science, R.D.A, Jeju, South Korea. It was slaughtered through exsanguination after electric stunning. Each tissue was collected from 18 regions of skeletal muscles in a hot-carcass state and immediately frozen to liquid nitrogen and preserved at −80 °C.

2.3. RNA Sequencing

Samples were obtained from a Jeju horse, a region-specific horse breed in Jeju, South Korea. Sample libraries from 18 different skeletal muscles were obtained from one Jeju horse. Ribosomal RNA was removed from total RNA using a RiboMinus Eukaryote kit for RNA-Seq (Thermo Fisher Scientific, Sunnyvale, CA, USA). The RNA-Seq library was prepared using a TruSeq RNA kit (Illumina; San Diego, CA, USA). Sampling and RNA sequencing were conducted by the National Institute of Animal Science, Rural Development Administration, and sequenced using an Illumina HiSeq 2000 (2 × 101 bp). The GEO accession number for this data set is GSE113147.

2.4. Data Processing

Quality assessment was conducted using FastQC version 0.11.5 software (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) to analyze sequence reads in the fastq file format. Using the NGSQCToolkit [22] with a Phred quality score < 20 and less than 50 bp in total length were removed with paired-read and adapter sequences were trimmed.
The reference RNA sequence FASTA files were mapped to the horse genome (EquCab2 79, Equus_caballus.EquCab2.DNA.chromosome.1~31, X, MT, nonchromosomal.fa) downloaded from the Ensembl FTP database. The alignment was performed using the STAR (v 2.5.2b) RNA-Seq aligner with a two-pass method [23]. Using this approach, the first alignment detected splice junctions based on transcript information, and the final alignment was performed using splice junctions as a guide. As a guide for transcript alignment, we used the relevant gtf file (Equus_caballus.EquCab2.87.gtf) as a reference. The alignment process provided read counts for a total of 26,841 genes in 18 samples.

2.5. Differential Expression Analysis

First, the differentially expressed genes (DEGs) between each Jeju horse skeletal muscle were identified by pairwise analysis using the negative binomial test of the DESeq R package [24]. Based on the reads per gene counts identified through sequence mapping, the pairwise DEGs (p-value < 0.05) of all 18 tissues were identified. All DEGs (n = 1292) identified by pairwise analysis and their read counts were used to implement clustering of samples by their differential gene expression.
Clustering was performed using the hclust and the k-means function of the R package, and the optimal k value (k = 2) was calculated by the NbClust package in the R. The k-means clustering divides 18 tissues into large two groups, and we subdivided two groups into three subgroups each. Based on the classified cluster information, we conducted a DE analysis between groups using entire genes (n = 26,841). DESeq2 was used for DE analysis [25].

2.6. IPA Analysis

Further analysis to determine the functions of the identified DEGs was conducted using Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Redwood, CA, USA) [26]. IPA was used to conduct enrichment analysis of the canonical pathway of the gene set and the inter-molecule network. To import the gene list and avoid omissions in the advanced analysis results as much as possible, the Ensembl Equus Caballus gene ID was converted into a well-annotated Ensembl human ID to identify DEGs.

2.7. Gene Ontology

Gene ontology (GO) enrichment analysis of DEGs was performed using the biological processes of the GeneOntology site (http://geneontology.org/, Powered by PANTHER) [27]. GO terms selected only the term of the main category, except subcategories.

3. Results

3.1. Read Alignments and Results

By mapping the raw data to the ensemble horse genome (EquCab2 79) using STAR-2PASS alignment, the raw reads were aligned to a total of 26,841 genes. The number of raw sequence reads per sample ranged from 55.7 to 83.1 M, with an average of 64.1 M and an average input read length of 202 bp (2 × 101 bp). Uniquely mapped reads averaged 53.1 M, representing 83.2% of the total average reads.

3.2. Differential Expression Analysis

3.2.1. Classification

A total of 1292 DEGs were identified through pairwise analysis. This number means that most of the DEGs shown in Figure 2 are overlapping. These are key features that can explain the difference between skeletal muscles among 26,841 genes prepared for clustering 18 samples. We classified samples into hierarchical clustering (using the value of log2 (count + 1)) and k-means clustering through the raw count data of these 1292 DEGs. The optimal k value calculated using the NbClust was set as 2 (Supplementary Materials Figure S1).
Samples were divided primarily into two groups, and the classification results were the same in both hierarchical clustering and k-means clustering (Figure 3). Each group classified by clustering was labeled as “A” or “B” (Figure 3A) and normalized using DESeq2 [25], and subsequently plotted on a PCA plot (Figure 3B). Groups A and B were further divided into subgroups based on the results of k-means clustering (k = 3). The k = 3 is a maximum value that allows two or more tissues to be included in one subgroup for DESeq2 analysis. The clustering results through k-means are presented on the PCA plot of Figure 3B–D.

3.2.2. Comparison of A vs. B

As a result of performing DE analysis on the total gene counts of groups classified into A and B using DESeq2, a total of 1264 DEGs (up = 597, down = 667, adjusted-p < 0.05, |log2Foldchange| > 0.58 (FC = 1.5)) were identified. The expression fold-change for A vs. B in DESeq2 results represents the value of A/B. Subsequent “upregulated/downregulated” means higher/lower expression in the former’s condition of A vs. B. Genes in the Ensemble gene id format were converted into HGNC gene symbol format for further analysis, and a total of 857 (up = 487, down = 370) were successfully converted.
The IPA and gene ontology were used for gene set enrichment analysis of the identified DEGs. Table 1 shows the pathways within the top 10 scores (–log (p-value); > 1.3) with z-score values (not 0) in the IPA canonical pathway results.
GO analysis was performed to elucidate the functional enrichment of upregulated DEGs of biological processes using PANTHER. Table S1 contains the top 10 GO terms based on the fold enrichment value. Directly related to the function of skeletal muscles are “regulation of muscle adaptation (GO:0043502)” and “regulation of striated muscle contraction (GO:0006942)”.

3.2.3. Group A

To identify the characteristics of the muscles in each group, we conducted DE for subgroups in groups A and B, respectively. In A and B, a total of three sets of 1 vs. 2, 1 vs. 3, and 2 vs. 3 were performed. In comparison of group A1 and A2, a total of 236 DEGs (up = 128, down = 108, adjust-P <0.05, |log2FC| > 1) were identified, of which 200 genes (up = 101, down = 99) were converted to the human gene symbol. As a result of analyzing the DEGs by IPA (Table 2), an increase in the glycolysis and gluconeogenesis canonical pathway was confirmed. This is similar in GO (Table S2), as we could find fast-twitch glycolytic muscle-related GO terms like “positive regulation of fast-twitch skeletal muscle fiber contraction (GO:0031448)”, “canonical glycolysis (GO:0061621)”, and “gluconeogenesis (GO:0006094)” in GO analysis of upregulated DEGs. Considering these results, the muscles of group A1 are more glycolytic than the muscles of A2 and are thought to be closer to the fast twitch.
A total of 383 DEGs (up = 3, down = 380, adjust-p < 0.05, |log2FC| > 1) were identified through DE analysis of groups A1 and A3 (301 genes (up = 1, down = 300) were converted to a gene symbol). In A1 vs. A3, all genes were down regulated except three genes, including the only annotated TLE1. The results of the IPA analysis could not identify canonical pathways directly related to metabolic processes, such as glycolysis, and most pathways were associated with immune responses, including T cell signaling pathways (Table 2). Additionally, in GO analysis, which was conducted on downregulated DEGs (Table S3), most of the GO terms related to T cell and immune response were confirmed, and there was no direct result on metabolism.
In group A2 vs. A3, a total of 363 DEGs (up = 98, down = 265, adjust-p < 0.05, |log2FC| > 1) were identified. Among them, 283 were converted to human gene symbols, 89 upregulated DEGs, and 194 downregulated DEGs. In the IPA canonical pathway, downregulation of glycolysis was confirmed, and also included a number of immune-related pathways, such as those found in A1 vs. A3 (Table 2). Likewise, in the GO of downregulated DEGs (Table S3), “canonical glycolysis (GO:0061621)” and “gluconeogenesis (GO:0006094)” were identified, and the GO term related to the immune response occupied the majority as in the IPA analysis. GO terms identified in upregulated DEGs (Table S2) include “cardiac myofibril assembly (GO:0055003)”, “detection of calcium ion (GO:0005513)”, “muscle system process (GO:0003012)”, “locomotion (GO:0040011)”, and “movement of a cell or subcellular component (GO:0006928)”.

3.2.4. Group B

The same as group A, organizations classified as group B were also divided into three subgroups by k-means clustering, and DE analysis was conducted for each group. In comparison with group B1 and group B2, a total of 416 DEGs (up = 184, down = 232, adjust-p <0.05, |log2Foldchanage| > 1) were identified, of which 355 (up = 148, down = 208) has been converted to human gene symbol. The upregulation of glycolysis, oxidative phosphorylation, and gluconeogenesis was found in the IPA canonical pathway for these genes (Table 3). This suggests that in tissues belonging to group B1, overall energy metabolism occurs more actively than muscles in group B2. This is the same in GO analyzed using upregulated gene (Table S4), “canonical glycolysis (GO:0061621)”, “gluconeogenesis (GO:0006094)”, “regulation of oxidative phosphorylation (GO:0002082)”, “mitochondrial electron transport, NADH to ubiquinone (GO:0006120)”, “aerobic respiration (GO:0009060)”, and other GO terms were also identified.
In group B1 vs. group B3, a total of 219 DEGs (up = 151, down = 68, |log2Foldchanage| > 1) were identified, of which 164 (up = 117, down = 46) were converted to gene symbols. In the IPA analysis, a pathway showing a certain tendency to metabolism was not found (Table 3), and pathways presumed to be upregulated are AMPK signaling, senescence pathway, synaptogenesis signaling pathway, factors promoting cardiogenesis in vertebrates, colorectal cancer metastasis signaling, adrenomedullin signaling pathway, cardiac hypertrophy signaling (enhanced), systemic lupus erythematosus in B cell signaling pathway, neuroinflammation signaling pathway, and hepatic fibrosis signaling pathway were identified. In GO of upregulated DEGs (Table S4), “response to nutrient levels (GO:0031667)”, “regulation of lipid metabolic process (GO:0019216)”, “positive regulation of developmental process (GO:0051094)”, and “positive regulation of multicellular organismal process (GO:0051240)” enrichment was confirmed for, and the GO term indicating a certain metabolism type was not identified.
In the DE of group B2 vs. group B3, a total of 164 DEGs (up = 81, down = 83, |log2Foldchanage| > 1) were identified. Of these, 141 (up = 72, down = 69) were converted to human gene symbols. These DEGs showed downregulation of glycolysis and gluconeogenesis in the IPA canonical pathway (Table 3).
The upregulated DEGs showed only association with “muscle contraction (GO:0006936)” in GO of upregulated DEGs (Table S4), but in analysis using downregulated DEGs (Table S3), they also identified the characteristics of the fast-twitch muscle, such as “positive regulation of fast-twitch skeletal muscle fiber contraction (GO:0031448)”, “canonical glycolysis (GO:0061621)”, and “gluconeogenesis (GO:0006094)”.

4. Discussion

In general, it is known that the difference in muscle fiber type and their metabolic properties is the trait that classifies the skeletal muscles. Many studies have demonstrated differences in biochemical-metabolic phenotypes depending on muscle use and fiber composition. In addition, the type and nature of the skeletal muscle fiber that makes up the muscle are transformable and can vary depending on the role the muscle plays. As a result of previous studies of equine exercise capacity, the classification, properties, and related genes of skeletal muscle tissue are well characterized. These prior studies demonstrated a large number of exercise-induced changes, with increased expression of genes associated with oxidative phosphorylation and mitochondrial function in skeletal muscle of individuals endurance-exercised continuously for long periods of time [8,11,15,27]. These results have been verified by molecular biology studies [8,9,15].
In the present study, biochemical analyses using RNA-Seq data from various regions of horse skeletal muscle tissue were conducted, with the objective of identifying whether a classification based on tissue characteristics was possible within very similar tissues. In our 18 skeletal muscle tissue data obtained from Jeju horse, there were no biological or technical replicates for each of the muscle tissues. To compensate for this limitation, the DE between groups was identified by considering the groups of k-means clustering as a replicate for a similar trait.
The DEGs obtained in the pairwise comparison classified 18 muscle tissues into two groups (hierarchical and k-means both), and these were subdivided into three subgroups each. When the expression of genes that determine the fiber type, and the accompanying genes [28], and the expression of genes involved in the functional enrichment of the related traits [29,30,31,32,33] were compared by listing on the basis of clustered samples as a form of the heatmap (Figure 4), it was found that the genes are known to be expressed together with the type of muscle fiber [28], usually increased and decreased here as well. However, the expression of the gene that represents the muscle fiber type and the type of gene that is expected to accompany were not completely consistent.
Among DEGs, the expression of the gene known to be muscle fiber specific also represented various values within the same group A and group B. The difference between the two main groups was more distinguishable in the gene expression of the functional enrichment group of GO. Among the identified DEGs on the heatmap (Figure 4), genes included in “regulation of muscle adaptation (GO:0043502)” were upregulated in group A and downregulated in group B, and genes in the “positive regulation of FA β-oxidation (GO:0032000)” represent similar patterns as well. The functional enrichment for “regulation of muscle adaptation (GO:0043502)” was also identified in the GO results of upregulated DEGs of group A (Table S1). These results suggest that the classification results were preferentially distinguished by differences in muscle adaptation-related functions and gene expressions involved therein.
The canonical pathways identified from group A vs. group B do not directly refer to the effect on skeletal muscle, but many studies have verified their metabolic functions in skeletal muscle and muscle adaptation. As shown in Figure S2, estrogen receptor (ER) signaling exists at the center of functions that DEGs are involved (Table 1). Estrogen is a type of sex hormone that circulates in the body and directly or indirectly affects many molecular pathways through ER, and exercise-induced skeletal muscle adaptation is also affected by ER. ER signaling increases muscle mass regulation and regeneration and is known to enhance ATP production and lipid metabolism by promoting fission fusion and β-oxidation of mitochondria [34,35,36,37,38]. In addition, it was confirmed that a number of canonical pathways included in the table were also affected by the control mechanisms of ER signaling or involved in changes in the types associated with exercise-indicated adaptation [39,40,41,42,43,44,45,46,47,48,49].
In the distribution of muscles shown in Figure 1, the muscles in group A are usually located in the anterior and the muscles in group B are mostly located in the posterior. Taking together these differences in group A vs. group B, this result suggests that the muscles located in the front and rear perform different main functions and thus exhibit different gene expression. Anatomical research of the muscle architecture revealed that the muscle of the forelimbs and hindlimbs have different main roles [21,50]. The proximal horse muscles of the hindlimbs provide energy for locomotion and the muscles of forelimbs act as stiff spring-like struts to support a greater proportion of the body mass [21]. Therefore, the gene expressions for adapting to continuous activity are stronger in the muscles of the forelimbs and anterior that supports the weight of the body. On the other hand, the gene expression for energy generation appears to be more dominant in the muscles of the hindlimbs and the posterior that provide power.
In skeletal muscle adaptation research conducted on porcine, it has been found that adaptation to endurance exercise occurs primarily in forelimb musculature [51]. As shown in Figure 1 and Figure 4, gene expressions for exercise adaptation are dominant in muscles of the anterior and it is thought that the same mechanism will work in equine also. According to Harrison et al. [52], muscle #9 (extensor carpi radialis) located in the forelimbs is activated in the swing state, unlike most of the other muscles in the forelimbs, which are activated in the state of stance. This report could explain why muscle #9 shows a similar expression to the muscles of group B, even though it is located in the anterior.
Unlike general knowledge, the expression of genes directly associated with the skeletal muscle fiber type and with a metabolic tendency (glycolysis, OXPHOS) caused by the fiber type is not unilaterally proportional to the associations of “slow-OXPHOS” and” fast-glycolytic”. It has complex patterns that cannot be explained by a single element. The sub-clusters of groups A and B are sets of tissues that have similarities in these complex expression patterns. The analysis of DEGs between sub-clusters could identify the direct differences of detailed traits, such as fiber-specific gene expression, glycolysis, and oxidative phosphorylation levels (Table 1, Table 2 and Table 3) (Figure 4). Here, the tissues included in the cluster are also mostly located in close or homologous positions (group B3) and are expected to play similar roles. As more tissues are added or replication is secured, more clusters are separated. Thus, differences in simple traits overlap and multiple expression types appear, making it possible to classify a narrower range of locations and functions with only in silico data.

5. Conclusions

In this study, we analyzed the function of muscles with gene expression values obtained from transcriptome sequencing data based on in silico data and literature. We identified that RNA-Seq expression data can be used to classify tissues according to their specific characteristics, even among highly similar tissues, such as different skeletal muscles. Skeletal muscles were categorized by their role and difference in gene expression caused by hormonal and cell signaling pathways. The muscles we analyzed were largely classified into two groups by muscle adaptation-related pathways, which reflected their main roles and location. The two groups of tissues classified can once again form clusters with those with similar properties. This sub-cluster is clustered by detailed types of gene expressions, such as fiber types of muscles or energy metabolic pathways.
If there is an opportunity to confirm the tissue composition and protein expression through analyses using classical methods, this would provide an opportunity to verify the results of the present study and to examine the utility of studying animal tissues with NGS analysis.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4425/11/11/1359/s1, Figure S1: Optimal number of clusters for k-means clustering, Figure S2: IPA Summary network of Group A vs. Group B, Table S1: GO biological process of up-regulated DEGs in Group A vs. Group B, Table S2: GO biological process of up-regulated DEGs in Group A sub-clusters, Table S3: GO biological process of down-regulated DEGs, Table S4: GO biological process of up-regulated DEGs in Group B sub-clusters.

Author Contributions

Conceptualization, H.-Y.L. and N.K.; Methodology, H.-Y.L., J.-Y.K.; Data curation, K.H.K., H.-Y.L.; Data analysis, H.-Y.L.; Writing—original draft preparation, H.-Y.L.; Writing—review and editing, J.-Y.K., K.H.K., S.J., Y.C.; Funding acquisition, N.K.; Supervision, N.K., J.-Y.K., S.J. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by grants from the National Research Foundation of Korea (NRF-2014M3C9A3064552), the KRIBB Initiative program and Agriculture Science & Technology Development (Project No. PJ010453) of Rural Development Administration.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Keel, B.N.; Zarek, C.M.; Keele, J.W.; Kuehn, L.A.; Snelling, W.M.; Oliver, W.T.; Freetly, H.C.; Lindholm-Perry, A.K. Rna-seq meta-analysis identifies genes in skeletal muscle associated with gain and intake across a multi-season study of crossbred beef steers. BMC Genom. 2018, 19, 430. [Google Scholar] [CrossRef]
  2. Kong, B.W.; Hudson, N.; Seo, D.; Lee, S.; Khatri, B.; Lassiter, K.; Cook, D.; Piekarski, A.; Dridi, S.; Anthony, N.; et al. Rna sequencing for global gene expression associated with muscle growth in a single male modern broiler line compared to a foundational barred plymouth rock chicken line. BMC Genom. 2017, 18, 82. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Ropka-Molik, K.; Bereta, A.; Zukowski, K.; Piorkowska, K.; Gurgul, A.; Zak, G. Transcriptomic gene profiling of porcine muscle tissue depending on histological properties. Anim. Sci. J. 2017, 88, 1178–1188. [Google Scholar] [CrossRef] [PubMed]
  4. Li, Y.; Xu, Z.; Li, H.; Xiong, Y.; Zuo, B. Differential transcriptional analysis between red and white skeletal muscle of chinese meishan pigs. Int. J. Biol. Sci. 2010, 6, 350–360. [Google Scholar] [CrossRef] [PubMed]
  5. Francisco, C.L.; Jorge, A.M.; Dal-Pai-Silva, M.; Carani, F.R.; Cabeco, L.C.; Silva, S.R. Muscle fiber type characterization and myosin heavy chain (myhc) isoform expression in mediterranean buffaloes. Meat Sci. 2011, 88, 535–541. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Ryu, Y.C.; Kim, B.C. The relationship between muscle fiber characteristics, postmortem metabolic rate, and meat quality of pig longissimus dorsi muscle. Meat Sci. 2005, 71, 351–357. [Google Scholar] [CrossRef] [PubMed]
  7. Hyytiainen, H.K.; Mykkanen, A.K.; Hielm-Bjorkman, A.K.; Stubbs, N.C.; McGowan, C.M. Muscle fibre type distribution of the thoracolumbar and hindlimb regions of horses: Relating fibre type and functional role. Acta Vet. Scand. 2014, 56, 8. [Google Scholar] [CrossRef] [Green Version]
  8. Ropka-Molik, K.; Stefaniuk-Szmukier, M.; K, Z.U.; Piorkowska, K.; Bugno-Poniewierska, M. Exercise-induced modification of the skeletal muscle transcriptome in arabian horses. Physiol. Genom. 2017, 49, 318–326. [Google Scholar] [CrossRef] [Green Version]
  9. Ropka-Molik, K.; Stefaniuk-Szmukier, M.; Zukowski, K.; Piorkowska, K.; Gurgul, A.; Bugno-Poniewierska, M. Transcriptome profiling of arabian horse blood during training regimens. BMC Genet. 2017, 18, 31. [Google Scholar] [CrossRef] [Green Version]
  10. Serrano, A.L.; Rivero, J.L. Myosin heavy chain profile of equine gluteus medius muscle following prolonged draught-exercise training and detraining. J. Muscle Res. Cell Motil. 2000, 21, 235–245. [Google Scholar] [CrossRef]
  11. Bryan, K.; McGivney, B.A.; Farries, G.; McGettigan, P.A.; McGivney, C.L.; Gough, K.F.; MacHugh, D.E.; Katz, L.M.; Hill, E.W. Equine skeletal muscle adaptations to exercise and training: Evidence of differential regulation of autophagosomal and mitochondrial components. BMC Genom. 2017, 18, 595. [Google Scholar] [CrossRef] [PubMed]
  12. Catoire, M.; Mensink, M.; Boekschoten, M.V.; Hangelbroek, R.; Muller, M.; Schrauwen, P.; Kersten, S. Pronounced effects of acute endurance exercise on gene expression in resting and exercising human skeletal muscle. PLoS ONE 2012, 7, e51066. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Nader, G.A.; von Walden, F.; Liu, C.; Lindvall, J.; Gutmann, L.; Pistilli, E.E.; Gordon, P.M. Resistance exercise training modulates acute gene expression during human skeletal muscle hypertrophy. J. Appl. Physiol. 2014, 116, 693–702. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Popov, D.V.; Lysenko, E.A.; Makhnovskii, P.A.; Kurochkina, N.S.; Vinogradova, O.L. Regulation of PPARGC1A gene expression in trained and untrained human skeletal muscle. Physiol. Rep. 2017, 5, e13543. [Google Scholar] [CrossRef] [Green Version]
  15. Perry, C.G.R.; Hawley, J.A. Molecular Basis of Exercise-Induced Skeletal Muscle Mitochondrial Biogenesis: Historical Advances, Current Knowledge, and Future Challenges. Cold Spring Harb. Perspect. Med. 2017, 8, a029686. [Google Scholar] [CrossRef] [Green Version]
  16. Stuart, C.A.; Stone, W.L.; Howell, M.E.A.; Brannon, M.F.; Hall, H.K.; Gibson, A.L.; Stone, M.H. Myosin content of individual human muscle fibers isolated by laser capture microdissection. Am. J. Physiol. Physiol. 2016, 310, C381–C389. [Google Scholar] [CrossRef] [Green Version]
  17. Schiaffino, S.; Reggiani, C. Fiber Types in Mammalian Skeletal Muscles. Physiol. Rev. 2011, 91, 1447–1531. [Google Scholar] [CrossRef] [Green Version]
  18. Armstrong, E.; Iriarte, A.; Nicolini, P.; Santos, J.D.L.; Ithurralde, J.; Bielli, A.; Bianchi, G.; Peñagaricano, F. Comparison of transcriptomic landscapes of different lamb muscles using RNA-Seq. PLoS ONE 2018, 13, e0200732. [Google Scholar] [CrossRef]
  19. Yoon, S.H.; Kim, J.; Cho, S.; Kwak, W.; Lee, H.-K.; Park, K.D.; Shin, D.; Kim, H. Complete mitochondrial genome sequences of Korean native horse from Jeju Island: Uncovering the spatio-temporal dynamics. Mol. Biol. Rep. 2017, 44, 233–242. [Google Scholar] [CrossRef] [PubMed]
  20. Watson, J.; Wilson, A.M. Muscle architecture of biceps brachii, triceps brachii and supraspinatus in the horse. J. Anat. 2007, 210, 32–40. [Google Scholar] [CrossRef]
  21. Crook, T.C.; Cruickshank, S.E.; McGowan, C.M.; Stubbs, N.; Wakeling, J.M.; Wilson, A.M.; Payne, R.C. Comparative anatomy and muscle architecture of selected hind limb muscles in the Quarter Horse and Arab. J. Anat. 2008, 212, 144–152. [Google Scholar] [CrossRef] [PubMed]
  22. Patel, R.K.; Jain, M. NGS QC toolkit: A toolkit for quality control of next generation sequencing data. PLoS ONE 2012, 7, e30619. [Google Scholar] [CrossRef] [PubMed]
  23. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef] [PubMed]
  24. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11, R106. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for rna-seq data with deseq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Kramer, A.; Green, J.; Pollard, J., Jr.; Tugendreich, S. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics 2014, 30, 523–530. [Google Scholar] [CrossRef] [PubMed]
  27. Thomas, P.D.; Muruganujan, A.; Ebert, D.; Huang, X.; Thomas, P.D. PANTHER version 14: More genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res. 2018, 47, D419–D426. [Google Scholar]
  28. Lynch, C.J.; Xu, Y.; Hajnal, A.; Salzberg, A.C.; Kawasawa, Y.I. RNA Sequencing Reveals a Slow to Fast Muscle Fiber Type Transition after Olanzapine Infusion in Rats. PLoS ONE 2015, 10, e0123966. [Google Scholar] [CrossRef]
  29. Musumeci, O.; Bruno, C.; Mongini, T.; Rodolico, C.; Aguennouz, M.; Barca, E.; Amati, A.; Cassandrini, D.; Serlenga, L.; Vita, G.; et al. Clinical features and new molecular findings in muscle phosphofructokinase deficiency (GSD type VII). Neuromuscul. Disord. 2012, 22, 325–330. [Google Scholar] [CrossRef]
  30. Scott, W.; Stevens, J.; A Binder–Macleod, S. Human Skeletal Muscle Fiber Type Classifications. Phys. Ther. 2001, 81, 1810–1816. [Google Scholar] [CrossRef]
  31. Peshavaria, M.; Day, I.N. Molecular structure of the human muscle-specific enolase gene (ENO3). Biochem. J. 1991, 275, 427–433. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Tsujino, S.; Sakoda, S.; Mizuno, R.; Kobayashi, T.; Suzuki, T.; Kishimoto, S.; Shanske, S.; DiMauro, S.; A Schon, E. Structure of the gene encoding the muscle-specific subunit of human phosphoglycerate mutase. J. Biol. Chem. 1989, 264, 15334–15337. [Google Scholar] [PubMed]
  33. Wei, B.; Jin, J.-P. TNNT1, TNNT2, and TNNT3: Isoform genes, regulation, and structure–function relationships. Gene 2016, 582, 1–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Ikeda, K.; Horie-Inoue, K.; Inoue, S. Functions of estrogen and estrogen receptor signaling on skeletal muscle. J. Steroid Biochem. Mol. Biol. 2019, 191, 105375. [Google Scholar] [CrossRef] [PubMed]
  35. Brown, E.L.; Hazen, B.C.; Eury, E.; Wattez, J.-S.; Gantner, M.L.; Albert, V.; Chau, S.; Sanchez-Alavez, M.; Conti, B.; Kralli, A. Estrogen-Related Receptors Mediate the Adaptive Response of Brown Adipose Tissue to Adrenergic Stimulation. iScience 2018, 2, 221–237. [Google Scholar] [CrossRef] [Green Version]
  36. Frank, A.P.; Palmer, B.F.; Clegg, D.J. Do estrogens enhance activation of brown and beiging of adipose tissues? Physiol. Behav. 2018, 187, 24–31. [Google Scholar] [CrossRef]
  37. Sánchez-Aguilera, A.; Arranz, L.; Martín-Pérez, D.; García-García, A.; Stavropoulou, V.; Kubovcakova, L.; Isern, J.; Martín-Salamanca, S.; Langa, X.; Skoda, R.C.; et al. Estrogen Signaling Selectively Induces Apoptosis of Hematopoietic Progenitors and Myeloid Neoplasms without Harming Steady-State Hematopoiesis. Cell Stem Cell 2014, 15, 791–804. [Google Scholar] [CrossRef] [Green Version]
  38. Radovick, S.; Levine, J.E.; Wolfe, A.M. Estrogenic Regulation of the GnRH Neuron. Front. Endocrinol. 2012, 3, 52. [Google Scholar] [CrossRef] [Green Version]
  39. Lefeuvre, R.A.; Rothwell, N.J.; Stock, M.J. Activation of brown fat thermogenesis in response to central injection of corticotropin releasing hormone in the rat. Neuropharmacology 1987, 26, 1217–1221. [Google Scholar] [CrossRef]
  40. Schiaffino, S.; Sandri, M.; Murgia, M. Activity-Dependent Signaling Pathways Controlling Muscle Diversity and Plasticity. Physiology 2007, 22, 269–278. [Google Scholar] [CrossRef] [Green Version]
  41. Solinas, G.; Summermatter, S.; Mainieri, D.; Gubler, M.; Montani, J.P.; Seydoux, J.; Smith, S.R.; Dulloo, A.G. Corticotropin-Releasing Hormone Directly Stimulates Thermogenesis in Skeletal Muscle Possibly through Substrate Cycling between de Novo Lipogenesis and Lipid Oxidation. Endocrinology 2006, 147, 31–38. [Google Scholar] [CrossRef] [PubMed]
  42. Ubaida-Mohien, C.; Lyashkov, A.; Gonzalez-Freire, M.; Tharakan, R.; Shardell, M.; Moaddel, R.; Semba, R.D.; Chia, C.W.; Gorospe, M.; Rosen, C.J.; et al. Discovery proteomics in aging human skeletal muscle finds change in spliceosome, immunity, proteostasis and mitochondria. eLife 2019, 8, 8. [Google Scholar] [CrossRef] [PubMed]
  43. Chen, M.; Feng, H.Z.; Gupta, D.; Kelleher, J.; Dickerson, K.E.; Wang, J.; Hunt, D.; Jou, W.; Gavrilova, O.; Jin, J.P.; et al. G(s)α deficiency in skeletal muscle leads to reduced muscle mass, fiber-type switching, and glucose intolerance without insulin resistance or deficiency. Am. J. Physiol. Cell Physiol. 2009, 296, C930–C940. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Hood, D.A.; Memme, J.M.; Oliveira, A.N.; Triolo, M. Maintenance of Skeletal Muscle Mitochondria in Health, Exercise, and Aging. Annu. Rev. Physiol. 2019, 81, 19–41. [Google Scholar] [CrossRef] [PubMed]
  45. Roth, L.; Koncina, E.; Satkauskas, S.; Crémel, G.; Aunis, D.; Bagnard, D. The many faces of semaphorins: From development to pathology. Cell. Mol. Life Sci. 2008, 66, 649–666. [Google Scholar] [CrossRef]
  46. Ubaida-Mohien, C.; Gonzalez-Freire, M.; Lyashkov, A.; Moaddel, R.; Chia, C.W.; Simonsick, E.M.; Sen, R.; Ferrucci, L. Physical Activity Associated Proteomics of Skeletal Muscle: Being Physically Active in Daily Life May Protect Skeletal Muscle From Aging. Front. Physiol. 2019, 10, 312. [Google Scholar] [CrossRef]
  47. Xiang, W.; Zhang, B.; Lv, F.; Feng, G.; Chen, L.; Yang, F.; Zhang, K.; Cao, C.; Wang, P.; Chu, M. The potential regulatory mechanisms of the gonadotropin-releasing hormone in gonadotropin transcriptions identified with bioinformatics analyses. Reprod. Biol. Endocrinol. 2017, 15, 46. [Google Scholar] [CrossRef]
  48. Krzeminski, K. The Role of Adrenomedullin in Cardiovascular Response to Exercise—A Review. J. Hum. Kinet. 2016, 53, 127–142. [Google Scholar] [CrossRef] [Green Version]
  49. Shimizu, I.; Yoshida, Y.; Moriya, J.; Nojima, A.; Uemura, A.; Kobayashi, Y.; Minamino, T. Semaphorin3E-Induced Inflammation Contributes to Insulin Resistance in Dietary Obesity. Cell Metab. 2013, 18, 491–504. [Google Scholar] [CrossRef] [Green Version]
  50. Rivero, J.-L.L. A Scientific Background for Skeletal Muscle Conditioning in Equine Practice. J. Vet. Med. Ser. A 2007, 54, 321–332. [Google Scholar] [CrossRef]
  51. McAllister, R.M.; Reiter, B.L.; Amann, J.F.; Laughlin, M.H. Skeletal muscle biochemical adaptations to exercise training in miniature swine. J. Appl. Physiol. 1997, 82, 1862–1868. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Harrison, S.M.; Whitton, R.C.; King, M.; Haussler, K.K.; Kawcak, C.E.; Stover, S.M.; Pandy, M.G. Forelimb muscle activity during equine locomotion. J. Exp. Biol. 2012, 215, 2980–2991. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Genes 11 01359 g001 550
Figure 1. Name and location of 18 skeletal muscle tissues of Equus caballus.
Figure 1. Name and location of 18 skeletal muscle tissues of Equus caballus.
Genes 11 01359 g001
Genes 11 01359 g002 550
Figure 2. Number of DEGs (Differentially expressed genes). Number of DEGs per tissues identified through pairwise tests of 18 skeletal muscles (min = 20, max = 198), the sample ID is the same as shown in Figure 1.
Figure 2. Number of DEGs (Differentially expressed genes). Number of DEGs per tissues identified through pairwise tests of 18 skeletal muscles (min = 20, max = 198), the sample ID is the same as shown in Figure 1.
Genes 11 01359 g002
Genes 11 01359 g003 550
Figure 3. Classification of samples by gene counts. The results of k-means clustering are shown on hierarchical clustering and PCA (Principal component analysis) plot. (A) The result of hierarchical clustering using read counts of total DEGs was the same as K-means clustering (k = 2). (B) Clustering results for all 18 groups on PCA plot. (C) Subgroups in group A by k-means clustering on the PCA plot. (D) Subgroups in group B by k-means clustering on the PCA plot.
Figure 3. Classification of samples by gene counts. The results of k-means clustering are shown on hierarchical clustering and PCA (Principal component analysis) plot. (A) The result of hierarchical clustering using read counts of total DEGs was the same as K-means clustering (k = 2). (B) Clustering results for all 18 groups on PCA plot. (C) Subgroups in group A by k-means clustering on the PCA plot. (D) Subgroups in group B by k-means clustering on the PCA plot.
Genes 11 01359 g003
Genes 11 01359 g004 550
Figure 4. Expression heatmaps in each group and sample for slow-twitch- and fast-twitch-specific genes. The list of genes is based on the (i) fiber-type-specific genes from the literature base (TNNC2 and TNNI1 are not included because they do not have annotation ids in the reference), (ii) DEGs in IPA canonical pathways and GO, (iii) GO database (positive regulation of FA β-oxidation). The color of the heatmap for gene counts is blue for the lower percentile than 50 and red for higher than 50 (for each row).
Figure 4. Expression heatmaps in each group and sample for slow-twitch- and fast-twitch-specific genes. The list of genes is based on the (i) fiber-type-specific genes from the literature base (TNNC2 and TNNI1 are not included because they do not have annotation ids in the reference), (ii) DEGs in IPA canonical pathways and GO, (iii) GO database (positive regulation of FA β-oxidation). The color of the heatmap for gene counts is blue for the lower percentile than 50 and red for higher than 50 (for each row).
Genes 11 01359 g004
Table
Table 1. Top 10 canonical pathway in IPA analysis of group A vs. group B.
Table 1. Top 10 canonical pathway in IPA analysis of group A vs. group B.
Ingenuity Canonical Pathways-log
(p-Value)		z-ScoreMolecules
Estrogen Receptor Signaling4.11.569ADCY1, ADCY5, ADCY9, CACNA1C, CARM1, CREB5, CREBBP, DDX5, FBXO32, FOXO4, GNAZ, GPS2, IGF2R, MAP2K2, MED4, MMP14, MMP15, MMP16, MMP2, NCOR2, NOS3, NOTCH1, PIK3CB, PLCB3, PRKACA, TRIM63
Endocannabinoid Cancer Inhibition Pathway3.89−1.807ADCY1, ADCY5, ADCY9, CREB5, CREBBP, DDIT3, MAP2K2, MAP2K7, MMP2, NOS1, NOS2, NOS3, NUPR1, PIK3CB, PRKACA
Semaphorin Neuronal Repulsive Signaling Pathway3.48−0.535CSPG4, ITGA3, MAP2K2, MAP2K7, MAPT, PAK4, PDE4A, PIK3CB, PLXNA1, PLXNA2, PLXND1, PRKACA, SEMA3E, SMC3
GNRH Signaling3.471.941ADCY1, ADCY5, ADCY9, CACNA1C, CACNA1G, CACNA1H, CREB5, CREBBP, HBEGF, MAP2K2, MAP2K7, MAP3K11, MMP2, PAK4, PLCB3, PRKACA
Corticotropin Releasing Hormone Signaling3.31.387ADCY1, ADCY5, ADCY9, CACNA1C, CACNA1G, CACNA1H, CREB5, CREBBP, MAP2K2, NOS1, NOS2, NOS3, PRKACA, SLC39A7
Gαs Signaling2.932.111ADCY1, ADCY5, ADCY9, ADD3, ADRB2, CREB5, CREBBP, MAP2K2, PRKACA, RAPGEF2, RYR1
Spliceosomal Cycle2.92−2.646DDX46, DHX15, MAGOH, PRPF18, RBM8A, SLU7, ZNF830
Adrenomedullin signaling pathway2.862.673ADCY1, ADCY5, ADCY9, KCNH2, KCNN3, MAP2K2, MAP2K7, MMP2, MYLK2, NOS3, PIK3CB, PLCB3, PRKACA, RAMP2, RXRA, SLC39A7
White Adipose Tissue Browning Pathway2.772.887ADCY1, ADCY5, ADCY9, CACNA1C, CACNA1G, CACNA1H, CREB5, CREBBP, LDHD, PPARA, PRKACA, RXRA
Calcium Signaling2.661.897ATP2A2, CACNA1C, CACNA1G, CACNA1H, CHRNG, CREB5, CREBBP, MICU1, MYH14, MYH7, MYH8, PPP3CB, PRKACA, RCAN1, RYR1, TP63
Table
Table 2. Top 10 IPA canonical pathways in group A.
Table 2. Top 10 IPA canonical pathways in group A.
Ingenuity Canonical Pathways-log
(p-Value)		z-ScoreMolecules
A1 vs. A2
Glycolysis I141.897ALDOA, ALDOC, ENO3, GPI, PFKL, PFKM, PGAM2, PGK1, PKM, TPI1
Gluconeogenesis I7.111.633ALDOA, ALDOC, ENO3, GPI, PGAM2, PGK1
Calcium Signaling5.711ATP2A1, ATP2B2, CASQ2, CREB3L4, MYH1, MYH11, MYL1, MYL6B, TNNI2, TNNT3, TPM1
Actin Cytoskeleton Signaling5.471.667ACTN3, DIAPH3, EGF, FGF9, HRAS, LIMK1, MYH1, MYH11, MYL1, MYL6B, MYLPF
Protein Kinase A Signaling4.961.265CREB3L4, MYL1, MYL6B, MYLPF, NAPEPLD, PGP, PHKB, PLCL1, PLCL2, PPP1R3D, PTPN3, PYGM, TNNI2, UBASH3B
Estrogen Receptor Signaling3.841.508BCL2, CREB3L4, EGF, HRAS, LIMK1, MYL1, MYL6B, MYLPF, PLCL1, PLCL2, SETD7
Apelin Cardiomyocyte Signaling Pathway3.651.633ATP2A1, MYL1, MYL6B, MYLPF, PLCL1, PLCL2
Synaptic Long Term Potentiation3.030.816CREB3L4, HRAS, PLCL1, PLCL2, PPP1R1A, PPP1R3D
Semaphorin Neuronal Repulsive Signaling Pathway3.020.816DPYSL2, LIMK1, MYL1, MYL6B, MYLPF, VCAN
PAK Signaling2.81.342HRAS, LIMK1, MYL1, MYL6B, MYLPF
A1 vs. A3
iCOS-iCOSL Signaling in T Helper Cells24.8−4.359CAMK4, CD247,CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, FCER1G, HLA-DOA, HLA-DOB, HLA-DRA, ICOS, IKBKE, IL2RA, IL2RB, IL2RG, ITK, LAT, LCK, LCP2, PIK3CG, PTPRC, VAV1, ZAP70
CD28 Signaling in T Helper Cells23.8−3.771ARPC1B, CAMK4, CARD11, CD247, CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, FCER1G, HLA-DOA, HLA-DOB, HLA-DRA, IKBKE, ITK, LAT, LCK, LCP2, PIK3CG, PTPN6, PTPRC, SYK, VAV1, ZAP70
Th2 Pathway19.8−3.638CCR1, CD247, CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, CXCR4, HLA-DOA, HLA-DOB, HLA DRA, ICOS, IKZF1, IL2RA, IL2RB, IL2RG, ITGB2, JAK3, PIK3CG, SPI1, TIMD4, VAV1
Th1 Pathway17.2−3.638CD247, CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, CD8A, CXCR3, HLA-DOA, HLA-DOB, HLA-DRA, ICOS, IL10RA, IL18R1, IRF1, ITGB2, JAK3, PIK3CG, VAV1
PKCθ Signaling in T Lymphocytes16.1−4.69CARD11, CD247, CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, FCER1G, HLA-DOA, HLA-DOB, HLA-DRA, IKBKE, LAT, LCK, LCP2, PIK3CG, RAC2, VAV1, VAV3, ZAP70
Role of NFAT in Regulation of the Immune Response15.7−4.796CAMK4, CD247, CD28, CD3D, CD3E, CD3G, CD4, CD80, CD86, FCER1G, FCGR2C, HLA-DOA, HLA-DOB, HLA-DRA, IKBKE, ITK, LAT, LCK, LCP2, PIK3CG, PLCB2, SYK, ZAP70
Calcium-induced T Lymphocyte Apoptosis12.9−3.606ATP2A3, CAMK4, CD247, CD3D, CD3E, CD3G, CD4, FCER1G, HLA-DOA, HLA-DOB, HLA-DRA, LCK, PRKCB, ZAP70
PD-1, PD-L1 cancer immunotherapy pathway11.13.873CD247, CD28, CD80, HLA-DOA, HLA-DOB, HLA-DRA, IL2RA, IL2RB, IL2RG, JAK3, LAT, LCK, LCP2, PIK3CG, ZAP70
Type I Diabetes Mellitus Signaling10.8−2.828CASP8, CD247, CD28, CD3D, CD3E, CD3G, CD80, CD86, FCER1G, HLA-DOA, HLA-DOB, HLA-DRA, IKBKE, IRF1, PRF1
B Cell Receptor Signaling10.4−3.051APBB1IP, CAMK4, DAPP1, FCGR2C, IGHE, IGHG4, IGHM, IKBKE, PIK3AP1, PIK3CG, PRKCB, PTK2B, PTPN6, PTPRC, RAC2, SYK, VAV1, VAV3
A2 vs. A3
iCOS-iCOSL Signaling in T Helper Cells12.3−3.464CD3D, CD3E, CD4, FCER1G, HLA-DOA, ICOS, IL2RA, IL2RG, ITK, LCK, LCP2, PIK3CD, PIK3CG, PTPRC, VAV1, ZAP70
B Cell Receptor Signaling10.8−2.496CD22, CREB3L4, DAPP1, FCGR2C, IGHE, IGHG4, IGHM, MAP2K6, PIK3AP1, PIK3CD, PIK3CG, PLCG2, PRKCB, PTPRC, RAC2, RASSF5, SYNJ2, VAV1
Phospholipase C Signaling10.1−2.138CD3D, CD3E, CREB3L4, FCER1G, FCGR2C, IGHG4, ITGA4, ITK, LCK, LCP2, MYL2, MYL6B, MYLPF, NAPEPLD, PLCB2, PLCG2, PLD4, PRKCB, RAC2, ZAP70
Glycolysis I9.21−1.414ALDOA, ALDOC, ENO3, PFKL, PFKM, PGAM2, PGK1, PKM
Actin Cytoskeleton Signaling8.71−1.604EGF, FGD3, FGF9, ITGA4, LIMK1, MYH1, MYH10, MYL2, MYL6B, MYLPF, NCKAP1L, PIK3CD, PIK3CG, RAC2, TIAM2, TMSB4Y, VAV1
CD28 Signaling in T Helper Cells8.51−3CD3D, CD3E, CD4, FCER1G, HLA-DOA, ITK, LCK, LCP2, PIK3CD, PIK3CG, PTPRC, VAV1, ZAP70
Calcium-induced T Lymphocyte Apoptosis8.11−3ATP2A1, ATP2A3, CD3D, CD3E, CD4, FCER1G, HLA-DOA, LCK, PRKCB, ZAP70
Table
Table 3. Top 10 IPA canonical pathway in group B.
Table 3. Top 10 IPA canonical pathway in group B.
Ingenuity Canonical Pathways-log
(p-Value)		z-ScoreMolecules
B1 vs. B2
Actin Cytoskeleton Signaling6.570.632ACTN3, ARHGAP24, EGF, FGF1, FGF10, FGF7, FGF9, HRAS, MYH1, MYH10, MYH3, MYL6B, MYLPF, PAK1, PFN2, TIAM2
Glycolysis I5.672.449ALDOA, ENO3, GPI, PGAM2, PGK1, PKM
Oxidative Phosphorylation5.173.162COX4I1, COX5A, COX7A1, CYB5A, MT-ND4L, NDUFA8, NDUFS7, NDUFS8, NDUFV1, UQCR11
Gluconeogenesis I4.392.236ALDOA, ENO3, GPI, PGAM2, PGK1
Sirtuin Signaling Pathway3.78−1.732ARG2, IDH2, LDHA, MT-ND4L, NDUFA8, NDUFS7, NDUFS8, NDUFV1, PFKFB3, PGAM2, PGK1, PPIF, SREBF1, TUBA8
Calcium Signaling3.46−1CAMKK2, CASQ2, CREB3L4, MYH1, MYH10, MYH3, MYL6B, SLC8A3, TNNI2, TNNT3, TPM1
FGF Signaling2.74−0.816CREB3L4, FGF1, FGF10, FGF7, FGF9, HRAS
Neuregulin Signaling2.450.447BTC, EGF, ERBIN, HRAS, RNF41, STAT5A
Bladder Cancer Signaling2.43−0.447EGF, FGF1, FGF10, FGF7, FGF9, HRAS
LPS/IL-1 Mediated Inhibition of RXR Function1.660.0357ABCB9, ACOX3, APOE, FABP4, MAOB, PLTP, SREBF1, TNFRSF11B
B1 vs. B3
AMPK Signaling3.340.816ACACB, ADIPOQ, EIF4EBP1, PCK2, PFKFB3, PIK3R6, TBC1D1
Senescence Pathway2.071.633ACVR1C, DHCR24, E2F8, EIF4EBP1, PIK3R6, TGFB3
Synaptogenesis Signaling Pathway1.821.633APOE, CDH15, EIF4EBP1, GSK3B, PIK3R6, SNCG
Factors Promoting Cardiogenesis in Vertebrates1.82ACVR1C, BMPR1B, GSK3B, TGFB3
Colorectal Cancer Metastasis Signaling1.632BAX, GSK3B, PIK3R6, PTGER3, TGFB3
Adrenomedullin signaling pathway1.421BAX, GSK3B, GUCY2C, PIK3R6
B2 vs. B3
Glycolysis I11.6−2.121ALDOA, ALDOC, ENO3, GPI, PFKM, PGAM2, PGK1, PKM
Calcium Signaling8.23−0.447ATP2A1, CAMKK2, CASQ2, CREB3L4, MYH1, MYH3, MYL2, MYL6B, PPP3CA, TNNI2, TNNT3, TPM3
Gluconeogenesis I7.99−1.633ALDOA, ALDOC, ENO3, GPI, PGAM2, PGK1
Semaphorin Neuronal Repulsive Signaling Pathway4.79−0.378DPYSL2, DPYSL3, MYL2, MYL6B, PAK1, PDE4D, PLXNA3
Actin Cytoskeleton Signaling4.2−1.134ACTN3, EGF, MYH1, MYH3, MYL2, MYL6B, PAK1, SSH2
PFKFB4 Signaling Pathway3.741CREB3L4, GPI, PFK M, TGFB3
HIF1α Signaling2.76−1.633ADM, EGF, GPI, LDHA, PKM, PPP3CA
Colanic Acid Building Blocks Biosynthesis2.49#NUM!GPI, UGP2
White Adipose Tissue Browning Pathway2.081CAMKK2, CREB3L4, FNDC5, LDHA
AMPK Signaling1.971CAMKK2, CREB3L4, GYS1, PFKM, ULK1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lee, H.-Y.; Kim, J.-Y.; Kim, K.H.; Jeong, S.; Cho, Y.; Kim, N. Gene Expression Profile in Similar Tissues Using Transcriptome Sequencing Data of Whole-Body Horse Skeletal Muscle. Genes 2020, 11, 1359. https://doi.org/10.3390/genes11111359

AMA Style

Lee H-Y, Kim J-Y, Kim KH, Jeong S, Cho Y, Kim N. Gene Expression Profile in Similar Tissues Using Transcriptome Sequencing Data of Whole-Body Horse Skeletal Muscle. Genes. 2020; 11(11):1359. https://doi.org/10.3390/genes11111359

Chicago/Turabian Style

Lee, Ho-Yeon, Jae-Yoon Kim, Kyoung Hyoun Kim, Seongmun Jeong, Youngbum Cho, and Namshin Kim. 2020. "Gene Expression Profile in Similar Tissues Using Transcriptome Sequencing Data of Whole-Body Horse Skeletal Muscle" Genes 11, no. 11: 1359. https://doi.org/10.3390/genes11111359

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop