ORIGINAL RESEARCH article

Front. Vet. Sci., 13 December 2021

Sec. Livestock Genomics

Volume 8 - 2021 | https://doi.org/10.3389/fvets.2021.794628

Screening and Identification of Muscle-Specific Candidate Genes via Mouse Microarray Data Analysis

  • SH

    Sayed Haidar Abbas Raza 1

  • CL

    Chengcheng Liang 1

  • WG

    Wang Guohua 1

  • SD

    Sameer D. Pant 2

  • ZM

    Zuhair M. Mohammedsaleh 3

  • AF

    Abdullah F. Shater 3

  • MA

    Mashael Alhumaidi Alotaibi 4

  • RK

    Rajwali Khan 5

  • NS

    Nicola Schreurs 6

  • GC

    Gong Cheng 1

  • CM

    Chugang Mei 1

  • LZ

    Linsen Zan 1,7*

  • 1. College of Animal Science and Technology, Northwest A&F University, Yangling, China

  • 2. School of Animal & Veterinary Sciences, Charles Sturt University, Wagga Wagga, NSW, Australia

  • 3. Department of Medical Laboratory Technology, Faculty of Applied Medical Sciences, University of Tabuk, Tabuk, Saudi Arabia

  • 4. Department of Biology, College of Science, Jouf University, Sakaka, Saudi Arabia

  • 5. Department of Livestock Management, Breeding and Genetic, The University of Agriculture Peshawar, Peshawar, Pakistan

  • 6. Animal Science, School of Agriculture and Environment, Massey University, Palmerston North, New Zealand

  • 7. National Beef Cattle Improvement Center, Northwest A&F University, Yangling, China

Abstract

Muscle tissue is involved with every stage of life activities and has roles in biological processes. For example, the blood circulation system needs the heart muscle to transport blood to all parts, and the movement cannot be separated from the participation of skeletal muscle. However, the process of muscle development and the regulatory mechanisms of muscle development are not clear at present. In this study, we used bioinformatics techniques to identify differentially expressed genes specifically expressed in multiple muscle tissues of mice as potential candidate genes for studying the regulatory mechanisms of muscle development. Mouse tissue microarray data from 18 tissue samples was selected from the GEO database for analysis. Muscle tissue as the treatment group, and the other 17 tissues as the control group. Genes expressed in the muscle tissue were different to those in the other 17 tissues and identified 272 differential genes with highly specific expression in muscle tissue, including 260 up-regulated genes and 12 down regulated genes. is the genes were associated with the myofibril, contractile fibers, and sarcomere, cytoskeletal protein binding, and actin binding. KEGG pathway analysis showed that the differentially expressed genes in muscle tissue were mainly concentrated in pathways for AMPK signaling, cGMP PKG signaling calcium signaling, glycolysis, and, arginine and proline metabolism. A PPI protein interaction network was constructed for the selected differential genes, and the MCODE module used for modular analysis. Five modules with Score > 3.0 are selected. Then the Cytoscape software was used to analyze the tissue specificity of differential genes, and the genes with high degree scores collected, and some common genes selected for quantitative PCR verification. The conclusion is that we have screened the differentially expressed gene set specific to mouse muscle to provide potential candidate genes for the study of the important mechanisms of muscle development.

Introduction

Muscles represent a crucial group of soft tissues derived from the mesoderm that are primarily responsible for locomotion, and movement in all animals butthe World Health Organization estimates musculoskeletal disorders cause the highest proportion of disabilities worldwide, affecting approximately 1.7 billion people (). Therefore, there is significant interest in characterizing the genetics that underpin muscular development, and any associated pathophysiology.

There are three major group of muscles i.e., skeletal, myocardium and smooth muscles. Skeletal muscles weigh about 40% of adult weight in humans, and represent the main subgroup of muscles that allow for locomotion in conjunction with the skeletal system. Apart from locomotion, skeletal muscles also have other important functions e.g., heat production, support and protection of other soft tissues, and participation in metabolic homeostasis (, ). Diseases that affect primary skeletal muscles or the neuromuscular junction frequently manifest in the form of pathological muscle weakness or reduced skeletal muscle mass, which weakens the body's ability to respond to stress and chronic diseases (). Moreover, amino acids released from muscles help maintain blood sugar levels during starvation. Therefore, diseases affecting skeletal muscles can result in wide ranging pathologies, and represent a key cause of morbidity and disability in human populations.

Skeletal muscles are multinucleated, and develop via the fusion of myogenic progenitor cells called myoblasts, into muscle fibers called myotubes, via a complex process known as myogenesis (, ). Several genes are known to play a crucial role either during myogenesis, or subsequently, in ensuring normal muscle physiology (). Some of the main genes involved in muscular development include transcription factors MYOD1 (myogenic differentiation 1), MYF5 (myogenic factor 5), MYOG (myogenin) and MRF (myogenic regulatory factor), MYF6 (herculin), PAX3 (paired box 3), PAX7 (paired box 7) and MEF2 (myocyte enhancer factor 2) family (). MYOD1 and MYF5 are involved in the early phases of skeletal muscle development by promoting the proliferation and differentiation of myogenic progenitor cells into myoblasts, while MYOG plays an important role in the latter phases of myogenesis that involve fusion of myoblasts into myotubes. The precise function of MYF6 remains unknown, though it is thought to regulate myogenesis, and is exclusively expressed in skeletal muscles ().

Apart from these widely known genes, several other genes that influence either skeletal muscle development or physiology remain unidentified and/or uncharacterized.

High-throughput gene chip technologies that provide large-scale gene expression data by measuring transcript abundance in various tissues or cells (), can be leveraged in combination with online gene expression databases (e.g., NCBI's GEO database) to identify such genes. Moreover, given that human populations are genetically heterogenous, inbred animals models can be very useful in identifying and characterizing key genes associated with muscle development and disease.

Therefore, the overall aim of the present study, was to identify genes that are differentially expressed in skeletal muscles of 10–12 week old C57BL/6 mice, by comparing skeletal muscle expression profiles against 16 non-muscle tissues. Genes identified as differentially expressed in muscles, were subsequently subjected to bioinformatic analyses including process and pathway enrichment analysis, protein-protein interaction (PPI) network construction and molecular compounding. Finally, the genes with partial height difference multiples were selected for validation via qPCR.

Materials and Methods

Ethics Statement

All procedures were approved by the Experimental Animal Center of Xi'an Jiaotong University. Animal care and use protocols (EACXU 172) were approved by the Institutional Animal Care and Use Committee of Xi'an Jiaotong University and Northwest A&F University, Yangling. All animal experiments were performed in adherence with the NIH Guidelines on the Use of Laboratory Animals.

Microarray Data

The microarray data was downloaded from NCBI's GEO (Gene Expression Omnibus) database (GEO accession number GSE9954). The downloaded dataset contained microarray expression data from 70 samples that collectively represent 22 tissues (including muscles). The microarray dataset was derived from 10–12 week old male C57BL/6 mice using the Affymetrix Mouse Genome 430 2.0 Array platform (GPL1261). After euthanasia, multiple organs and tissues were taken for microarray analysis (). In this study, microarray data from 18 out of the 22 available tissues were selected. The selected tissues included muscles, adipose tissue, adrenal gland, bone marrow, brain, eye, heart, kidney, liver, lung, pituitary gland, placenta, salivary gland, seminal vesicle, small intestine, spleen, testis, and thymus.

Differential gene expression analysis was performed on the downloaded microarray data using the R project for statistical computing (version 3.5.2; https://www.r-project.org/) packages “limma” (http://www.bioconductor.org/packages/3.5/bioc/html/limma.html) (), and “impute” (http://www.bioconductor.org/packages/2.7/bioc/html/impute.html) (). Screeningfor differentially expressed genes (DEGs) was performed by comparing mouse muscle expression profiles against the remaining 17 tissues using a P-value threshold of <0.05, and log2(fold change) threshold of ≥2.

Process and Pathway Enrichment Analyses

Genes identified as differentially expressed in the initial screening, were subjected to Gene Ontology (GO) analysis via the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/home.jsp) using the Mus musculus genome annotation as background. Three aspects of the GO database were targeted in the GO enrichment analyses i.e., cellular component (CC), molecular function (MF), and biological process (BP). Similarly, KEGG pathway enrichment analysis was performed using DAVID and KOBAS (KEGG Orthology-Based Annotation System–http://kobas.cbi.pku.edu.cn/). The R package “ggplot2” (version: 3.1.0; http://ggplot2.tidyverse.org) was used for data visualization.

PPI Network and Module Analysis

Protein-protein interaction network and module analysis was performed using online tools and the String database (https://string-db.org/). Genes identified as differentially expressed in muscles were used to construct a PPI network map (), and the MCODE (Molecular Complex Detection) plug-in of Cytoscape software (version: 3.6.0; Java version: 1.8.0_201) was subsequently used to identify interconnected clusters within the PPI network using a node cutoff score of >3.0. The top 30 proteins with the highest number of degrees (i.e., edges) were represented in the form of a bar graph (); and network modules identified via MCODE (score >3.0) were also represented diagrammatically ().

Animals and Tissues Collection

The animals used in this study were obtained from the Experimental Animal Center of the Medical College of Xi'an Jiaotong University. As per approved animal use protocols 12-week-old female C57B/L mice were euthanized with 5% chloral hydrate, and tissue samples (heart, liver, spleen, lung, kidney, muscles and adipose) were subsequently collected surgically. All surgical instruments used in the experiment were put into 0.1% DEPC solution overnight, and then autoclaved and dried for use. Collected tissue samples were rinsed in pre-chilled Phosphate Buffer Saline (PBS), put into RNase-free centrifuge tubes, and immediately snap frozen in liquid nitrogen. Total RNA was extracted from these tissue samples after transportation to the laboratory.

Total RNA Extraction and cDNA

Frozen tissue samples were homogenized prior to RNA extraction using enzyme-free centrifuge tubes containing Trizol (TakaraBio, Dalian, China), as per manufacturer's instructions. The concentration of the extracted total RNA was determined via nanodrop quantification. Finally, extracted RNA samples were reverse transcribed into cDNA using a Prime Script RT Reagent Kit (TakaraBio, Dalian, China) for subsequent quantitative PCR.

Primer Information

Intron spanning primers were designed using Primer Premier ver. 5.0 (PREMIER Biosoft, http://www.premierbiosoft.com/). The primer sequences, as well as annealing temperatures are described in Table 1.

Table 1

GenesLOCUSSequenceAnnealing temperatureLength
GAPDHNM_001289726.1GTTTCCTCGTCCCGTAG54106
AATCTCCACTTTGCCACT
β-actinNM_007393.5CGTTGACATCCGTAAAGACC58281
AACAGTCCGCCTAGAAGCAC
18s RNANR_003278.3GAAACGGCTACCACATCC56167
ACCAGACTTGCCCTCCA
TNNT3NM_001360157.1TCATTGACAGCCACTTTGAAGC59144
CCAATCTGTTCTGGCGTTCC
PYGMNM_011224.2GAGAGTGGACACGGATGGTG60168
CGGGTATGGGGTCTGATGTG
CMYA5NM_023821.3ACTGGAGCGTGAACAAGGAG60178
GCCATCACCCACACTTGGTA
ENO3NM_001276285.1TCCACGGGTATCTATGAAGCA59133
AGAGCAGGACCTAGAGTCTTGTTG

Primer sequences and annealing temperatures used for quantitative PCR.

qRT-PCR

The CFX-96 (BIO-RAD, US) was used to carry out real-time fluorescence quantitative polymerase chain reaction (qRT-PCR) using a commercially available kit (TB Green Premix Ex Taq II, Tli RNaseH Plus, TakaraBio, Dalian, China). Three reference genes (18s rRNA, GAPDH and β-actin) were tested as internal controls via homogeneity checks. Subsequently, the geometric mean values of 18s rRNA and β-actin was decided to be used as internal reference for qRT-PCR. The final raw data was analyzed via the delta-delta Ct (2−ΔΔCT) calculation method (), and graphical analysis of data was performed via GraphPad Prism 6.0 (https://www.graphpad.com/scientific-software/prism/).

Results

Differentially Expressed Genes Analysis

Microarray data was normalized prior to differentially expressed genes (DEGs) analysis. The gene expression profiles prior to normalization, and after normalization are presented in Figures 1A,B respectively. A volcano plot showing differentially expressed genes (DEGs) that were upregulated or downregulated when contrasting muscle gene expression profiles against a combination of the remaining 17 tissues, is presented in Figure 1C.

Figure 1

Genes that were significantly up or downregulated (P < 0.05) in muscles, with a log2(fold change) ≥2, were identified as differentially expressed genes (DEGs). Gene expression profile of muscles was first contrasted against each of the 17 control tissues used in this study. The number of differentially expressed genes (DEGs) identified in each of these individual contrasts are noted in Table 2. Comparing, gene expression in muscle to the expression profiles of all 17 control tissues combined, allowed for the identification of 260 DEGs that were upregulated, and 12 differentially expressed genes (DEGs) that were downregulated in muscle. Some of the DEGs that were found to be upregulated include myocyte differentiation markers myosin light chain 1 (Myl1), myosin heavy chain 4 (Myh4), myosin heavy chain 2 (Myh2), and inositol protein (Myot). Myosin heavy chain 1 (Myh1) was found to have the highest log2 (fold change) of 6.675, which is indicative of a more than 100-fold greater expression in muscles relative to the combination of the 17 control tissues used in this study. Amongst the genes that were downregulated, Cytochrome C Oxidase Subunit 6A1 (Cox6a1) was found to have the lowest log2 (FC) of−2.472, which is equivalent to an approximately 5-fold reduction in gene expression. A complete list of the top 20 upregulated DEGs, and all of the 12 downregulated DEGs, is presented in Table 3.

Table 2

Control groupTreatment groupUpregulatedDownregulated
Adipose tissueMuscle9871,107
Adrenal glandMuscle9381,079
Bone marrowMuscle1,5221,471
BrainMuscle1,3971,754
EyeMuscle1,1561,316
HeartMuscle519416
KidneyMuscle1,0911,358
LiverMuscle1,3781,316
LungMuscle1,2421,399
Pituitary glandMuscle1,2491,407
PlacentaMuscle1,5061,663
Salivary glandMuscle1,2561,188
Seminal vesicleMuscle1,1851,251
Small intestineMuscle1,4351,572
SpleenMuscle1,4021,588
TestisMuscle2,2402,180
ThymusMuscle1,3871,442
17 tissuesMuscle26012

The number of upregulated and downregulated DEGs identified by contrasting muscle expression profiles against the expression profiles of different tissues.

Table 3

Top 20 Upregulated genesDownregulated genes
Symbollog2FCP-ValueCorrected P-ValueSymbollog2FCP-ValueCorrected P-Value
Tnnt36.67504715.40E-101.66E-07Cox6a1−2.4719250230.0022683880.027740401
Myl16.493218187.14E-087.40E-06Hsp90aa1−2.3181518720.0030972120.035113629
Tnnc26.406027441.91E-082.65E-06Cfl1−2.2658620370.0349294720.18835978
Mylpf6.3565492396.98E-101.96E-07Pgam1−2.2371325290.0398705870.205379883
Myh46.2686546054.99E-101.60E-07Arl6ip1−2.2158692663.95E-050.001211235
Mybpc26.2409637892.08E-111.54E-08Pgrmc1−2.201367293.60E-050.001125671
Myh26.1179351281.97E-108.96E-08Id2−2.1781701830.0004757350.008708711
Actn36.0895122861.14E-105.79E-08Stmn1−2.1473377257.44E-060.000308821
Myot6.0495369382.50E-095.61E-07Spint2−2.1094659737.72E-087.77E-06
Neb5.9709715871.82E-111.47E-08Cks2−2.0290167770.0012170920.017392641
Atp2a15.9257549144.63E-098.87E-07Krt8−2.0180445160.0152057780.107434766
Tnni25.9103193621.74E-108.34E-08Krt18−2.005627755.63E-074.07E-05
Acta15.8342826393.93E-060.000194769
Pvalb5.6710850915.60E-091.01E-06
Myh15.5861125771.43E-093.63E-07
Ckmt25.5376882083.51E-060.000177242
Tcap5.450156195.61E-091.01E-06
Asb55.4436720594.29E-101.49E-07
Pygm5.435759771.21E-081.85E-06
Ckm5.4142145442.96E-083.72E-06

Fold change and statistical significance the top 20 upregulated, all 12 downregulated DEGs when comparing muscle expression profile against the combined expression profile of remaining 17 tissues.

Process and Pathway Enrichment Analysis

Enrichment analysis targeting GO terms identified a total of 752 GO annotations that were significantly enriched (P < 0.01) in differentially expressed genes DEGs identified within this study. Of the total 752 GO terms, 548 GO terms represented biological processes amongst which, the most significantly enriched processes included muscle system process, muscle structure development, myofibril assembly and muscle cell development. A further 103 GO terms representing cellular components were identified, of which the most significantly enriched GO terms including myofibrils, contractile fibers, sarcomeres and contractile fibers. Finally, a total of 101 GO terms representing molecular functions were identified, of which, the most significantly enriched GO terms included cytoskeletal protein binding, actin binding, and structural molecular functions. A summary of the top 10 GO terms identified in each of the three GO aspect categories is presented in Table 4. Top GO terms identified through enrichment analysis are also presented diagrammatically in Figure 2. Complete enrichment analysis results are presented in Supplementary Table 2.

Table 4

OntologyGO IDGO termCountP-ValueFDR
MFGO:0008092Cytoskeletal protein binding712.56E-361.68E-33
GO:0003779Actin binding453.55E-291.16E-26
GO:0008307Structural constituent of muscle131.00E-192.20E-17
GO:0051015Actin filament binding258.69E-191.42E-16
GO:0005515Protein binding1911.36E-171.79E-15
GO:0051371Muscle alpha-actinin binding99.88E-151.08E-12
GO:0042805Actinin binding126.72E-146.30E-12
GO:0051393Alpha-actinin binding103.12E-122.56E-10
GO:0031432Titin binding71.17E-118.55E-10
GO:0005523Tropomyosin binding75.98E-113.92E-09
BPGO:0003012Muscle system process631.40E-504.81E-47
GO:0061061Muscle structure developments711.33E-472.30E-44
GO:0030239Myofibril assembly331.44E-461.66E-43
GO:0055001Muscle cell development451.81E-451.56E-42
GO:0055002Striated muscle cell development433.36E-442.32E-41
GO:0006936Muscle contraction516.66E-443.83E-41
GO:0051146Striated muscle cell differentiation481.24E-396.10E-37
GO:0010927Cellular component assembly involved in morphogenesis331.77E-387.61E-36
GO:0042692Muscle cell differentiation512.04E-387.83E-36
GO:0007517Muscle organ development492.17E-357.49E-33
CCGO:0030016Myofibril784.81E-961.69E-93
GO:0043292Contractile fiber797.03E-961.69E-93
GO:0030017Sarcomere741.06E-921.51E-90
GO:0044449Contractile fiber part751.25E-921.51E-90
GO:0031674I band515.63E-625.42E-60
GO:0030018Z disc432.59E-502.08E-48
GO:0099512Supramolecular fiber831.01E-496.93E-48
GO:0099081Supramolecular polymer831.53E-498.90E-48
GO:0099080Supramolecular complex831.67E-498.90E-48
GO:0005865Striated muscle thin filament211.45E-356.96E-34

GO (Gene Ontology) enrichment analysis of identified DEGs.

Figure 2

Functional enrichment analysis performed using DAVID identified 29 KEGG pathways significantly associated with muscle specific differentially expressed genes DEGs (P < 0.01). The top 20 of these KEGG pathways are presented in the form of a bubble chart in Figure 3, which demonstrates that the identified differentially expressed genes DEGs are highly relevant in cardiac function and pathophysiology. Other KEGG pathways identified via DAVID analysis included AMPK, cGMP-PKG and calcium signaling pathways; Glycolysis / Gluconeogenesis, Carbon metabolism, Arginine and Proline metabolism etc. (Figure 3).

Figure 3

PPI Network Analysis and Module Screening

Network analysis focused on protein-protein interactions performed using the online STRING (https://string-db.org/) database, identified a total of 247 Nodes and 2,813 Edges (score >0.4). Results from the PPI network analysis are presented diagrammatically in Figure 4A, which shows upregulated DEGs in red, and downregulated DEGs in blue, with the color intensity corresponding to fold changes (darker colors reflecting higher fold changes). These results clearly indicate the presence of a large highly correlated network of upregulated muscle specific genes. Network analysis was further performed via Cytoscape software to compute the number of connections of each individual node (i.e., node degrees), and these results are presented in Supplementary Table 4. The top 30 nodes with the highest number of connections (i.e., degrees), presented in Figure 4B, were comprised by Titin (Ttn, 103 degrees); Actinin Alpha 2 (Actn2, 86 degrees); Creatine Kinase, Mitochondrial 2 (Ckmt2, 83 degrees); LIM Domain Binding 3 (Ldb3, 83 degrees); Muscle Creatine Kinase (Ckm, 81 degrees); Obscurin, Cytoskeletal Calmodulin and Titin-Interacting RhoGEF (Obscn, 80 degrees); and Titin-Cap (Tcap, 80 degrees), in addition to many other muscle-specific DEGs.

Figure 4

Network module analysis performed via MCODE plug-in of Cytoscape, identified a further five modules (module score > 3.0), which are presented in Figures 4C–G. Module 1 (Figure 4C) has the highest MCODE score of 36.905, and included 43 interacting proteins (nodes) with 775 interactions (edges). The second module (Figure 4D) was considerably smaller with an MCODE score of 9.769, including 27 nodes and 127 edges. As evident in Figures 4C–G, most of the DEGs represented in these modules were upregulated (Red nodes indicate upregulated nodes, and blue indicates downregulated nodes).

qRT-PCR Validation of Identified DEGs

Differentially expressed genes DEGs that were identified in this study were also annotated for tissue specific expression using the online DAVID (https://david.ncifcrf.gov/) database, and this identified 55 genes with tissue specific expression in skeletal muscles (Table 5). A Venn diagram constructed to compare these 55 genes against the list of top 30 genes identified via Cytoscape network analysis, identified 11 genes shared in common (Figure 5). These genes included: Actinin Alpha 2 (ACTN2), LIM Domain-Binding Protein 3 (LDB3), Small Muscle Protein X-Linked (SMPX), Caveolin 3 (CAV3), Troponin T3, Fast Skeletal Type (TNNT3), Myozenin 2 (MYOZ2), Glycogen Phosphorylase, Muscle Associated (PYGM), Cardiomyopathy Associated 5 (CMYA5), Enolase 3 (Beta, Muscle) (ENO3), Sarcalumenin (SRL) and Actinin Alpha 3 (ACTN3).

Table 5

TissuesCountP-ValueGenes
Skeletal muscle552.15E-57PRKAG3, PDLIM3, ANKRD2, RTN2, ART3, JSRP1, PVALB, SH3BGR, JPH1, RBFOX1, MYH1, MYH2, TNMD, LDB3, MYH4, ACTN2, MYH7, FBP2, ACTN3, CACNG1, TACC2, TNNT3, TNNT1, SGCG, CFL2, ITGB1BP2, RYR1, SMPX, SGCA, SGCB, CAV3, FHL1, PHKA1, SRL, TPM2, KCNA7, ENO3, DDIT4L, HRC, MYF6, MUSTN1, MYOZ2, TRIM63, TNNI1, SLC16A3, TUBA8, NEB, NRAP, MAPK12, PYGM, GYG, ATP2A1, CMYA5, SYNM, VLDLR
Heart588.55E-19APOBEC2, LDHA, TNNC1, PGAM2, ANKRD1, TXLNB, TTN, ART1, LMOD2, PPP1R14C, USP13, HSP90AA1, SLC25A4, CRYAB, LDB3, ACTN2, MYH7, GMPR, IRS1, TRDN, MURC, HSPB6, CFL2, HSPB8, HSPB7, SMPX, TCAP, MYL3, SMTNL2, ASB14, TPM1, ASB15, KCNA7, MYOM2, CKMT2, MLIP, POPDC3, HRC, LPL, ADSSL1, ACTC1, CACNA2D1, COX8B, ALPK3, YIPF7, PDK4, ATP1A2, TRIM63, CSRP3, CACNA1S, IDH3A, FSD2, ABCC9, PTP4A3, FABP3, COX6A2, KLHL30, VLDLR
Muscle123.50E-11RBFOX1, ADSSL1, SLC25A4, SLC2A4, PHKG1, CMYA5, HSPB7, KY, SYNPO2, MYOM1, MYOT, SNTA1
Heart muscle52.21E-06ART3, TRIM54, NRAP, XIRP1, MYL2
Bone151.54E-05PRKAG3, FSD2, SLC2A4, PHKB, MYPN, PHKA1, SRL, MYLK2, TMOD4, ATP1A2, TXLNB, STAC3, ASB15, NMRK2, MB

Tissue specific expression annotations of DEGs identified via DAVID analysis.

Figure 5

qRT-PCR was performed to determine mRNA expression levels in seven different murine tissues including the heart, liver, spleen, lung, kidney, muscle and adipose tissue to validate muscle specific expression of selected genes. Results from qRT-PCR (Figure 6), confirmed high levels of expression of TNNT3, PYGM, ENO3, CMYA5 in muscles, reaffirming the validity of the findings in this study. We used 18S rRNA, β-actine and GAPDH as housekeeping genes for the mRNA expression analysis of DEGs in the target tissues. Although GAPDH is not considered a very suitable option for using as a reference gene (), however, we used triple reference genes for the expression of mRNA levels in all target tissues.

Figure 6

Discussion

The overall aim of the current study was to identify genes specifically expressed in skeletal muscles via bioinformatic analyses of publicly available microarray data, followed by qRT-PCR to validate muscle specific expression of selected genes in an independent set of samples. Bioinformatic analysis of publicly available microarray data resulted in the identification of 272 DEGs with at least a 4-fold expression level relative to expression in 17 other mice tissues. The majority of these genes were upregulated (n = 260), and a very small proportion of genes were downregulated (n = 12). These results suggest that upregulation of key genes is more crucial for muscle physiology and development, relative to downregulation of specific genes. A previous study () which describes gene expression profiles of different tissues (kidney, liver, lung, heart, muscle, and adipose tissue), also reported that key genes (e.g., Myot, Tnnc2, Tnni2, Tnnt3, Actn3, Mybpc1, Mybpc2, Myoz1) are highly upregulated in both human and murine muscles. In our study, we have also found almost all of the above genes to also be highly upregulated (Table 3) in muscles, and therefore our results align with these previous findings. A number of these DEGs have also previously been reported to be involved in muscle development. Some of these genes include Tnnt3 (, ), Myh1, Myh2, Myh4 () and Actn3 (), Pvalb () Ckmt2 () Cox8b (). However, several other genes that were identified to be differentially expressed in muscles have not yet been reported to have a role in muscle development or physiology (, ).

Enrichment and pathway analyses performed identified several GO annotations terms (n = 752) to be significantly enriched in the 272 genes identified to be differentially expressed in skeletal muscles. The GO enrichment demonstrated significant involvement of ontologies relevant to muscle development, physiology and function; which in turn accords with findings from DEG analyses. Pathway analysis also identified 29 KEGG pathways, several of which were relevant to muscle development. However, one of the more interesting findings here was that the top four pathways identified were all associated with cardiac muscle physiology and pathology. This could suggest that some genes specifically expressed in skeletal muscles, could be involved in cardiac myopathies.

While the skeletal muscle samples used in this study were derived from 10–12 week old mice, these samples include both satellite cells and skeletal muscle-derived stem cells, which together form the pool of cells required for myogenesis (31). When satellite cells are activated (e.g., due to muscle injury), they get induced to undergo myogenic differentiation, which in turn requires highly specific temporal and spatial expression patters of different transcription factors and proteins (32) that is consistent with findings in this study. Similarly, skeletal muscle stem cell proliferation and muscle differentiation can also be triggered in adults under the influence of hormones like IGF1 (33), which in turn activates a number of downstream pathways including MAPK, PI3K-AKt-mTOR-P70S60K and PI3K-AKt-mTOR-GSKβ signaling pathways (3437). Therefore, the identification of several genes, ontologies and pathways associated with muscle development is not surprising.

To affirm the findings from differentially expressed genes DEGs, enrichment and pathway analysis, we constructed a PPI protein interaction network map, consisting of 2,813 edges (interactions) between 247 nodes (proteins). The PPI network map identified several structural proteins and enzymes as core nodules (e.g., TMOD4, MYL1, MYBPC2, ATP2A1). When degree scores were computed via Cytoscape network analysis, the top nodes were also mainly comprised of structural genes, and genes involved in muscle physiology and function (e.g., TNN, ACTN2, LDB3, CKMT2). Network module analysis via the MCODE plug-in also identified 7 interaction modules. The largest of these modules was comprised of a total of 43 nodes and 775 edges, and included several structural myosin-related (Myh7, Myl2, Myl3, Myh2, Myh4, Myh1) and actin-related (Acta1, Actc1, Actn2, Actn3) proteins. Overall, findings from enrichment, pathway and network analysis were in accord and reaffirmed the involvement of identified DEGs in muscle structure and physiology.

Module analysis of the PPI network identified several genes that have been previously reported to be involved in muscle development (e.g., Module 1, Figure 4C). However, several interesting candidates, whose roles in muscle development are yet to be characterized, were also identified. Examples of such genes include Tripartite motif-containing 54 (Trim54), Creatine kinase, mitochondrial 2 (Ckmt2), cardiac disease associated 5 (Cmya5) and Leiomodin 2 (Lmod2). Future research aimed at characterizing the function of these genes could offer novel insights into mechanistic aspects of muscle development and associated pathophysiology.

Finally, we used qRT-PCR to validate the expression patterns of selected genes that were identified as specifically expressed in skeletal muscles (via DAVID analysis), and were also identified within the top 30 genes of the PPI network (i.e., those having the highest scores). The obtained results are consistent with the results from microarray DEG analyses, which reaffirms the findings from bioinformatic analyses of the microarray data.

Conclusion

In conclusion, 272 genes with muscle-specific expression profiles were identified in this study, which included several genes widely known to be involved in muscle development and function. Downstream enrichment and pathway analysis identified several muscle specific ontologies and pathways reaffirming findings of differentially expressed genes DEG analysis. Validation of results in an independent set of samples via qRT-PCR also reaffirmed muscle specific expression of selected DEGs. Several of the 272 differentially expressed genes DEGs identified in this study are yet to be functionally characterized in context of muscle development and physiology. Once characterized, these candidate genes could offer new targets for development of mutant mouse models of human muscle associated diseases and disorders. Therefore, future research aimed at investigating the role of these candidate genes in the context of muscle development and physiology is warranted.

Funding

The present study was supported by grants from the National Key Research and Development Program of China (2018YFD0501700), National Natural Science Foundation of China (31972994), Key Research and Development Program of Ningxia Province (2019BEF02004), and National Beef and Yak Industrial Technology System (CARS-37).

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Statements

Data availability statement

All datasets used in this article are public and sources cited accordingly. The screening data respectively, source is available at: https://david.ncifcrf.gov; https://david.ncifcrf.gov/home.jsp; http://kobas.cbi.pku.edu.cn/; http://ggplot2.tidyverse.org; https://www.r403project.org/; http://www.bioconductor.org/packages/3.5/bioc/html/limma.html.

Ethics statement

All procedures were approved by the Experimental Animal Center of Xi'an Jiaotong University. Animal care and use protocols (EACXU 172) were approved by the Institutional Animal Care and Use Committee of Xi'an Jiaotong University and Northwest A&F University, Yangling. All animal experiments were performed in adherence with the NIH Guidelines on the Use of Laboratory Animals.

Author contributions

SR conceived, conceptualization, and designed the experiments. SR and WG performed the experiments and data analysis. CL and WG contributed to data curation. GC and CM contributed to methodology. ZM, AS, and MA contributed to the investigation, methodology, and validation. SP and RK provided constructive suggestions for the discussion and validation and contributed in drafting, editing, and review of manuscript. NS revised manuscript critically for content and grammar. LZ contributed to project administration, supervision of the overall study, and provided the necessary resources. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

All authors acknowledge and thank their respective institutes and universities.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2021.794628/full#supplementary-material

Supplementary Table 1

Screening of DEGs in mouse muscle tissue.

Supplementary Table 2

Details of GO analysis.

Supplementary Table 3

Details of KEGG pathway analysis.

Supplementary Table 4

Analysis of the degree of PPI network of DEGs.

References

Summary

Keywords

muscle development, microarray analysis, differential genes, bioinformatics, biotechnology

Citation

Raza SHA, Liang C, Guohua W, Pant SD, Mohammedsaleh ZM, Shater AF, Alotaibi MA, Khan R, Schreurs N, Cheng G, Mei C and Zan L (2021) Screening and Identification of Muscle-Specific Candidate Genes via Mouse Microarray Data Analysis. Front. Vet. Sci. 8:794628. doi: 10.3389/fvets.2021.794628

Received

13 October 2021

Accepted

22 November 2021

Published

13 December 2021

Volume

8 - 2021

Edited by

Adriana Mércia Guaratini Ibelli, EMBRAPA Swine and Poultry, Brazil

Reviewed by

Ali Raza Jahejo, Shanxi Agricultural University, China; Paolo Zambonelli, University of Bologna, Italy

Updates

Copyright

*Correspondence: Linsen Zan

This article was submitted to Livestock Genomics, a section of the journal Frontiers in Veterinary Science

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics