Single-Cell RNA Sequencing (scRNA-seq)

Placental samples were minced on ice into <1-mm³ pieces, followed by enzymatic digestion using trypsin. Subsequently, the solution was centrifuged at 300 rcf for 30 sec at room temperature and the supernatant was removed. 1 PBS (calcium- and magnesium-free) containing 0.04% weight/volume BSA (400 g/ml) was added to the supernatant, followed by centrifugation at 300 rcf for 5min. The cell pellet was then resuspended in 1ml red blood cell lysis buffer and incubated for 10min at 4C. After red blood cell lysis, samples were resuspended in 1ml PBS containing 0.04% BSA and filtered over SciencewareFlowmi 40-m cell strainers (VWR). Finally, cell concentration and viability were determined by hemocytometers and Trypan Blue staining.

The scRNA-seq libraries were prepared using Chromium Single Cell 3 Reagent v3 Kits. Single-cell suspensions were loaded on a Chromium Single Cell Controller Instrument (10X Genomics) to generate single-cell gel beads in emulsions (GEMs). Briefly, about 16,00020.000 cells were added to each channel, with a targeted cell recovery estimate of 5,0008,000 cells. After generation of GEMs, single-cell RNA-seq libraries were prepared using the Chromium Single Cell 3Library& Cell Bead Kit (10X Genomics) according to the manufacturers protocol. Libraries were sequenced with an IlluminaNovaseq6000 using high-output 75-cycle kits with apreviously reported read length configuration (21).

The Cell Ranger software pipeline (version 3.0) provided by 10X Genomics was used to demultiplex cellular barcodes, map reads to the genome and transcriptome using the STAR aligner, and down-sample reads as required to generate normalized aggregate data across samples, producing a matrix of gene counts versus cells. We processed the unique molecular identifier (UMI) count matrix using the R package Seurat (ver. 2.3.4). Low-quality cells (UMI/gene numbers out of the limit of mean value +/- 2 fold of standard and >10% mitochondrial genes) were excluded.

The top variable genes across single cells were identified using the method described in Macosko et al. (22). Briefly, the average expression and dispersion were calculated for each gene, and genes were subsequently placed into 22 bins based on expression. Principal component analysis (PCA) was performed to reduce the dimensionality in the log-transformed gene-barcode matrices of the top variable genes. Cells were clustered using a graph-based clustering approach and visualized in two dimensions using tSNE. We used likelihood ratio tests that simultaneously assessed changes in mean expression and in the percentage of expressed cells to identify significantly DEGs between clusters. We used the R package SingleR, a novel computational method for unbiased cell type recognition of scRNA-seq, with two reference transcriptomic datasets from the Human Primary Cell Atlas (23) to infer the cell of origin of each of the single cells independently, and to identify cell types.

Differentially expressed genes (DEGs) were identified using the FindMarkers function (test.use = MAST) in Seurat (24). P value < 0.05 and log2foldchange > 0.58 were set as the threshold for significantly differential expression. GO enrichment and KEGG pathway enrichment analysis of DEGs were performed using R software (R Development Core Team, Vienna, Austria) based on the hypergeometric distribution.

We determined the developmental pseudotime trajectory using the Monocle2 package (25). The raw count was first converted from a Seurat object into CellDataSet object with the importCDS function in Monocle2. We used the differential GeneTest function of the Monocle2 package to identify ordering genes (qval < 0.01) that were likely to be informative for the ordering of cells along the pseudotime trajectory. Dimensional reduction clustering analysis was performed with the reduce Dimension function, followed by trajectory inference with the order Cells function using default parameters. Gene expression was plotted with the plot gene in pseudotime function to track changes over pseudo-time.

We performed RNA velocity analysis using the R package velocyto.R (26) v0.6. The RNA velocity was calculated based on spliced and unspliced transcript reads and estimated using a gene-relative model. The resulting velocity estimates were projected onto the t-SNE embedding obtained in Seurat and the pseudotime space produced by Monocle 2.

We used CellPhoneDB (27) (v2.0) to identify biologically relevant ligand-receptor interactions from single-cell transcriptomics (scRNAseq) data. We considered a ligand or receptor to be expressed in a particular cell type if 10% of the cells of that type had non-zero read counts for the ligand/receptor encoding gene. Statistical significance was then assessed by randomly shuffling the cluster labels of all cells and repeating the above steps, which generated a null distribution for each LR pair in each pairwise comparison between two cell types. After running 1,000 permutations, P-values were calculated with the normal distribution curve generated from the permuted LR pair interaction scores. To define cell-cell communication networks, we paired ligand-receptor-expressing cell types. The R packages igraph and circlize were used to display the cell-cell communication networks.