FENG Xi-kang ,XIE Chun-di ,LI Yong-yao ,WANG Zi-shuai# ,BAI Li-jing#

1 School of Software,Northwestern Polytechnical University,Xi’an 710072,P.R.China

2 Institute of Animal Sciences,Chinese Academy of Agricultural Sciences,Beijing 100193,P.R.China

3 Shenzhen Branch,Guangdong Laboratory for Lingnan Modern Agriculture,Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs,Agricultural Genomics Institute at Shenzhen,Chinese Academy of Agricultural Sciences,Shenzhen 518000,P.R.China

Abstract Skeletal muscle regeneration is a complex process where various cell types and cytokines are involved. Single-cell RNA-sequencing (scRNA-seq) provides the opportunity to deconvolute heterogeneous tissue into individual cells based on their transcriptomic profiles. Recent scRNA-seq studies on mouse muscle regeneration have provided insights to understand the transcriptional dynamics that underpin muscle regeneration. However,a database to investigate gene expression profiling during skeletal muscle regeneration at the single-cell level is lacking. Here,we collected over 105 000 cells at 7 key regenerative time-points and non-injured muscles and developed a database,the Singlecell Skeletal Muscle Regeneration Database (SCSMRD). SCSMRD allows users to search the dynamic expression profiles of genes of interest across different cell types during the skeletal muscle regeneration process. It also provides a network to show the activity of regulons in different cell types at different time points. Pesudotime analysis showed the state changes trajectory of muscle stem cells (MuSCs) during skeletal muscle regeneration. This database is freely available at https://scsmrd.fengs-lab.com.

Keywords: scRNA-seq,skeletal muscle regeneration,database,regulon network,pseudotime

1.Introduction

Skeletal muscle comprises 40-50% of the human body mass (Huardet al.2002) and constitutes the largest metabolic and endocrine organ of the body. Skeletal muscle fitness and plasticity is an important determinant of human health and disease (Kent-Braunet al.1995).Previous studies have shown that the metabolic adaptive remodeling of skeletal muscle in response to environmental stimuli is an important regulator of the occurrence of metabolic diseases (Brass 1996;Stumpet al.2006;Argileset al.2016). Skeletal muscle dysfunction is closely related to more than 500 human diseases including muscular dystrophy,atrophy,type 2 diabetes,and aging-related sarcopenia (Ganet al.2018). Moreover,a recent study showed that proteasome stress in skeletal muscle can cause a long-term protective response,thereby delaying the aging of the retina and brain (Raiet al.2021).

In general,the turnover rate of adult skeletal muscle is about 1-2% of myonuclei per week (Schmalbruch and Lewis 2000). In daily life,skeletal muscle is susceptible to various injuries including mechanical trauma,thermal stress,myotoxic agents,ischaemia,neurological damage,and other pathogenic conditions (Beiner and Jokl 2001)which causes the influx of extracellular calcium (Jrvinenet al.2000) and eventually leads to the degradation of muscle proteins and necrosis (St Pierre and Tidball 1994). Then,muscle regeneration mediated by muscle stem cells is switched on. Skeletal muscle regeneration is a complex process involving multiple steps and various cell types (Turner and Badylak 2012). Therefore,a deeper mechanistic understanding of the gene regulatory networks during muscle regeneration might provide targets for regenerative medicine strategies,which can enhance or inducede-novoformation of functional skeletal muscle as a treatment for congenital absence or traumatic loss of tissue.

Meat production is an important productivity indicator and economic trait for livestock including pig,cattle,and sheep. The basic biological principles for meat production are the growth and development of skeletal muscle tissue.skeletal muscle tissue growth is a complex biological process,including the formation,proliferation,apoptosis,migration,fusion,and differentiation of myoblasts(Chal and Pourquie 2017). Gene expression changes during regeneration are likely to replicate the analogous events during initial embryonic development,suggesting regeneration is an ideal model to study skeletal muscle developmentin vivo(Goldman and Poss 2020). Thus,identifying important genes and specific regulatory networks that regulate myogenesis during skeletal muscle regeneration can provide new theoretical support for analyzing the molecular mechanism of meat production traits in livestock.

An individual skeletal muscle is a complex structure,composed of large multinucleated mature muscle cells(myofibers),muscle stem cells (MuSCs),fibro-adipogenic progenitors (FAPs),connective cells,nerve cells,immune cells,and the vasculature (Payloret al.2011). Previous studies have emphasized the importance of MuSCs (Relaix and Zammit 2012) and the need of other cell types and intercellular communication networks for successful muscle regeneration (Arnoldet al.2007;Villaltaet al.2011;Yang and Hu 2018). However,these research did not provide a full picture of the events and their dynamics. Recently,scRNA-seq has provided insights about the regeneration program of skeletal muscle and identified dynamic changes in the muscle stem cell,fibroblast,and immune cell populations occurring during muscle regeneration(Dell’Orsoet al.2019;Giordaniet al.2019;De Micheliet al.2020;Oprescuet al.2020). Therefore,given the importance of the multitude of cells and extrinsic factors for skeletal muscle regeneration,a database to visualize the single-cell transcriptional profile and dynamic changes in cell populations that may influence each other to effectively and rapidly regenerate muscle tissue is needed.

Here,we developed the single-cell skeletal muscle regeneration database (SCSMRD),the first web-based public database for accessing the skeletal muscle regeneration process at the single-cell level. The current version of SCSMRD contains 105 500 cells and 8-time points from no injury to 21 days post-injury. We have performed a series of analyses including cell clustering,cell type annotation,cell cycle identification,cell-typespecific regulons identification,and pseudotime analysis of skeletal muscle cells. These results are freely available for visualization and download at https://scsmrd.fengs-lab.com.

2.Materials and methods

The SCSMRD provides visualization of our analyzed results based on gene expression profiling during skeletal muscle regeneration at the single-cell level inside three perspectives,gene expression changes in different cell clusters along with different time points,cell-type-specific regulons,and gene expression changes in muscle cells along with pseudotime trajectory. The workflow of SCSMRD construction is demonstrated as Fig.1-A and the detailed description is presented in the following sections.

Fig.1 An overview of the SCSMRD construction. A,the 4 major steps for SCSMRD construction: (i) sample collection;(ii) data analysis;(iii) search function development;(iv) visualization &download. B,the UMAP visualization of cells grouped by dataset,cell type,cell date,and cell stage. C,two examples of the skeletal muscle cell-type-specific regulon network (MYOD1 and MYOG).D,the result of pseudotime analysis of skeletal muscle cells in our database.

2.1.Sample information

Our data set integrates single cells of skeletal muscle regeneration tissues at 8-time points including no injury(D0),0.5,2,3.5,5,7,10,and 21 days post-injury from previously published data (De Micheliet al.2020;Oprescuet al.2020),providing a total of 105 554 cells.All expression data were obtained from the NCBI Gene Expression Omnibus (GEO) with accession numbers:GSE143437,GSE143435 and GSE138826.

2.2.Data pre-processing and data analysis

Data pre-processing analysis including batch effects removal,normalization,PCA,and cell clustering,were mainly performed using the R package LIGER(https://github.com/welch-lab/liger) (Welchet al.2019)by following the standard guidelines. Briefly,several standard preprocessing steps with default parameters were firstly performed to normalize expression data to account for differences in sequencing depth and efficiency between cells,identify variably expressed genes,and scale the data so that each gene has the same variance. Subsequently,integrative non-negative matrix factorization was run on the normalized and scaled datasets with default parameters. Then,the resulting factors were used to jointly cluster cells with a resolution of 0.4. Gene markers for all clusters were identified using the “runWilcoxon” function,and markers which have padj(Benjamini-Hochberg adjustedP-value) less than 0.05 and logFC (log fold change between observations in groupvs.out group) larger than 3 were selected for cell type annotation analysis. Further,we manually assigned cell population identity based on cell-type-specific markers as previous described (Wosczyna and Rando 2018;Giordaniet al.2019). Cell cycle score were calculated as previously described (Fig.1-B). Specifically,we collected the gene sets reflecting 5 phases of the cell cycle (G1/S,S,G2/M,M,and M/G1) used ccording to Hanet al.(2008). We next calculated the correlation between the expression level of each gene and the average expression level of all genes in that phase-specific gene set and excluded genes with a correlation less than 0.3.The mean expression value of the remained gene sets was used as the score for that phase. Then,the phasespecific scores were normalized twice. First,each phase scores were centered and divided by their standard deviation. Second,we normalized the phase score across all the phases within each cell by centering and normalizing. The cells were assigned to a cell phase by their maximal phase scores. Cell-type-specific regulons were identified using IRIS3 (Maet al.2020),a web server for cell-type-specific regulons inference from scRNA-seq data for human and mouse (https://bmbl.bmi.osumc.edu/iris3/) (Fig.1-C). Pseudotime analysis of skeletal muscle cells were performed using Python package Scanpy Version 1.8 (https://scanpy.readthedocs.io/en/stable/tutorials.html#trajectory-inference) (Wolfet al.2018) with default parameters (Fig.1-D). And the linear regression model in webpage of SCSMRD was performed using the “jsregression.LinearRegression” model of JavaScript regression library with the default parameters (https://github.com/chen0040/js-regression). General codes for computational analysis follow the instructions of the respective software and customized modifications are available at https://github.com/xikanfeng2/SCSMRD.

2.3.Database implementation

In summary,the SCSMRD website was built combing the Apache HTTP server,Python programming language,HTML,JavaScript,and the MySQL Database. Clean single-cell gene expression profiling data were processed with Python scripts and finally was stored into the MySQL database and as text files. The front-end interface was developed based on the Bootstrap open-source toolkit with the version of 3.4.1 (https://getbootstrap.com/). The serverside application was written in Python 3.8,and the Django Python web framework with the version 3.2.3 (https://www.djangoproject.com) was selected as the model-viewcontroller (MVC) framework for this whole database. The web interactive visualization graphs were implemented using the Plotly Javascript open-source graphing library with the version 1.58.5 (https://plotly.com/javascript/).SCSMRD was published using the Apache2 HTTP server and is accessible at https://scsmrd.fengs-lab.com.

3.Results

To facilitate the utilization of SCSMRD,we designed a user-friendly interface to allow users to perform various operations. We provided four main functionalities: (1)Gene expression level search in different cell populations along with different time points;(2) cell-type-specific regulon network search;(3) dynamic gene expression changes along with the continuous myogenic cell-states search;and (4) database resource download.

3.1.Gene expression level search

For the “gene expression level search” function,our database allows users to visualize the search result in an interactive Boxplot (Fig.2-A and B),where the data can be grouped by cell type or date. Users can view the gene expression level of a gene of interest in different cell clusters by giving a gene symbol (e.g.,PAX7) in the search input box. Specifically,users can click the“Example” button to view the default gene expression level in SCSMRD. Furthermore,an interactive Boxplot of gene expression level grouped by different cell clusters(e.g.,cell type or cell date) is displayed by clicking the“Search” button. Users can select different grouping rules on the menu in the upper left corner to switch the Boxplot between different clustering rules. Clicking on the cluster name listed in the graph legend (e.g.,“D0” in Fig.2-A or “BT cells” in Fig.2-B) can remove the Boxplot data of this cluster. Moreover,double-clicking on a cluster name allows users to view the detailed gene expression value of this cluster. These interactive functions allow users to compare gene expression levels only for their clusters of interest. Finally,the Boxplot can be downloaded in PNG format for further usage. As demonstrated in Fig.2-A and B,we investigated the gene expression level of the genePAX7as an example for this function. Our results showed thatPAX7is specifically expressed in Musc (muscle stem cells) cells and has the least expression level at 2 days after injury. Given thatPAX7is a marker gene of Musc cells,our result is consistent with previous studies (Turner and Badylak 2012;Baghdadi and Tajbakhsh 2018).

Fig.2 The web interface of SCSMRD search functions. A and B,the result of “gene expression level search” page under the selection of PAX7 gene and the boxplot was grouped by cell type (A) or cell date (B). C,the result of “cell-type-specific regulon search” page under the selection of KLF4 regulator. D,the result of “gene dynamic expression along with the pseudotime trajectory search” page under the selection of SMOC2 gene.

3.2.Cell-type-specific regulon search

For “cell-type-specific regulon search” function,we designed an interactive Sankey diagram to present the cell-type-specific regulon network (Fig.2-C). Users can search regulator factors of interest by inputting the gene symbol in the search input box. As demonstrated in Fig.2-C,the result Sankey diagram is composed of 3 columns that are connected by lines. The middle column of the Sankey diagram is a group of genes controlled by the searched regulator. The left and right columns connected to the middle column are the total expression level of cells grouped by different cell types and different cell time points,respectively. Specifically,users can drag the items in each column to reorder the Sankey graph (e.g.,reordering the column according to the total expression level). Also,the expression level is reflected by the thickness of the connected lines and the detailed expression value can be viewed by moving the cursor on the specific line. In Fig.2-C,we took the regulator,KLF4,as an example for this function.KLF4is a transcription factor involved in regulating diverse cellular processes including cell growth,proliferation,and differentiation in multiple tissues (Ghaleb and Yang 2017). Our results showed thatKLF4regulon is activated in non-injured (D0)and after 5 days post-injury (from D5 to D21) stem cells including FAP and Muscs. Considering that D5 is a key time point for muscle regeneration where various cells are proliferating and differentiating,our result suggests thatKLF4may play important roles in skeletal muscle regeneration.

3.3.Gene dynamic expression along with the pseudotime trajectory search

For “gene dynamic expression along with the pseudotime trajectory search” function,we represented the expression level changes of genes along the continuous myogenic cell-states with an interactive Scatter plot including the linear regression line (Fig.2-D). Users can view the gene expression changes in muscle cells along with pseudotime trajectory. As presented in Fig.2-D,we used the gene,SMOC2,as an example to demonstrate this function. Our visualization shows that along the pseudotime trajectory the expression level ofSMOC2increases in the early stage and stabilizes in the end stage.

3.4.Database resource download

For the “database resource download” function,SCSMRD provides the download function of all processed data in the database,including: (1) the whole gene-cell expression level matrix which contains 105 554 cells from 8-time points;(2) the cell annotation matrix which contains the detailed information of each cell (e.g.,cell type,cell stage,and cell time point);and (3) the regulon network matrix which contains all the regulator and its regulated genes. These materials provide the possibility for researchers to conduct in-depth studies on single-cell skeletal muscle regeneration.The related data for database construction can be downloaded at https://scsmrd.fengs-lab.com/download/.

4.Discussion

One of the most fascinating questions in regenerative biology is why regenerative capacity is restricted to a subset of tissues. Skeletal muscle has a remarkable capacity to regenerate even after repeated traumas,making it a major focus of regeneration studies during evolution (Baghdadi and Tajbakhsh 2018). In addition,skeletal muscle is also an ideal tissue to investigate mechanisms underlying successful regeneration. Given skeletal muscle is a complex tissue structure composed of hundreds of cell populations and sub-populations.Therefore,our database provides a resource to investigate cell-type and gene-expression dynamics of skeletal muscle regeneration at single-cell resolution.

Recent studies of regenerating tissues have revealed that transcriptional networks of regeneration pathways are regulated in a context-specific manner to control key gene expression programs. A regulon is a group of genes controlled by the same repressor or activator gene (Jankyet al.2014). A clear assessment of regulons and the transcription factors (TFs) that control them during skeletal muscle regeneration is an effective strategy to pinpoint crucial and heterogeneous regulatory mechanisms encoded in diverse cell types,and those responsible for the injuries and degenerative diseases. In this database,we identified cell-type-specific regulons during skeletal muscle regeneration and visualized the activity of these TFs genes networks in different cell types and different time points.

Skeletal muscle regeneration is mediated by MuSCs.After the injury,MuSCs activate,proliferate,differentiate,and fuse together to repair damaged myofibers (Payloret al.2011). Thus,a clear dynamic gene regulatory program along the continuum of myogenic cellstates is critical to further understand the mechanisms underlying successful regeneration toward improving stem-cell-based therapies. In this study,we selected the MuSCs,myoblasts,and mature myocytes from our data and subsequently performed pseudotime analysis.Pseudotime analysis allows us to order cells from quiescent MuSCs to cycling and committed progenitors to mature muscle cells. Consequently,our database provides insights for the unbiased study of dynamic gene regulatory programs along the continuous myogenic cellstates during muscle regeneration processes.

Previous study have showed that expression patterns of tissue-specific genes including protein-coding RNA(Brawandet al.2011;Merkinet al.2012) and noncoding RNAs (both miRNA and lncRNA) (Meunieret al.2013;Necsulea and Kaessmann 2014) are highly conserved among mammals. The gene regulatory network constructed based on these gene expression data is also highly conserved among mammals (Guschanskiet al.2017;Chenet al.2019). Therefore,although the dynamic gene expression profiles and skeletal muscle specific network of regulons provided in SCSMRD were constructed based on mice data,it also provides a reference resource for studies focused on skeletal muscle development of livestock.

5.Conclusion

Our database provides insights into the cellular and molecular underpinning of skeletal muscle regeneration.We believe that this database will facilitate researchers within the fields of regenerative biology,skeletal muscle development biology,and meat production of livestock to zoom in on their study at the single-cell level.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (31972539 and 32102513),the Science,Technology,and Innovation Commission of Shenzhen Municipality,China (JCYJ20180306173644635),the Fundamental Research Funds for the Central Universities,China (G2020KY05109),the Natural Science Basic Research Program of Shaanxi Province,China(2022JQ-644),and the Basic Research Programs of Taicang,China (TC2021JC14).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Ethical approval

This study does not contain any studies with animal subjects performed by any authors.