Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells

Abstract

Recent molecular studies have shown that, even when derived from a seemingly homogenous population, individual cells can exhibit substantial differences in gene expression, protein levels and phenotypic output, with important functional consequences. Existing studies of cellular heterogeneity, however, have typically measured only a few pre-selected RNAs or proteins simultaneously, because genomic profiling methods could not be applied to single cells until very recently. Here we use single-cell RNA sequencing to investigate heterogeneity in the response of mouse bone-marrow-derived dendritic cells (BMDCs) to lipopolysaccharide. We find extensive, and previously unobserved, bimodal variation in messenger RNA abundance and splicing patterns, which we validate by RNA-fluorescence in situ hybridization for select transcripts. In particular, hundreds of key immune genes are bimodally expressed across cells, surprisingly even for genes that are very highly expressed at the population average. Moreover, splicing patterns demonstrate previously unobserved levels of heterogeneity between cells. Some of the observed bimodality can be attributed to closely related, yet distinct, known maturity states of BMDCs; other portions reflect differences in the usage of key regulatory circuits. For example, we identify a module of 137 highly variable, yet co-regulated, antiviral response genes. Using cells from knockout mice, we show that variability in this module may be propagated through an interferon feedback circuit, involving the transcriptional regulators Stat2 and Irf7. Our study demonstrates the power and promise of single-cell genomics in uncovering functional diversity between cells and in deciphering cell states and circuits.

Figure 1. Single-cell RNA-Seq of LPS-stimulated BMDCs reveals extensive transcriptome heterogeneity

a–c, Correlations of transcript expression levels (x & y-axes: log-scale TPM+1) between two 10,000-cell population replicates (a), two single cells (b), and the ‘average’ single cell and a population (c). d,e, RNA-Seq read densities in single cells (blue) and population replicates (grey) for three non-variable genes (d) and four variable ones (e). f–g, RNA-FISH of representative transcripts. Optical micrographs (cell boundaries; grey outlines) and maximum-normalized distributions of expression levels from a RNA-FISH co-staining (n = 3,193 cells) for Il6 (yellow) and Cxcl1 (magenta).

Figure 2. Bimodal variation in expression levels across single cells

a, Relationship between average expression level in single cells (μ, X axis) and standard deviation (σ, Y axis) for 6,313 genes (Supplementary Table 2). Blue dashed line: maximum theoretical σ for an average expression level (SI); Grey dashed line: constant Fano factor (σ/μ = 0.25). Magenta: immune response genes; Green: housekeeping genes; light blue shaded region: single-cell average TPM < 250. b, Cellular heterogeneity for the 522 most highly expressed genes (single cell average; Supplementary Table 3). Each row represents a discretized histogram for a single gene (sorted by the Fano factor from low to high (top to bottom)). Color represents the number of cells (yellow: 18 cells; black: 0) that express the gene at the noted level. Grey dashed line denotes the constant Fano factor (0.25) highlighted in (a). c, Averaged expression density distributions for the 281 low-variability genes (top) and the 241 highly variable genes (bottom).

Figure 3. Variation in isoform usage between single cells

a, RNA-Seq read densities in single cells (blue) and population replicates (grey) for two illustrative loci, each with two different isoforms (bottom). b, Distributions of exon inclusion (PSI scores, X axis) for alternatively spliced exons of highly expressed genes (single-cell TPM > 250) in individual cells (blue histogram, top) and populations (grey histogram, bottom). c, Left: RNA-Seq read densities for Irf7 (only cells where the transcript is expressed are shown). Colored boxes mark exons analyzed by RNA-FISH. Right: RNA-FISH images from simultaneous hybridization with probes for two constitutive (‘Con’) regions of the transcript (A: cyan (C); B: magenta (M)) and one alternatively spliced exon (‘Specific’: orange (O)). White arrows (middle panel) highlight two cells with high levels of Irf7, but opposite preferences for the alternatively spliced exon. Histograms showing global abundance ratios for isoform-specific and constitutive probes (cells with less than 5 constitutive counts have been excluded; n = 490 cells; bottom histogram deviates from 0.5 due to probe design, see SI).

Figure 4. Analysis of co-variation in single-cell mRNA expression levels reveals distinct maturity states and an antiviral cell circuit

a, PCA of 632 LPS-induced genes. Contributions of each cell (points) to the first two principal components. b, Clustered correlation matrix of induced genes. Left: the Pearson correlation coefficients (r) between single-cell expression profiles of every pair of 632 LPS-induced genes (rows, columns). Right: the projection score (green: high; blue: low) for each gene (row) onto PC1 (left) and PC2 (right). c, Confirmation of correlations for Irf7-Stat2 (n = 655 cells) and Irf7-Ifit1 (n = 934 cells) by RNA-FISH. d–f, Expression levels for 16 genes in single BMDCs (columns), measured using single-cell qRT-PCR, in wild type (WT) (n = 36) (d), Irf7 −/− (n = 47) (e), and Ifnr −/− (n = 18) (f) at 4h after LPS stimulation (SI).