Integrative single-cell RNA-seq and ATAC-seq analysis of myogenic differentiation in pig

scRNA-seq identified the major cell types in developing pig somites

To gain a comprehensive view of the cell populations present during pig skeletal muscle ontogeny, we used scRNA-seq to evaluate somite and myotome tissues of ZZ and DZ from embryonic day (E) 16 to E28, which covers the transition from Pax3+ progenitors to myocytes. Single cells from 8 samples (E16-ZZ, E16-DZ, E18-ZZ, E18-DZ, E21-ZZ, E21-DZ, E28-ZZ, and E28-DZ) were processed for scRNA-seq using a Chromium system (10× Genomics) (Fig. 1A). Overall, 70,201 cells passed quality control (QC) with an average of 1892 genes per cell and 5276 unique molecular identifiers (UMIs) per cell (Additional file 1: Figure S1A, B). To exclude technical batch effects, the datasets from all samples and tissues were merged using autoencoders (AEs) and applied the batch-balanced k nearest neighbors (BBKNN) approach [26] to the latent space [27] (Additional file 1: Figure S1C). Dimensional reduction and unsupervised clustering for all 70,201 cells identified 31 distinct clusters (Additional file 1: Figure S1D). Based on differential expression analysis and the expression of selected marker genes from the literature, we manually annotated 12 distinct populations (Additional file 2: Table S1 and Fig. 1B, C), including mesenchymal cells, fibroblasts, epithelial cells, neural stem cells, myogenic progenitors/myoblasts, osteogenic cells, neurons, neurogliocytes, endothelial cells, myocytes, chondrocytes, and muscle cells. Bubble plots of marker gene expression demonstrated the accuracy of the cell annotations (Fig. 1D). Gene Ontology (GO) analysis of the differentially expressed genes (DEGs) for each cell type verified the characteristics and functions of different cells (Additional file 3: Table S2 and Additional file 1: Figure S2). We next quantified changes in the cell type percentages during somite development. As shown in Fig. 1E, the patterns of the cell populations vary considerably across development stages, with the proportions of low differentiation cells, such as epithelial cells, endothelial cells, and neural stem cells, decreased with development, while the percentages of highly differentiated cells, such as fibroblasts, osteogenic cells, and myocytes, increased gradually with development, suggesting a rapidly somite development during E16 and E28.

Fig. 1figure 1

scRNA-seq identified major cell types in developing pig somites. A Experimental workflow schematic. Somite tissues [a mixture of embryos (n=5) from the same sow] of ZZ and DZ at E16, E18, E21, and myotome tissues [a mixture of embryos (n=5) from the same sow] at E28 were isolated. Tissue samples were dissociated into a single-cell solution and then single-cell transcriptomes were captured and analyzed using 10×Genomics. The minimum scale of the ruler in the embryo photograph is 1 mm. B t-Stochastic neighbor embedding (tSNE) plots showing the distribution of the main cell populations using the scRNA-seq. Using marker genes, cells were annotated as mesenchymal cells, fibroblasts, epithelial cells, neural stem cells, myogenic progenitors/myoblasts, osteogenic cells, neurons, neurogliocytes, endothelial cells, myocytes, chondrocytes, or muscle cells. Colors indicate cell types. Each dot represents one cell. C Heatmap showing the top 20 markers for each of the 12 cell populations. D Dot plot of the mean expression of canonical marker genes for 12 cell populations. E Bar plot showing the percentage of different cell types within each sample. F Visualization of myogenic cell (including myogenic progenitors/myoblasts, myocytes, and muscle cells in Fig. 1B) sub-clusters via t-SNE by developmental stage (left) and sub-cluster number (right). G tSNE plots showing the cell identities of myogenic cell sub-clusters. H Violin plots showing feature gene expression in each cell sub-cluster. Colors represent sub-clusters described in Fig. 1G

To further dissect the cellular heterogeneity and transcriptional landscape of developing myogenic cells, we extracted myogenic progenitors/myoblasts, myocytes, and muscle cells for further clustering. The myogenic cells were further divided into 8 distinct sub-clusters with increased resolution and annotated as Pax3+ progenitors, myogenic progenitors, myoblasts, myocytes, cardiac muscle cells, and other cells (Fig. 1F, G). Among them, Pax3+ progenitors were characterized by the highest expression level of Pax3, whereas myogenic progenitors and myoblasts were characterized by the expression of the muscle stem cell marker Pax7 and the myogenic regulatory factor MyoD. Both myocytes and cardiomyocytes expressed the muscle cell marker ACTC1, but the differentially expressed skeletal myogenic cell markers MYOG and MYL1 as well as the cardiomyocyte-specific marker gene MYBPC3 [28] accurately distinguished the two cell types (Fig. 1H).

Reconstruction of the myogenic differentiation trajectory of Pax3+ progenitors

To investigate the molecular processes underlying skeletal muscle development, the cells were ordered in a pseudotime manner using Monocle 2 [29]. Pseudotime trajectory analysis revealed seven different cell states (states 1~7) and presented the distributions of cell states along with pseudotime flows (Fig. 2A). An organized, branched progression of cells from Pax3+ progenitors to differentiated myocytes was shown by labeling individual cells using the cell population annotations from the unified atlas in Fig. 1E (Fig. 2B, C). Unexpectedly, the first small branch (state 2) was almost entirely enriched by myocytes, and the cells at the end of pseudotime trajectory (state 7) also belonged to myocytes (Fig. 2A, C). The differential gene expression analysis indicated that many muscle development-related genes were highly expressed in both states of myocytes (e.g., MYMK, FNDC5, MEF2C, and TNNI1). However, the myocytes in state 2 expressed the cardiomyocyte-specific marker MYL9, while those in state 7 expressed the skeletal muscle cell-specific markers MYOD1 and MYOG (Fig. 2D). In addition, the DEGs with high levels in state 2 were largely involved in the regulation of biological processes such as “heart process” and “cardiac cell development,” while the DEGs with highly levels in state 7 were involved in “skeletal muscle tissue development” and “skeletal muscle cell differentiation” (Fig. 2E). These results indicated that the myocytes in state 2 were actually cardiomyocytes which were excluded from subsequent analysis (Fig. 2F).

Fig. 2figure 2

Reconstruction of the myogenic differentiation trajectory in a pseudotime manner. A Pseudotime analysis of myogenic cells (including Pax3+ progenitors, myogenic progenitors, myoblasts, and myocytes in Fig. 1E) was performed by Monocle 2 and revealed seven different cell states (states 1~7). The distributions of cell states were presented along with pseudotime flows. Each dot is a cell. B Visualization of myogenic differentiation trajectory by cell origins (left) and developmental stages (right). C Visualization of myogenic differentiation trajectory by cell identity. D Violin plots showing feature gene expression in each cell cluster. MYMK, FNDC5, MEF2C, and TNNI1 are muscle development-related genes. MYL9 is a cardiomyocyte-specific marker. MYOD1 and MYOG are skeletal muscle cell-specific markers. E Gene Ontology (GO) analysis of the differentially expressed genes with high levels in myocytes (state 2) and myocytes (state 7). F Visualization of myogenic differentiation trajectory by cell state, with cardiac cells distinguished from differentiating skeletal muscle cells. G Bar plot showing the percentage of cells within each sample assigned to the annotated myogenic cell types. H Immunofluorescence staining for Pax7 and MyoD on somite cross sections of ZZ and DZ at E21 and E28. Scale bar = 100μm

Then, the percentages of Pax3+ progenitors, myogenic progenitors, myoblasts, and myocytes in each sample were calculated to assess myogenesis progression at different stages. Although the proportion of total myogenic cells in all cells of somite tissue did not change obviously (Fig. 1E), the percentages of four types of myogenic cells in different developmental stages changed significantly during E16~E28. At E16, there were almost no myoblasts and myocytes in the somites, with pax3+ progenitors accounting for about 90% of the four types of cells. At E18, some progenitors were committed to become myoblasts and further differentiated into myocytes. Subsequently, the proportion of pax3+ progenitors decreased markedly, and myogenic progenitors, myoblasts, and myocytes accounted for about 90% of the four types of cells at E21 and E28 (Fig. 2G). Dual immunostaining of Pax7 and MyoD showed that MyoD+ cells appeared at E21 and increased significantly at E28, indicating that E18~28 was a critical period for the establishment of skeletal muscle lineage (Fig. 2H). As expected, cells from E16 and E18 embryos tended to be distributed at the root of the trajectory suggesting versatile progenitor properties, whereas those from E21 and E28 embryos were distributed in the later part of the trajectory indicating the decreased proportion of progenitor cells and increased proportion of differentiating myogenic cells during skeletal muscle development (Fig. 2B).

In addition to differentiation, proliferation is also a major biological event during early skeletal muscle development. The cell cycle phase of each cell was predicted to evaluate the proportion of proliferating cells at different states of myogenic progression. We observed a shift in the transcriptomically defined cell cycle state accompanying the change in cell type representation (Fig. 2C and Additional file 1: Figure S3A). In the early part (states 1, 3, and 4) of the pseudotime trajectory, most of the cells are progenitors, so they are primarily predicted to be proliferating (S and G2M phases), with only a small fraction being predicted to be in G1 phase (non-proliferating). However, in the middle stage of the trajectory (State 5), proliferating cells decreased to 73.9%, while non-proliferating cells increased significantly to 26.1%. At the end of the trajectory (States 6 and 7), most of the cells have differentiated into myocytes that no longer proliferate, so they are primarily in G1 phase (75.3–85.1%) (Additional file 1: Figure S3B). Consistently, the variation in the expression of cell cycle-related genes and myogenic genes suggested that proliferating cells decreased with the differentiation trajectory (Additional file 1: Figure S3C, D).

Transcriptome dynamics of Pax3+ progenitor differentiation

To gain insights into the gene expression dynamics along the trajectory, the expression changes of the 1700 top DEGs among the four types of myogenic cells (including Pax3+ progenitors, myogenic progenitors, myoblasts, and myocytes) were analyzed and clustered into five major categories of transcriptional gene clusters (Additional file 4: Table S3 and Fig. 3A). Genes in cluster 1, with the highest expression in Pax3+ progenitors, were gradually downregulated from the beginning of programming and were largely involved in the regulation of biological processes such as “ATP metabolic process” (e.g., ACADS, ATIC, and MEIS1) (Fig. 3B). Subsequently, gene clusters 2 and 3 were transiently upregulated but finally downregulated, representing two temporary transcriptional waves. Cluster 2 genes highly expressed in myogenic progenitors were involved in “mitotic cell cycle process” (e.g., PCLAF, CENPA, and HMGB2) indicating the strongest proliferation capacity of myogenic progenitors during myogenesis. Cluster 3 genes highly expressed in myogenic progenitors and myoblasts are involved in “extracellular matrix” and “response to growth factor” (e.g., HES1, HEY1, and SIX1) (Fig. 3B). Concurrent with cluster 3 activation, cluster 4 genes highly expressed in myoblasts and myocytes were upregulated and maintained at high expression levels until the final stage with enrichment of the GO terms “response to endogenous stimulus” and “striated muscle tissue development” (e.g., TCF12, FOS, and EGR1). Finally, cluster 5 genes were activated at the end of the trajectory with predominant involvement in the GO term “muscle cell differentiation” (e.g., KLF2, RHOB, and SOX6) (Fig. 3B). These data illustrated the trajectory of myogenic differentiation and revealed the ordered activation of transcriptional waves throughout this process.

Fig. 3figure 3

Transcriptome dynamics of the myogenic differentiation. A Heatmap showing the expression changes of the 1700 top differentially expressed genes (DEGs) in a pseudotemporal order, with the DEGs, were cataloged into 5 five major clusters in characterized patterns (right). The GO analysis was performed for each gene cluster, and the representative enriched biological process (BP) terms are presented (left). B The expression dynamics of DEGs in gene clusters. Thick lines indicate the average gene expression patterns in each cluster (left). Gene signatures and expression dynamics of representative genes in each gene cluster (right)

Single-cell chromatin accessibility profiling of pig skeletal muscle ontogeny

To further investigate the regulatory events in developing myogenic cells, the single-cell chromatin accessibility landscape was analyzed (Additional file 1: Figure S4A-C). Using a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm, we obtained 15 distinct clusters of differentially accessible peaks (Additional file 5: Table S4 and Additional file 1: Figure S4D). Clusters 4 and 8 were annotated as myogenic cells for their high accessibility of marker genes associated with myogenic lineage (Figure. S4E, F). Then, the myogenic cells were further divided into 7 distinct sub-clusters with increased resolution (Fig. 4A).

Fig. 4figure 4

Single-cell chromatin accessibility analysis of pig myogenic cells. A The myogenic cells in the scATAC-seq dataset are shown in the Uniform Manifold Approximation and Projection (UMAP) space, colored by cluster. B Top: Bar plot showing the average accessibility of 13 selected marker genes from our scRNA-seq data considering all myogenic cells. Bottom: dot plot of the standardized accessibility of the marker genes (gene body ± 2 kb) in each of the seven clusters. For each gene, the minimum accessibility value is subtracted, and the result is divided by its maximum accessibility value. The dot size indicates the percentage of cells in each cluster in which the gene of interest is accessible. The standardized accessibility level is indicated by color intensity. C UMAP visualization of the myogenic cells in the scATAC-seq dataset, colored by cell identity. D Percentage distribution of open chromatin elements in each scATAC-seq sample. E Percentage distribution of open chromatin elements in scATAC-seq myogenic cell types. F Heatmap showing cell type-specific differentially accessible peaks (DAPs) (yellow: open chromatin, purple: closed chromatin). G Distribution of open chromatin elements among DAPs in myogenic cell types. H Number of shared and unique peaks among snATAC-seq cell types

To explore the chromatin accessibility profiles across the seven clusters, we examined the accessibility of selected marker genes from our scRNA-seq data (Fig. 4B). In clusters 2 and 4, we observed greater accessibility of marker genes associated with myocytes (e.g., MYOG, MYH3, MYL1, and CKM) and lower accessibility of genes associated with progenitor cells (e.g., PAX3 and PAX7) (Additional file 6: Table S5 and Fig. 4B). In contrast, cells in clusters 0, 3, and 5 showed greater accessibility for marker genes of progenitor cells and lower accessibility for marker genes of myocytes. Clusters 1 and 6 had mixed signatures, with greater accessibility of PAX7, MSC, MYF5, and MYOG. Based on these observations, we manually annotated the seven clusters as Pax3+ progenitors, myogenic progenitors/myoblasts, and myocytes (Fig. 4C).

To characterize different genomic elements captured by scATAC-seq data, the genome was stratified into promoters, exons, 5′ and 3′ untranslated regions, introns, and distal regions using the GENCODE annotation [30] (Additional file 1: Figure S5A, B). There was little difference in the proportions of different regions between samples or cell types, with exons accounting for about 50%, promoters, introns, and distal regions accounting for about 15% each, and 5′ and 3′ untranslated regions accounting for about 5% (Fig. 4D, E). To study the open chromatin heterogeneity across cell types and developmental stages, we derived a cell type-specific chromatin accessibility landscape by conducting pairwise Fisher’s exact test for each peak between every cluster. In total, we identified 6422 differentially accessible open chromatin peaks (DAPs) across the 3 cell types, which separated the three cell types perfectly (Additional file 7: Table S6 and Fig. 4F). Among these peaks, most were in regions characterized as distal elements or introns, while relatively few (<10%) were in the promoters or 5′ and 3′ untranslated regions (Fig. 4G and Additional file 1: Figure S5C), indicating a critical role for enhancer elements in skeletal muscle development. In addition to the cell type-specific peaks, some cell type-independent open chromatin areas (present across Pax3+ progenitors, myogenic progenitors/myoblasts, and myocytes) also were found, likely consisting of basal housekeeping genes and regulatory elements (Fig. 4H). The overlapping peaks between Pax3+ progenitors and myogenic progenitors/myoblasts, and between myogenic progenitors/myoblasts and myocytes, were more than that between Pax3+ progenitors and Myocytes, which is consistent with their biological similarities and differentiation process (Fig. 4H). The common peaks of the three cell types were much more than those of the other groups, revealing their close lineage relationship.

Cell type-specific gene regulatory landscape of embryonic skeletal muscle in pigs

Cell type-specific chromatin opening and closing events associated with TF binding changes establish the cell type-specific regulatory landscape, resulting in cell-type specification and development. Therefore, the motif enrichment analysis was performed on the cell type-specific open chromatin regions using 10× Genomics. The full list of cell type-specific TF binding motifs is shown in Additional file 8: Table S7. We next correlated the motif enrichment with scRNA-seq TF expression (Additional file 9: Table S8). Using this combined motif enrichment and gene expression approach, the pig skeletal muscle cell type-specific TF landscape was defined (Fig. 5A, B). Correlation of RNA expression and chromatin accessibility in individual single cells revealed two characteristic patterns of ATAC–RNA pairs: (i) RNA expression of TFs directly matches accessibility of corresponding TF bindings sites as exemplified for ARID3A, MEIS2, MEIS1, HOXB4, and HOXD4 in Pax3+ progenitors, and MYOG, MYOD1, KLF2, and SOX6 in myocytes, suggesting that these TFs actively regulate their respective target genes at the specific developmental stage; (ii) RNA expression of TFs precedes the increase in accessibility of the corresponding TF binding sites. This scenario was apparent for MYF6, MYF5, SNAI1, and TGIF1, which reached their highest expression levels in myogenic progenitors/myoblasts but showed the strongest motif enrichment in myocytes, suggesting that additional epigenetic regulation could occur before TFs take action (Fig. 5A, B).

Fig. 5figure 5

Cell type-specific gene regulatory landscape of pig embryonic skeletal muscle. A Heatmap showing motif enrichment analysis on the cell type-specific open chromatin regions using 10× Genomics (full results are shown in Additional file 7: Table S6). B TF expression z score heatmap that corresponding to the motif enrichment in each cell type. C Heatmap of cell type-specific regulons, as inferred by the SCENIC algorithm. Regulon activity was binarized to “on” (black) or “off” (white). D tSNE depiction of regulon activity (“on-blue”, “off-gray”), TF gene expression (red scale), and expression of predicted target genes (purple scale) of exemplary regulons for Pax3+ progenitors (MEIS1), myogenic progenitors (EZH2 and HDAC2), myoblasts (EGR1), and myocytes (MYOG). Examples of target gene expression of the TFs (PAX3, PCNA, SRSF7, RHOB, and SPG21) are shown in purple scale. Additional examples are given in Figure S6. The full list of regulons and their respective predicted target genes can be found in Additional file 11: Table S10

To study the putative target genes of TFs, single-cell regulatory network inference and clustering (SCENIC) was performed to examine TF regulon activity [31]. The activity of each regulon in each cell was quantified and then binarized to “on” or “off” based on activity distribution across cells. The SCENIC results indicated strong enrichment of HOXB9 and MEIS1 regulon activity in Pax3+ progenitors; RAD21, EZH2, and CTCF activity in myogenic progenitors; TCF12, EGR1, and FOSB activity in myoblasts; and MYOD1, MEF2C, MYOG, and SOX6 activity in myocytes (Additional file 10: Table S9 and Fig. 5C). Although the activity of MEF2A and MEF2D, which belong to the same family as MEF2C, was elevated in myocytes, it was not as significant as MEF2C, indicating that they did not play a leading role in myogenesis. SCENIC also successfully inferred multiple downstream target genes. The complete list of regulons and their respective predicted target genes can be found in Additional file 11: Table S10. The scaled and binarized regulon activity is available in Additional file 12: Table S11. Examples of regulon activity, corresponding TF expression, and predicted target gene expression are depicted in Fig. 5D and Additional file 1: Figure S6A. TFs such as SOX6, TEAD4, and FOXO1 were predicted as targets of skeletal muscle-specific TF MYOG, and FOS was predicted as a target of FOSB, indicating a critical transcriptional hierarchy of skeletal muscle development. Corresponding chromatin accessibility in scATAC data for these TFs and predicted target genes are shown in Figure. S6B.

Chromatin dynamics of Pax3+ progenitor differentiation

The pseudotime trajectory in the scATAC-seq dataset was evaluated, which resulted in a similar cellular differentiation trajectory to scRNA-seq dataset (Fig. 6A). We integrated the scRNA-seq and scATAC-seq datasets to correlate and cross-validate gene expression profiles and the chromatin accessibility landscape in the myogenic cells using the Harmony algorithm [32]. The coembedded Uniform Manifold Approximation and Projection (UMAP) plots with cell type assignment for the scATAC-seq and scRNA-seq data suggested that the changes in chromatin accessibility and corresponding gene transcript expression in most myocytes occurred in a synchronous manner, whereas in other cell types, it was not fully synchronized, suggesting that other regulation is involved (Fig. 6B, C). Correlations between cell types of scATAC-seq and scRNA-seq cell types were computed with scRNA-seq variable genes (Fig. 6D). These results obtained from two independent approaches demonstrate that our two datasets are highly concordant and cross-validated.

Fig. 6figure 6

Integrated analysis of scATAC-seq and scRNA-seq data. A The pseudotime trajectory in the scATAC-seq dataset. B UMAP representation of scATAC-scRNA integration results. Cells are colored by technology (scATAC=red, scRNA=blue). C UMAP representation of scATAC-scRNA integration results. Cells are colored by cell type assignment. D Dot plot showing scATAC-scRNA integration cell type assignment using confusion matrix. Each column represents the original cell type assignment of scRNA-seq data, and each row represents the cell type assignment predicted after integration with scATAC-seq data. The size of the dots represents the absolute value of the correlation, and the red and gray dots represent the positive and negative correlations, respectively

We next performed pseudotime ordering of the chromatin accessibility-based TF motif enrichment of individual cells and correlated changes of the TF motif patterns with TF expression (Fig. 7A, B). To this end, we also investigated TFs and target genes differentially expressed over scRNA-seq pseudotime. We noticed good concordance of time-dependent changes in TF and predicted target gene expression along with motif enrichment, suggesting that a set of TFs cooperatively regulates myogenic differentiation (Fig. 7C and Additional file 1: Figure S7).

Fig. 7figure 7

Activity and RNA expression dynamics of TFs along the pseudotime trajectory. A Heatmap showing the activity of TFs along the differentiation trajectory. B TF expression heatmap corresponding to the motif enrichment along the differentiation trajectory. C Pseudotime-dependent chromatin accessibility and gene expression changes along the myogenic lineages. The first column shows the dynamics of the 10× Genomics TF enrichment score. The second column shows the dynamics of TF gene expression values, and the third and fourth columns represent the dynamics of SCENIC-reported target gene expression values of corresponding TFs. Error bars denote 95% confidence intervals of local polynomial regression fitting. Additional examples are given in Additional file 1: Figure S7

EGR1 and RHOB play critical roles in myogenesis

Although the roles of several identified TFs in myoblast specification and differentiation have been established, the expression and functions of many other genes that present specific expression profiles over scRNA-seq pseudotime have not been well studied. In consideration of the possibility that pig is used as a model of muscle disease in the future, and in order to identify the conserved functional genes related to muscle development among mammals, we compared the expression of these genes in skeletal muscle from wild-type and Duchenne muscular dystrophy (DMD) mice using an RNA-Seq dataset (GSE162455). Consistent with the patterns of classical myogenic genes, most of the genes upregulated along the pseudotime trajectory were expressed at higher levels in DMD mice than that in wild-type mice (Additional file 4: Table S3 and Additional file 1: Figure S8A, B), suggesting that these genes may be induced by muscle regeneration in muscular dystrophy mice and play an important role in myogenesis. The expression levels of genes downregulated along the pseudotime trajectory were also downregulated in DMD mice, suggesting that they are involved in muscular dystrophy and may not contribute to terminal myoblast differentiation (Additional file 4: Table S3 and Additional file 1: Figure S8C). The scRNA analysis highlighted that the expression of EGR1 and its predicted target gene RHOB gradually increased along the pseudotime trajectory, and scATAC analysis showed that EGR1 reached its strongest motif enrichment in myocytes (Fig. 7C). The peaks in the vicinity of EGR1 and RHOB, which were located at +3979 bp and +11,950 bp of the EGR1 and RHOB transcriptional start sites (TSSs) respectively, were the most accessible in myocytes (Additional file 1: Figure S9). Based on gene expression correlation and binding motif analysis, we identified a total of 215 high-confidence annotated TF-target pairs between 12 TFs and 75 target genes, and constructed a putative gene regulatory network associated with skeletal muscle development (Fig. 8A). Among them, the number of links between transcription factor EGR1 and myogenesis-related target genes is second only to the classical myogenic transcription factor MYOD1 (Fig. 8B). RHOB is an important regulator of cell and tissue morphology and function, acting mainly through the cellular cytoskeleton [33, 34]. A few studies have revealed that RHOB is a key mediator during diverse cellular and physiological processes like cell division, cell migration in smooth muscle cells [35, 36]. Totally, we speculated that EGR1 and RHOB are likely to play positive roles during myogenic differentiation in pigs, so they were selected for in vitro functional validation.

Fig. 8figure 8

Functional analysis of EGR1 and RHOB in myogenesis. A To explore the connection network between TFs and targets in skeletal muscle development, 172 genes associated with muscle development (http://wiki.geneontology.org/index.php/Muscle biology) were extracted as target genes. Then, 215 high-confidence annotation TF-target pairs were selected from regulon activity network to construct the regulatory network. B Target genes counts of TFs. C The mRNA levels of EGR1, RHOB, and myogenic markers during the differentiation of pig primary myogenic cells (PPMCs) at several indicated time points. When the cells were cultured in growth medium (GM) at sub-confluent density, it was defined as day 0 (day 0); when the cells reached 100% confluence, GM was changed to differentiation medium (DM). D The mRNA levels of EGR1, Myf5, and MyoD in proliferating C2C12 cells at 36 h after transfection with the pcDNA3.1-EGR1 vector. E Immunofluorescence staining for MyHC in PPMCs after transfection with the pcDNA3.1-EGR1 vector and differentiation induction for 5 days. The fusion index (the percentage of nuclei in fused myotubes out of the total nuclei) was calculated. F The mRNA levels of myogenic differentiation markers in C2C12 cells after transfection with the pcDNA3.1-EGR1 vector and induction of differentiation for 3 days. G The mRNA levels of RHOB, Myf5, and MyoD1 in proliferating C2C12 cells at 36 h after transfection with the pcDNA3.1-RHOB vector. H Immunofluorescence co-staining for Pax7 and MyoD in C2C12 cells transfected with the pcDNA3.1-EGR1 vector and cultured in growth medium for 36 h. I Statistical analysis was performed to quantify the percentages of the three myogenic cell populations. J Immunofluorescence staining for MyHC in PPMCs after transfection with the pcDNA3.1-RHOB vector and induced differentiation for 5 days. The fusion index was calculated. *P < 0.05, **P < 0.01, ***P < 0.001

To confirm the dynamic expression of EGR1 and RHOB during myogenic differentiation, quantitative PCR (qPCR) was performed on porcine primary myogenic cells (PPMCs) at several differentiation points (day 0, day 2, day 4, day 6, and day 8). It was found that EGR1 and RHOB had the same expression pattern as the well-known myogenic differentiation makers (Fig. 8C). To validate the effect of EGR1 and RHOB on myogenesis progression, they were overexpressed in PPMCs and mouse C2C12 myoblasts. Ethynyl-2′-deoxyuridine (EdU) incorporation and immunofluorescence assays showed that EGR1 overexpression did not influence cell proliferation but promoted myogenic differentiation, inducing a significant increase in the fusion index (Fig. 8D, E, and Additional file 1: Figure S10). In line with this, the expression of myogenic differentiation markers increased when EGR1 was overexpressed (Fig. 8F). Consistent with EGR1, RHOB overexpression did not change the proliferation ability of myoblasts, but the expression level of MYOD1 was significantly upregulated (Additional file 1: Figure S11A-C, and Fig. 8G). This prompted us to further verify whether RHOB regulates myoblast fate commitment. Immunofluorescence co-staining of PAX7 and MYOD showed that RHOB accelerated the transformation of PAX7+/MYOD− myogenic progenitors into PAX7+/MYOD+ myoblasts (Fig. 8H, I). After induction of differentiation, the RHOB overexpression group formed more myotubes with a significantly increased the fusion index (Fig. 8J, and Additional file 1: Figure S11D). These findings revealed that EGR1 and RHOB are critical regulators of pig embryonic myogenesis.

Cell–cell communications

To predict the cellular communications involved in pig skeletal muscle ontogeny, we evaluated the potential cell–cell interactions by using CellPhoneDB [37]. First, the cellular interactions between myogenic cells were analyzed, and it was found that the interactions involving Pax3+ progenitor were predicted to be more significant than those involving myogenic progenitor, myoblast, and myocyte (Fig. 9A). This result illustrated that the stronger the characteristics of stem cells, the greater the possibility of their interaction with other cells. Then, all somite cell populations, including mesenchymal cells, fibroblasts, epithelial cells, neural stem cells, osteogenic cells, neurons, neurogliacytes, endothelial cells, chondrocytes, and myogenic cells, were included in the analysis. Interestingly, the interaction between cells of the same cell type was weaker than that between different cell types. Myogenic cells tended to communicate wi

留言 (0)

沒有登入
gif