Bovine satellite cells were isolated from the longissimus muscle of a 2-week-old Holstein bull calf following euthanasia. Satellite cells were isolated using a procedure involving pronase digestion and differential centrifugation as described before (Hathaway et al., 1991). Satellite cells were cultured in growth medium consisting of Dulbecco’s modified eagle medium, 10% fetal bovine serum, 2 mM L-glutamine, and 1% antibiotics-antimycotics (Thermo Fisher Scientific, Waltham, MA, United States) for a week to remove dead cells and increase the number of viable cells prior to scRNA-seq. The animal-related procedure was approved by the Virginia Tech Institutional Animal Care and Use Committee.

Cells cultured above were detached from the culture plate with trypsin-EDTA (0.25%) and washed with resuspension buffer (phosphate-buffered saline, 0.04% bovine serum albumin). Cells were filtered through a 40 µm strainer to remove cell clumps. Cell viability was determined to be 90% by trypan blue staining and hemocytometer counting. Two scRNA-seq libraries were constructed from satellite cells to increase the number of cells sequenced. The scRNA-seq libraries were constructed using the Chromium Single Cell 3ʹ Gel Bead-in-Emulsion (GEM), Library & Gel Bead Kit v3 (10x Genomics, Pleasanton, CA, United States). We adjusted cell concentration to 1,000 cells/μl with resuspension buffer, and loaded 12.8 µl of diluted cell suspension with a master mix of reverse transcription (RT) reagent, template switch oligo, and RT enzyme into a Chromium Chip B (10x Genomics). Single-cell GEM generation, barcoding, and reverse transcription were achieved by running the Chromium Chip B on the Chromium Controller (10x Genomics). Specifically, single cells, RT reagents, gel beads containing barcoded oligonucleotides, and oil were combined on a microfluidic chip to form reaction nanovesicles. Within each reaction nanovesicle, a single cell was lysed, the gel bead was dissolved to free the identically barcoded RT oligonucleotides, and polyadenylated mRNA was then reverse transcribed into cDNA. Thus, all cDNAs from the same cell would have the same barcode, which would allow the sequencing reads to be traced back to their single cells of origin. cDNA was amplified for 11 cycles. Based on an Agilent High Sensitivity TapeStation analysis, both scRNA-seq libraries contained a single peak of DNA between 300 and 700 bp, with an average fragment size of 440 bp.

Each scRNA-seq library was pair-end sequenced in a single cell flow lane on an Illumina HiSeq system at the Novogene-UC Davis Sequencing Center (Novogene, Sacramento, CA, United States). Sequencing reads were processed and analyzed using the 10x Genomics Cell Ranger 3.0.2 software, which was composed of different analysis pipelines (Zheng et al., 2017). Specifically, sequencing reads were de-multiplexed using the cellranger mkfastq pipeline and aligned to the bovine genome and transcriptome (Bos_taurus.ARS-UCD1.2) at default parameters using the cellranger count pipeline. Uniquely mapped sequences from each library were used for unique molecular identifiers (UMI) counting using the cellranger count pipeline. The output files from the cellranger counting of reads from two scRNA-seq libraries were combined using the cellranger aggr pipeline.

Cells were initially clustered by expression similarity using the graph-based clustering algorithm, which consisted of building a sparse nearest-neighbor graph followed by Louvain Modularity Optimization, an algorithm that sought to find highly-connected modules in the graph (Blondel et al., 2008). Cells were further clustered by an additional cluster-merging step, which included hierarchical clustering on the cluster-medoids in principal components analysis (PCA) space and merging pairs of sibling clusters if there were no genes differentially expressed between them (with B-H adjusted p-value below 0.05). The hierarchical clustering and merging was repeated until there were no more cluster-pairs to merge. Differential gene expression between cell clusters was identified using the quick and simple method sSeq and edgeR (Robinson and Smyth, 2007; Yu et al., 2013). All these analysis pipelines and algorithms were part of the Cell Ranger software and run using parameters recommended by 10X Genomics.

Cell clusters were annotated based on significantly enriched expression (with B-H adjusted p-value below 0.05) of marker genes of cell types. Cell clusters and gene expression data were visualized in the Loupe Cell Browser (10x Genomics). Gene ontology analysis was performed with the PANTHER classification system (Mi et al., 2013). Gene ontology terms with corrected p-values less than 0.05 were considered significantly enriched.

Trajectory inference was performed using the Monocle algorithm (version 2.20.0) (Trapnell et al., 2014; Qiu et al., 2017). Briefly, the output files of the cellranger count pipeline were read in to generate the count matrix using the Read10X function of Seurat 4 (Hao et al., 2021). A Seurat object was then built from the count matrix using the CreateSeuratObject function of Seurat. The Seurat object was converted to a CellDataSet modeled with the negative binomial distribution using the newimport function of Monocle. Genes with mean expression levels >0.1 were selected to define the state of cells. Following a dimensionality reduction using the reduceDimension function, cells were ordered along the trajectory using the orderCells function of Monocle.

Scripts and codes used to run the Cell Ranger, Seurat, and Monocle programs in this study were adopted from the manuals and websites for these programs and are listed in the Supplementary File S1.

Sequencing of the two scRNA-seq libraries of bovine satellite cells in culture generated approximately 430 million reads per library. More than 85 and 50% of these reads were uniquely mapped to the bovine genome and transcriptome (bos_taurus.ars-ucd1.2), respectively (Table 1). Based on the Cell Ranger analyses, these reads represented the transcriptomes in more than 9,000 cells for each library (Table 1). The mean number of reads detected per cell was more than 43,000 for both libraries, and the median number of genes detected per cell was over 3,500 for both libraries (Table 1). These numbers met or exceeded the minimum requirements for a quality scRNA-seq analysis (Handley et al., 2015; Haque et al., 2017; Ziegenhain et al., 2017). Because the two scRNA-seq libraries were prepared from the same cells, sequencing data from the two libraries were combined to increase the total number of cells analyzed. Clustering the transcriptomes of 19,096 cells combined from the two scRNA-libraries revealed 15 cell clusters that differed in gene expression pattern (Figure 1, Supplementary File S2). These 15 clusters contained 300–2,500 cells, or 2–13% total cells analyzed, with cluster 1 being the largest cluster and cluster 15 being the smallest cluster (Table 2).

Based on the expression levels of marker genes (MYOD1, MYF5, and DES) of myoblasts (Pownall et al., 2002), clusters 1, 2, 3, and 12 were characterized as subsets of myoblasts, which are activated and proliferating satellite cells (Figures 1A, 2A). Expression of MYF5 mRNA in cluster 1, expression of MYOD1 mRNA in cluster 2, and expression of both MYF5 and MYOD1 mRNAs in cluster 3 were significantly enriched (Supplementary File S2, Figure 2A). PAX7 and PAX3 are markers of satellite cells (Maroto et al., 1997; Seale et al., 2000). Expression of PAX7 mRNA was significantly enriched in cluster 4, and PAX3 mRNA expression was significantly enriched in clusters 5 and 8 (Supplementary File S2). Because they were enriched with PAX7 or PAX3, clusters 4, 5, and 8 were determined to be subsets of satellite cells (Figure 1A).

MYOG is a master transcriptional regulator of myoblast differentiation (Pownall et al., 2002). MYOG mRNA expression was significantly enriched in clusters 3 and 12 (Supplementary File S2, Figure 2A). Besides MYOG, many muscle-specific genes such as MB, MYH3, MYL1, NEB, and STAC3 and several myoblast differentiation and fusion regulatory genes such as MEF2A, MEF2D and MYMK (Millay et al., 2013; Estrella et al., 2015) were upregulated in clusters 3 and 12 (Supplementary File S2, Figure 2A). These two clusters clearly contained differentiating or differentiated myoblasts, i.e., myocytes. Between clusters 3 and 12, more muscle-specific genes were upregulated in cluster 12 than in cluster 3, and the same muscle-specific genes were expressed at greater levels in cluster 12 than in cluster 3 (Supplementary File S2, Figure 2A). These differences suggest that myoblasts in cluster 12 were more terminally differentiated than those in cluster 3.

Gene ontology analyses of genes upregulated in each cluster indicated that cells in different clusters differed in function. For example, gene ontology analyses of genes upregulated in clusters 2 and 4 indicated that many of these genes were involved in the biological processes, cellular components, and molecular function related to DNA synthesis and cell cycle (Supplementary File S3), suggesting that cells in these clusters were undergoing active proliferation. Gene ontology analyses of genes upregulated in clusters 3 and 12 indicated that many of these genes were involved in the biological processes, cellular components, and molecular function related to mature skeletal muscle structure and contraction (Supplementary File S3), suggesting that cells in these two clusters were differentiating into functional muscle cells.

A trajectory inference analysis using the Monocle program (Trapnell et al., 2014) revealed the potential lineage relationships between the seven myogenic cell clusters, i.e., clusters 1–5, 8, and 12 (Figure 2B). This analysis suggested two trajectories along which the myogenic progenitors PAX3+ satellite cells transitioned toward myogenic differentiation. On one trajectory, PAX3+ satellite cells in cluster 8 committed to MYF5+ myoblasts in cluster 1, and these MYF5+ myoblasts then differentiated into MYOG+ myocytes in cluster 3, which then differentiated further into MYOG+ myocytes in cluster 12 (Figure 2B). On the other trajectory, PAX3+ satellite cells in cluster 8 first transitioned to a population of satellite cells in cluster 5 that had a different gene expression pattern from cells in cluster 8 but were still PAX3 positive; PAX3+ satellite cells in cluster 5 then became PAX7+ satellite cells in cluster 4, which were subsequently activated to become MYOD1+ myoblasts in cluster 2 (Figure 2B). Overall, this trajectory analysis further supports the conclusion earlier that cultured bovine satellite cells are composed of subsets of myogenic cells that differ in transcriptional and myogenic state.

Platelet derived growth factor receptor alpha (PDGFRA) is established as a marker of the progenitor cells for intramuscular adipocytes, i.e., intramuscular fibro-adipogenic progenitor (FAP) cells, in mice and humans (Uezumi et al., 2010; Uezumi et al., 2014; Joe et al., 2010). The expression of PDGFRA mRNA was significantly higher in clusters 7 and 9 than in other clusters (Supplementary File S2 and Figure 2A). It is interesting to note that cells in cluster 9 were also enriched with peroxisome proliferator activated receptor gamma (PPARG) and zinc finger protein 423 (ZNF423) mRNAs, two transcriptional regulators of early adipogenesis (Tontonoz and Spiegelman, 2008; Gupta et al., 2010; Lefterova et al., 2014). Gene ontology analyses of genes upregulated in clusters 7 and 9 indicated active lipogenesis in cluster 9 (Supplementary File S3). Cells in neither cluster 7 nor cluster 9 expressed markers of mature adipocytes such as leptin (LEP) and adiponectin (ADIPOQ) (Supplementary File S2). These data together indicated that clusters 7 and 9 were FAPs, or intramuscular preadipocytes, with cluster 9 appearing to be more developed preadipocytes than cluster 7. Because cells in clusters 7 and 9 were not enriched with MYOD1, MYF5, DES, or MYOG, markers of myogenic cells, it remains to be determined if these FAP cells were derived from satellite cells during culture or accidently co-isolated with satellite cells from skeletal muscle.

Compared to other clusters (Supplementary File S2), clusters 6, 10, 11, 13, 14, and 15 did not express significantly higher levels of markers of myogenic cells such as PAX3, PAX7, MYOD1, MYF5, and MYOG (Pownall et al., 2002); thus, these clusters were not myogenic cells. Skeletal muscle contains not only myogenic cells but also nonmyogenic cells such as endothelial cells, pericytes, smooth muscle cells, fibroblasts, and glial cells. However, none of clusters 6, 10, 11, 13, 14, and 15 appeared to be these nonmyogenic cells based on the expression levels of marker genes such as CDH5 and PECAM1 for endothelial cells (Elmentaite et al., 2021), NOTCH3 and MCAM for pericytes (Elmentaite et al., 2021), ACTA2 and MYH11 for smooth muscle cells (Kumar et al., 2017), COL1A and S100A4 for fibroblasts (Kumar et al., 2017), and FOXD3 and SOX10 for glial cells (Elmentaite et al., 2021) in these clusters (Supplementary File S2). It is possible that clusters 6, 10, 11, 13, 14, and 15 represented novel cell types in bovine skeletal muscle co-isolated with satellite cells or that they were inaccurately clustered by the computational program used.