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

Additional file 1: Figure S1. Quality control and batch effect correction in scRNA-Seq, related to Figure 1 A. Violin plots showing the number of expressed genes, the number of reads uniquely mapped against the reference genome, and the fraction of mitochondrial genes compared to all genes per cell in scRNA-Seq data. B. Box plot showing the number of genes (left) and the number of uniquely mapped reads (right) per cell in each identified cell type in scRNA-Seq data. C. tSNE plot visualization of the sample source for all 70,201 cells. Each dot is a cell. Different colors represent different samples. D. tSNE plot visualization of unsupervised clustering analysis for all 70,201 cells based on scRNA-Seq data after quality control, which gave rise to 31 distinct clusters. Figure S2. Gene Ontology (GO) analysis of the DEGs for each cell type was performed and the representative enriched GO terms are presented, related to Figure 1. Figure S3. Expression of selected marker genes along the differentiation trajectory, related to Figure 2 A. tSNE plot demonstrating cell cycle regression (left). Visualization of myogenic differentiation trajectory by cell cycle phases (G1, S, and G2/M) (right). B. Donut plots showing the percentages of cells in G1, S, and G2M phase at different cell states. C. Expression levels of cell cycle-related genes in the myogenic cells organized into the Monocle trajectory. D. Expression levels of muscle related genes in the myogenic cells organized into the Monocle trajectory. Figure S4. Unsupervised clustering analysis for all cells in scATAC-Seq data and myogenic-specific scATAC-seq peaks, related to Figure 4 A-C. tSNE plot visualization of the sample source for all 48514 cells in scATAC-Seq. Each dot is a cell. Different colors represent different pigs (A), different embryonic stages (B), or different samples (C). D. tSNE plot visualization of unsupervised clustering analysis for all 48514 cells after quality control in scATAC-Seq data, which gave rise to 15 distinct clusters. E. tSNE plot visualization of myogenic cells and other cells. Clusters 4 and 8 in Figure S4D were annotated as myogenic cells due to their high levels of accessibility of marker genes associated with myogenic lineage. F. Genome browser view of myogenic-specific peaks at the TSS of MyoG and Myf5 for myogenic cells and other cells in the scATAC-seq dataset. Figure S5. Percentage distribution of open chromatin elements in scATAC-Seq data, related to Figure 4 A. Distribution of open chromatin elements in each snATAC-seq sample. B. Distribution of open chromatin elements in snATAC-seq of myogenic cell types. C. Percentage distribution of open chromatin elements among DAPs in myogenic cell types. Figure S6. Integrative analysis of transcription factors and target genes, related to Figure 5 A. tSNE depiction of regulon activity (“on-blue”, “off-gray”), TF gene expression (red scale), and expression of predicted target genes (purple scale) of MyoG, FOSB, and TCF12. B. Corresponding chromatin accessibility in scATAC data for TFs and predicted target genes are depicted. Figure S7. Pseudotime-dependent chromatin accessibility and gene expression changes, related to Figure 7. 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 the SCENIC-reported target gene expression values of corresponding TFs, respectively. Figure S8. Myogenesis related gene expression in DMD (Duchenne muscular dystrophy) mice. Comparison of RNA-seq data of flexor digitorum short (FDB), extensor digitorum long (EDL), and soleus (SOL) in DMD and wild-type mice including 2- month and 5-month age. A. The expression levels of myogenesis related genes (Myod1, Myog, Myf5, Pax7). B. The expression levels of related genes that were upregulated during porcine embryonic myogenesis (EGR1, RHOB, KLF4, SOX8, NGFR, MAX, RBFOX2, ANXA6, HES6, RASSF4, PLS3, SPG21). C. The expression levels of related genes that were downregulated during porcine embryonic myogenesis COX5A, HOMER2, BNIP3, CNCS). Data were obtained from the GEO database (GSE162455; WT, n = 4; DMD, n = 7). Figure S9. Genome browser view of differentially accessible peaks at the TSS of EGR1 and RHOB between myogenic cells in the scATAC-seq dataset, related to Figure 8. Figure S10. Functional analysis of EGR1 in myogenesis, related to Figure 8 A-B. EdU assays for the proliferation of pig primary myogenic cells (A) and C2C12 myoblasts following EGR1 overexpression. C. qPCR analysis of the mRNA levels of cell cycle regulators in C2C12 cells following EGR1 overexpression. D. Immunofluorescence staining for MyHC in C2C12 cells following EGR1 overexpression and differentiation for 3 d. Then, the fusion index was calculated. Figure S11. Functional analysis of RHOB in myogenesis, related to Figure 8 A-B. EdU assays for proliferation of pig primary myogenic cells (A) and C2C12 myoblasts following RHOB overexpression. C. qPCR analysis of the mRNA levels of cell-cycle regulators in C2C12 cells following RHOB overexpression. D. Immunofluorescence staining for MyHC in C2C12 cells following RHOB overexpression and differentiation for 3 d. Then, the fusion index was calculated.

附加文件1：图S1。单细胞RNA测序（scRNA-Seq）中的质量控制与批次效应校正，对应图1。A. 小提琴图展示单细胞RNA测序数据中每个细胞的表达基因数、唯一比对至参考基因组的读段数，以及线粒体基因占所有基因的比例。B. 箱线图展示单细胞RNA测序数据中每种已鉴定细胞类型的每个细胞的基因数（左）和唯一比对读段数（右）。C. tSNE图可视化70201个全部细胞的样本来源，每个点代表一个细胞，不同颜色代表不同样本。D. 经质量控制后的单细胞RNA测序数据中，对70201个全部细胞进行无监督聚类分析的tSNE可视化结果，最终得到31个不同的聚类簇。
图S2。对每种细胞类型的差异表达基因（differentially expressed genes, DEGs）进行基因本体（Gene Ontology, GO）富集分析，展示代表性富集GO条目，对应图1。
图S3。分化轨迹中选定标记基因的表达情况，对应图2。A. tSNE图展示细胞周期回归（左图）；按细胞周期时相（G1、S和G2/M）可视化肌源性分化轨迹（右图）。B. 环形图展示不同细胞状态下G1、S和G2M时相的细胞占比。C. 按Monocle轨迹排列的肌源性细胞中，细胞周期相关基因的表达水平。D. 按Monocle轨迹排列的肌源性细胞中，肌肉相关基因的表达水平。
图S4。单细胞转座酶可及性测序（scATAC-Seq）数据的无监督聚类分析及肌源性特异性scATAC峰，对应图4。A-C. tSNE图可视化scATAC-Seq中48514个全部细胞的样本来源，不同颜色分别代表不同猪只（A）、不同胚胎发育阶段（B）或不同样本（C）。D. scATAC-Seq数据经质量控制后，对48514个全部细胞进行无监督聚类分析的tSNE可视化结果，最终得到15个不同的聚类簇。E. tSNE图可视化肌源性细胞与其他细胞：由于图S4D中的簇4和簇8具有肌源性谱系相关标记基因的高染色质可及性，因此将其注释为肌源性细胞。F. scATAC-seq数据集的肌源性细胞与其他细胞中，MyoG和Myf5基因转录起始位点（TSS）处肌源性特异性峰的基因组浏览器视图。
图S5。scATAC-Seq数据中开放染色质元件的百分比分布，对应图4。A. 每个单细胞核scATAC-seq样本的开放染色质元件分布。B. 肌源性细胞类型的单细胞核scATAC-seq开放染色质元件分布。C. 肌源性细胞类型中差异可及峰（differentially accessible peaks, DAPs）的开放染色质元件百分比分布。
图S6。转录因子与靶基因的整合分析，对应图5。A. tSNE图展示调控网络活性（“蓝色激活”“灰色沉默”）、MyoG、FOSB和TCF12的转录因子（transcription factor, TF）基因表达（红色标尺）及预测靶基因的表达（紫色标尺）。B. 展示转录因子及预测靶基因在scATAC数据中的对应染色质可及性。
图S7。拟时间依赖性染色质可及性与基因表达变化，对应图7。第一列展示10× Genomics转录因子富集得分的动态变化；第二列展示转录因子基因表达值的动态变化；第三、第四列分别展示对应转录因子经SCENIC分析得到的靶基因表达值的动态变化。
图S8。杜氏肌营养不良症（Duchenne muscular dystrophy, DMD）小鼠的肌发生相关基因表达。对比DMD小鼠与野生型小鼠（涵盖2月龄和5月龄）的指浅屈肌（FDB）、趾长伸肌（EDL）及比目鱼肌（SOL）的RNA-seq数据。A. 肌发生相关基因（Myod1、Myog、Myf5、Pax7）的表达水平。B. 猪胚胎肌发生过程中上调的相关基因（EGR1、RHOB、KLF4、SOX8、NGFR、MAX、RBFOX2、ANXA6、HES6、RASSF4、PLS3、SPG21）的表达水平。C. 猪胚胎肌发生过程中下调的相关基因（COX5A、HOMER2、BNIP3、CNCS）的表达水平。数据来源于GEO数据库（GSE162455；野生型组，n=4；DMD模型组，n=7）。
图S9。scATAC-seq数据集中肌源性细胞间EGR1和RHOB基因转录起始位点处差异可及峰的基因组浏览器视图，对应图8。
图S10。EGR1在肌发生中的功能分析，对应图8。A-B. EGR1过表达后，猪原代肌源性细胞（A）和C2C12成肌细胞增殖的EdU检测实验。C. EGR1过表达后，C2C12细胞中细胞周期调控因子的mRNA水平qPCR分析。D. EGR1过表达并诱导分化3天后，对C2C12细胞进行MyHC免疫荧光染色，随后计算融合指数。
图S11。RHOB在肌发生中的功能分析，对应图8。A-B. RHOB过表达后，猪原代肌源性细胞（A）和C2C12成肌细胞增殖的EdU检测实验。C. RHOB过表达后，C2C12细胞中细胞周期调控因子的mRNA水平qPCR分析。D. RHOB过表达并诱导分化3天后，对C2C12细胞进行MyHC免疫荧光染色，随后计算融合指数。