Identification of DEMs in denervated muscle
Based on the analysis of the GSE81914 dataset, the DEMs list was obtained with the | logFC | ≥1 and p-value < 0.05 as the threshold. Compared with innervated muscle samples, a total of 21 DEMs were detected in denervated muscle samples, including 16 upregulated DEMs and 5 downregulated DEMs. As shown in the volcano plot and cluster heatmap, we extracted a total of 21 DEMs from 1242 miRNAs in the GEO expression matrix and obtained the distribution of upregulated and downregulated DEMs (Figure 1).
Prediction and validation of upstream transcription factors of DEMs
TransmiR v2.0 database was used to predict the potential upstream transcription factors of candidate DEMs. Of the 21 obtained DEMs submitted to TransmiR v2.0, 11 DEMs were predicted for 14 transcription factors (Figure 2a). One transcription factor can act on multiple miRNAs, and one miRNA can also be regulated by more than one transcription factor. For the downregulated DEM, mir-495, the transcription factors were Mlxipl, while for the upregulated DEMs, the transcription factors were Brd4, Hnf4a, Nfkb1, Rela, Gh1, Myod1, Camta1, Mecp2, Nkx2-2, Sox10, Ep300, and Nr3c1. Hnf4a acted on 6 DEMs such as mir-21, mir-132, mir-203, mir-326, mir-511, and mir-1949 simultaneously, while mir-132 was be regulated by 5 transcription factors (Nkx2-2, Mlxipl, Mecp2, Hnf4a, Camta1).
As shown in TF-miRNA regulatory network, Hnf4a, Nkx2-2, Mlxipl, and Med1 can simultaneously modulate more than one DEM, so we consider that they play a more important role in regulating denervated muscle atrophy. The qRT-PCR was applied to verify the expression trend of transcription factors in the two groups of muscle tissues. The results showed that compared with the sham group, Med1 was significantly up-regulated in the denervated muscles, while Hnf4a, Nkx2-2, and Mlxipl had no significant statistical difference between the two groups (Figure 2b).
Prediction and construction miRNA-mRNA regulatory network
Among the candidate DEMs submitted to the miRNet database, a total of 59 target genes were predicted by 11 DEMs, of which 10 up-regulated DEMs predicted 55 target genes, and 1 down-regulated DEM predicted 5 target genes. For better visualization, we constructed a miRNA-mRNA regulatory network of DEMs and their target genes through Cytoscape software (Figure 3). In addition, the target genes predicted by each DEM are listed in Table 1.
GO and KEGG analysis of target genes
To clarify the biological characteristics of target genes related to denervation-induced muscle atrophy, the DAVID was used to perform (Gene Ontology) GO analysis. The obtained results of the enriched GO terms were shown in figure 4. GO functional enrichment results include gene ontology biological process (BP), molecular function (MF), and cellular components (CC). Our results show that target genes mainly involve biological processes such as cellular response to amino acid stimulus, endodermal cell differentiation, positive regulation of apoptotic process, cell adhesion, cell-matrix adhesion, and cell migration. In the MF group, target genes were mainly enriched in extracellular matrix structural constituent, protein binding, fibronectin binding, protein kinase binding, identical protein binding, kinase binding, and glycoprotein binding. As shown in figure 4, target genes mainly involve the following cellular components, such as extracellular matrix, collagen trimer, basement membrane, cytosol, collagen type V trimer, cytoplasm, and extracellular space. These data indicated that in the denervation-induced muscle atrophy model, cells in the extracellular matrix, cytoplasm, and collagen trimer perform molecular functions such as protein binding, fibronectin binding, protein kinase binding, and extracellular matrix structural constituent, resulting in cell apoptosis, cell adhesion, and cell migration.
The (Kyoto Encyclopedia of Genes and Genomes) KEGG analysis was applied to explore enriched pathways of target genes. The enriched categories of KEGG pathways are shown in figure 4d. Results showed that the PI3K-Akt signaling pathway, Protein digestion and absorption, Focal adhesion, ECM-receptor interaction, and FoxO signaling pathway were enriched in denervated muscles, which provided molecular pathway mechanisms for studying denervated muscle atrophy.
Module analysis and hub genes expression validation
The target genes of denervated gastrocnemius muscle were constructed to construct protein-protein interaction (PPI) networks by using the STRING database. A complex PPI network with 49 nodes and 175 edges was constructed after removing the isolated nodes. Two significant modules with high corresponding degrees, all of which have a score ≥ 5, were screened by using MCODE in Cytoscape. GO and KEGG Enrichment analysis of module 1 and module 2 were displayed in figure 5. As shown in figure 5, the target genes of module 1 are regulated by mir-29b and mainly participate in biological processes such as cellular response to amino acid stimulus, endodermal cell differentiation, collagen fibril organization, collagen-activated tyrosine kinase receptor signaling pathway, and cell adhesion. The target genes module 2 are mainly enriched in the biological process of cell apoptosis. The most significant pathway in module 1 and module 2 was simultaneously enriched in the PI3K-Akt signaling pathway.
Among the 49 nodes, there are 18 central nodes were filtered with the criterion of degree score≥10. The 18 most crucial genes that shared high interaction were Col1a1, Mmp9, Col4a2, Col3a1, Pten, Col4a1, Mmp2, Col18a1, Col5a2, Vegfa, Ctgf, Foxo3, Col5a1, Col8a1, Col16a1, Col7a1, Col5a3 and Col12a1. Among them, 12 encode collagen-related genes, the target genes of mir-29b, are completely located in module 1(Figure 5a). Therefore, we consider that the collagen-related genes regulated by mir-29b play an important role in muscle atrophy. Since COL1A1 has the highest degree value among hub genes and is predominantly expressed in muscle, we selected COL1A1 to represent collagen-related genes for experimental verification. Then, we used qRT-PCR to verify the expression trend of COL1A1, Mmp9, Pten, Mmp2, Vegfa, Ctgf, Foxo3(Figure 6). For the up-regulated DEM (mir-29b, mir-132, mir-212, and mir-203), unlike the decreased expression of COL1A1, there was no significant change in the expression of Pten, while the expression of Mmp9, Mmp2, FOXO3, and Vegfa increased in contrast. For the down-regulated DEM (mir-133a), the expression of Ctgf was significantly increased. Thus, mir-29b-COL1A1 and mir-133a-Ctgf were identified as potential regulatory pathways in denervated muscle.