Single-cell RNA sequencing data analysis

Fastq files were processed with Cell Ranger (version 3.1.0) (https://support.10xgenomics.com) on the DolphinNext pipeline builder [83] to generate the merged gene-barcode matrix. Matrix was loaded into R (version 3.6.3) and analyzed mainly with Seurat package ([84], version 3.1.5). Low-quality cells or outliers were filtered out from the data by the following threshold: 200 < total Gene count < 7000, 500 < total UMI count < 30000. Filtered matrices were normalized with SCTransform [85] and reduced the dimension with principal component analysis. The batch effect was corrected by harmony [86]. The significant principal components (PCs) are used to generate the UMAP plot for visualization. Clustering was performed by the Louvain algorithm based on K-nearest neighbor graph. Hurdle model tailored was used to perform differential expression analysis of the individual clusters [87]. Gene Ontology (GO) analysis was performed with clusterProfiler (version 3.18.0) [66].

For trajectory inference analysis, the cells from vegetative, slug, and fruiting stages were down-sampled and merged. The uninformative genes were filtered out by keeping the genes expressing 3 UMIs per cell for at least 10 cells. The filtered matrices were normalized with full quantile normalization and reduced the dimension with PCA. The first 3 PCs were used for k-means clustering and subsequent analysis. Slingshot wase performed to predict the multicellularity lineage [88]. The generalized additive model wse used to identify the temporally expressed gene through the pseudotime.