Quartz-Seq2–based scRNA-seq

We prepared two biological replicates for scRNA-seq analysis (fig. S9A). A total of 170 to 195 hypocotyl explants for each genotype (Col and wox13-2), which have been incubated on CIM for 4 days and then additionally incubated on SIM for 7 days, were treated with digestion buffer [600 mM mannitol, 2 mM MgCl₂, 0.1% bovine serum albumin, 2 mM CaCl₂, 2 mM MES, 10 mM KCl, 1.5% Cellulase RS, and 0.1% Pectolyase (pH 5.5)] for 1 hour at room temperature as described before (47). The protoplast solution was strained through a 70-μm filter. The filtered solution was centrifuged at 200g for 6 min and then washed once and lastly strained through a 40-μm filter. The protoplast pellet was resuspended in 0.5 ml of resuspension buffer. Cell viability was confirmed by PI staining. Alive singlet cells were captured by BD FACSAria II SORP (BD Biosciences, NJ, USA; fig. S9B) or cellenONE (Cellenion, Lyon, France) and collected into eight 384-well plates for each platform. Sequence library preparation of Quartz-Seq2 was performed for scRNA-seq analysis according to the previous study (26), and the libraries were sequenced by an Illumina NextSeq 500 High Output Kit v2.5 (75 cycles; Illumina, CA, USA). In calculation for digital expression matrix, we used genome file and annotation gtf file, which was provided by B. Cole (48). Flow cytometer information obtained from FACS and cell analysis data from cellenONE were analyzed together with the digital expression matrix. Cell analysis data obtained from cellenONE includes diameter information per isolated cell. Seurat v3.1.4 was used for the downstream data analyses including filtering, dimension reduction, cell clustering, and identification of marker genes. The UMI expression matrix with cells expressing at least 2000 genes was loaded into Seurat (Fig. 4A). A total of 3987 cells (1898 cells for WT and 2089 cells for wox13 mutant) with an average of 8373 genes were used for the following analyses. The dimensions of the expression matrix were then reduced by the RunPCA function, and the top 15 dimensions were used for FindNeighbors function and UMAP analysis. The cell clusters were identified by the FindClusters function with a resolution of 0.7. Marker genes for each cluster were identified using FindAllMarkers function. Because there was no obvious batch effect between biological replicates (fig. S9C), we did not perform the integration process.