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

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.

To characterize the extent of expression variability on a genomic scale and decipher its functional implications, we used single-cell RNA-Seq to profile a temporal snapshot of the BMDC response to LPS. This is an attractive model system for single-cell analyses for several reasons. First, LPS, a component of gram-negative bacteria and a ligand of Toll-like receptor 4, strongly synchronizes cellular responses and mitigates temporal phasing[11]. Second, LPS activation evokes a robust transcriptional program that has been extensively investigated at the population level[12]. Third, LPS stimulation should increase the correlation between mRNA and protein levels for induced genes, thus reducing a potentially confounding factor[13]. Lastly, differentiated BMDCs are post-mitotic, largely removing cell cycle-dependent transcriptional variation[3].

We stimulated BMDCs with LPS and harvested single cells after four hours[12] (Supplementary Information (SI)). Using SMART-Seq[9], we constructed cDNA libraries from 18 single BMDCs (S1–S18), three replicate populations of 10,000 cells, and two negative controls (empty wells), and sequenced each to an average depth of 27-million read-pairs. Negative control libraries failed to align (<0.25%) to the mouse genome, and were discarded from all further analyses. Library quality metrics, such as genomic alignment rates, rRNA contamination, and 3′ or 5′ coverage bias, were similar across all libraries (Supplementary Table 1). We estimated expression levels for all UCSC-annotated genes using RSEM[14] (Supplementary Table 2, SI) and discarded genes that were not appreciably expressed (transcripts per million (TPM) > 1) in at least three individual cells, retaining 6,313 genes for further analysis.

While the gene expression levels of population replicates were tightly correlated with one another (Pearson r > 0.98, log-scale, Fig. 1a), there were substantial differences in expression between individual cells (0.29 < r < 0.62, mean: 0.48, Fig. 1b, Supplementary Fig. 1). Despite this extensive cell-to-cell variation, expression levels for an “average” single cell correlated well with the population samples (0.79 < r < 0.81, Fig. 1c, Supplementary Fig. 1).

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).

We used RNA-FISH, an amplification-free imaging technique[2], to verify that heterogeneity in our single-cell expression data reflected true biological differences, rather than technical noise associated with the amplification of small amounts of cellular RNA. For 25 genes, selected to cover a wide range of expression levels, the variation in gene expression detected by RNA-FISH closely mirrored the heterogeneity observed in our sequencing data (Fig. 1d–g, Supplementary Fig. 2). For example, expression of housekeeping genes (e.g., Beta-Actin (Actb), Beta-2-microglobulin (B2m)) matched a log-normal distribution in both single-cell RNA-Seq and RNA-FISH measurements, consistent with previous studies[1]. In contrast, many genes involved in the LPS response, although highly expressed on average, showed significantly greater levels of heterogeneity, with expression levels deviating ~1,000 fold between individual cells in extreme cases (Fig. 1e–g).

More generally, we observed that single cell variability existed across a wide range of population expression levels (Fig. 2a). Of the 522 most highly expressed genes (single-cell average TPM > 250, Fig. 2a: unshaded region, Supplementary Table 3), 281 had low cell-to-cell variability (coefficient of variation (CV) < 0.25, SI) and were well described by log-normal distributions (RNA-Seq: Fig. 2b,c top, RNA-FISH (Actb, B2m): Supplementary Fig. 2). These 281 genes were enriched for housekeeping genes, encoding ribosomal and other structural proteins (Supplementary Table 2 & 3, Bonferroni-corrected p=1.5×10−6), consistent with previous findings in yeast[15] and mammalian cells[1].

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).

Surprisingly, however, 185 of the remaining 241 (coefficient of variation (CV) > 0.25, SI) highly expressed genes had bimodal expression patterns (Fig. 2b,c bottom): mRNA levels for these genes were high in many of the cells, but were at least an order of magnitude lower (often very low or undetectable) than the single-cell average in three or more cells. We independently verified this disparity by RNA-FISH (e.g., Cxcl1, Cxcl10, Ifit1, and others: Fig. 1f,g & Supplementary Fig. 2), confirming that it was not a result of technical noise. This variable set included both antiviral and inflammatory response genes, and was highly enriched for genes whose expression was increased by at least two-fold upon LPS stimulation at the population level[16] (p = 2.7×10−7; hypergeometric test; Supplementary Table 2). Still, bimodal expression was not a universal feature of immune response transcripts; some key chemokines and chemokine receptors (Ccl3, Ccl4, Ccrl2), cytokines (Cxcl2), and signaling molecules (Tank) were highly expressed in every cell (Supplementary Fig. 3), indicating that all cells were indeed activated by LPS.

This degree of variation in expression for highly expressed (on average) transcripts has not been observed in previous reports[7-10]. For example, examination of published single-cell RNA-Seq datasets of human embryonic stem cells[9] (Fig. 2a), mouse embryonic stem cells, and terminally differentiated fibroblasts[10] (Supplementary Fig. 4) revealed far less heterogeneity in expression for highly abundant (population average) genes. Similarly, studies of protein expression in mid-log yeast cells and dividing human cell lines[15,17] did not find such bimodality in (on average) highly expressed genes. We thus hypothesized that widespread variability in single-cell gene expression may reflect functionally important differences in the stimulated BMDC population.

Furthermore, we found that splicing patterns also showed previously unobserved levels of heterogeneity across single cells. Specifically, for genes that have multiple splice isoforms at the population level, individual cells predominantly expressed one particular isoform. We calculated the frequency (percent spliced in, PSI) of previously annotated splicing events in each of our samples using MISO[18], a Bayesian framework for calculating isoform ratios (Supplementary Table 4). Although the population-derived estimates were highly reproducible, single cells exhibited significant variability in their exon-inclusion frequencies (Fig. 3a,b).

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).

We considered the possibility that PCR amplification (intrinsic to the library preparation process) could potentially produce an overestimation of isoform regulation variability, particularly for weakly expressed transcripts[19]. However, even when we limited our analysis to 89 alternatively spliced exons (0.2 < population PSI < 0.8) that were very highly expressed within a single cell (single cell TPM > 250, SI), we still observed the same variability in splicing patterns amongst individual cells, with highly skewed expression towards a single splice variant (Fig. 3b). We obtained similar results when we generated three additional single-cell cDNA libraries using a slightly modified SMART-Seq protocol (SI) in which a four nucleotide barcode was introduced onto each RNA molecule during reverse transcription[19], enabling us to estimate the number of unique RNA transcripts that existed prior to PCR (Supplementary Fig. 5 & 6 and SI).

To the best of our knowledge, single-cell variation in splicing patterns has rarely been studied for individual genes, and never been analyzed on a genomic scale. One recent report[20] used RNA-FISH to study variation in alternative isoforms in two genes, and observed lower levels of isoform variability across single cells (the levels of heterogeneity differed in different cell types). Another study that used fluorescent reporters to quantify single-cell exon inclusion levels for one gene discovered highly variable and bimodal splicing patterns[21].

To independently verify the existence of extensive differences in isoform ratios between cells, we designed RNA-FISH probes targeting constitutive and isoform-specific exons in two genes (Irf7 and Acpp, Fig. 3c and Supplementary Fig. 7 & 8)[20]. We found substantial expression variability in overall Irf7 levels between individual cells (as reflected by the ‘constitutive’ probes, Fig. 3c, bottom and top panels), mirroring our single-cell sequencing results (and further explored below). Additionally, within each Irf7-expressing cell, we observed a bias towards either the inclusion or exclusion of the cassette exon (Fig. 3c, Supplementary Fig. 7, middle panel, e.g., compare ‘high’ and ‘low’ marked cells). We obtained comparable results for Acpp using two probes designed to detect mutually exclusive alternative final exons (Supplementary Fig. 8).

We next explored the sources and functional implications of expression variability. Bimodality amongst highly expressed immune response genes may reflect either the presence of distinct cellular subtypes or stochastic differences in the activation of signaling circuits[11]. We performed a principal components analysis (PCA, Fig. 4a) on our single-cell expression profiles, focusing on the 632 genes that were induced at least twofold in the population-wide response to LPS[16] (Supplementary Table 5). We found two distinct subpopulations, clearly distinguishable by the first principal component (PC1, 15% of the total variation, Fig. 4a). One group of fifteen cells expressed a core set of antiviral and inflammatory defense cytokines (including: Tnf, Il1a, Il1b, and Cxcl10) at extremely high levels (TPM > 1,000), while the remaining three cells expressed them at far weaker levels (TPM < 50). Some cell surface proteins (Ccr7, Cd83) and chemokines (Ccl22), which are known markers of BMDC maturation, showed the opposite expression pattern (Fig. 4b, Supplementary Fig. 9).

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).

During maturation, BMDCs switch from antigen-capturing to antigen-presenting cells that prime the adaptive immune system[22]. Maturation can occur either in response to pathogen-derived ligands (pathogen-dependent maturation), such as LPS, or when clusters of BMDCs are disrupted in culture[22] (pathogen-independent maturation). Both processes lead to induction of maturation markers, but only pathogen-dependent maturation results in co-expression of defense cytokines.

Examining the expression of maturation markers and defense cytokines (Supplementary Fig. 9) suggested that our 18 cells represent two distinct maturity states: (1) fifteen cells that were in the early stages of pathogen-dependent maturation (Fig. 4a, ‘maturing’, triangles; grey triangles, the two cells furthest along in this process); and, (2) three cells that likely matured during the culturing process (Fig. 4a, ‘mature’, squares; pathogen-independent). We further verified the existence of these sub-populations via RNA-FISH (Supplementary Fig. 10), single-cell quantitative reverse transcription polymerase chain reaction (qRT-PCR; Supplementary Fig. 11, SI, Supplementary Table 6), and cell sorting based on surface markers identified from the RNA-Seq data (Supplementary Fig. 12, SI). These results highlight that single-cell RNA-Seq can sensitively distinguish between closely related, yet distinct, developmental states, even within the same cell type.

Since differences in cell state explain only a small portion of the observed heterogeneity, we next examined the variation that might arise from the differential activity of regulatory circuits. We reasoned that co-variation across single cells between the mRNA levels of a transcription factor and its targets would represent a potential regulatory interaction, and, furthermore, would suggest that heterogeneity in the regulator’s expression may underlie the variability of its targets. Such a correlative approach has successfully identified regulatory connections from population-level transcription profiles measured in different conditions[12,23]. Here, we attempted to apply it to multiple single cells in the same condition.

To this end, we calculated the correlation in expression profiles between every pair of induced genes across all single cells, and identified a cluster of 137 genes that varied in a correlated way and were strongly discriminated by the second principal component (PC2, 8% of the variation, Fig. 4a,b). The cluster’s genes included the known antiviral master regulators Irf7 and Stat2, and were highly enriched for members of the antiviral response[12] (60 of 137 genes, p = 2.5×10−3, hypergeometric test, Supplementary Table 5), as well as STAT2 targets[16] (73/137 genes, p = 4.5×10−5, hypergeometric test). Most (100/137) of the cluster’s genes were bimodally expressed across single cells (Fig. 2c, bottom) despite being strongly expressed at the population level (13 genes TPM > 250; 53 genes TPM > 50). We independently validated a subset of these correlations using single-cell qRT-PCR and RNA-FISH (Fig. 4c,d). Moreover, single-cell qRT-PCR analysis of additional time points demonstrated that these correlations persisted at 6h as well (Supplementary Discussion in SI, Supplementary Fig. 13).

We hypothesized that bimodal variation in the expression of the cluster’s genes may be related to differences in the levels and activities of Stat2 and Irf7. To test this hypothesis, we measured expression of a set of antiviral genes by single-cell qRT-PCR in LPS-stimulated BMDCs from Irf7 knockout (Irf7 −/−) mice (SI). As expected, this perturbation ablated expression of most of the variable antiviral transcripts in our signature, while leaving non-variable antiviral transcripts relatively unaffected (Fig. 4e). However, Stat2 expression and variability levels were unaffected by the Irf7 knockout, implying that Stat2 may act either upstream or in parallel to Irf7 during the response[24] (Supplementary Fig. 14). As both Stat2 and Irf7 are targets of the interferon-signaling pathway, we stimulated and profiled BMDCs from interferon receptor knockout (Ifnr −/−) mice. In these cells, we found drastically reduced expression for both Stat2 and Irf7, as well as all other measured cluster genes (Fig. 4f).

Our analysis provides a proof-of-concept demonstrating how co-variation between transcripts across seemingly homogeneous single cells can help to identify and assemble regulatory circuits. Specifically, in our variable circuit (Supplementary Fig. 14) interferon signaling is required for induction of Stat2 and Irf7, which, in turn, act to induce our variable antiviral cluster genes. Our experiments do not definitively determine, however, which component of the circuit causes the observed heterogeneity per se. One compelling possibility is that upstream noise is propagated from the interferon-signaling pathway first to Stat2 and Irf7 and then to the target genes[25,26]. This hypothesis is supported by the variation we observed in STAT protein levels and nuclear localization (Supplementary Discussion in SI, Supplementary Fig. 15 & 16). However, since temporal snapshots of RNA and protein are not always directly comparable (Supplementary Discussion in SI, Supplementary Fig. 15 & 16), new strategies for tracing the spatiotemporal dynamics of both proteins and RNA in single living cells are needed to fully test this hypothesis[11].

A similar approach could potentially be used to explore the consequences of bimodality in splicing. Even looking at just 18 cells, we witnessed interesting examples of bimodal splicing patterns for genes whose isoforms have distinct functional consequences. For example, the splicing regulators Srsf3 and Srsf7 are each known to contain a “poison cassette exon” that, when included, targets the RNA for degradation via nonsense-mediated decay[27] (Supplementary Fig. 17). Meanwhile, splicing differences in other regulatory genes may further enhance expression diversity: for example, proteins encoded by different isoforms of Irf7 (Fig. 3c) differentially activate interferon-responsive genes in vitro[24]. These examples suggest that heterogeneity in splicing may represent another layer of response encoding.

In conclusion, our study reveals extensive bimodality in the transcriptional response of BMDCs to LPS, reflected in gene expression, alternative splicing, and regulatory circuit activity. While some variation in expression reflects differences in developmental state, other bimodal patterns reflect the differential activity of an antiviral regulatory circuit in this temporal snapshot. These phenomena allowed us to treat each cell as a “perturbation system” for reconstructing cell circuits[28], even with relatively few cells.

Moreover, our results demonstrate how co-variation across single cells can help dissect and refine gene modules that may be indistinguishable in population-scale measurements. For instance, in a recent population-scale study[16], we identified a large cluster of 808 “late-induced” LPS genes that was enriched for both maturation genes and STAT-regulated antiviral genes. These two subsets could not be separated by population-level expression profiles alone[16], but our single-cell data from a single timepoint clearly distinguishes them. Similarly, the unexpected and prevalent skewing we discovered in alternative splicing between single cells revises our molecular view of this process.

Finally, although many of our analyses focused on highly expressed genes to reduce the potential influence of amplification noise, our data also revealed substantial bimodality amongst more moderately expressed transcripts, such as large non-coding RNAs (lincRNAs, Supplementary Fig. 18). This suggests that the low population-level expression of these transcripts[29] may sometimes reflect high expression in a small subset of cells as opposed to uniform levels of low expression. While further technical improvements will be necessary to disentangle these two hypotheses (Supplementary Fig. 5), single-cell measurements should help facilitate the discovery and annotation of lincRNAs.

Comparing our results to other single-cell RNA-Seq data sets (e.g., Fig. 2a, Supplementary Fig. 4) indicates that the source of the analyzed tissue (in vitro vs. ex vivo), the biological condition of the individual cells (steady state vs. dynamically responding), and the cellular microenvironment all likely influence the extent of single-cell heterogeneity within a system. When applied to complex tissues – such as unsorted bone marrow, developing embryos, tumors, and other rare clinical samples – the variability seen through single-cell genomics may help determine new cell classification schemes, identify transitional states, discover previously unrecognized biological distinctions, and map markers that differentiate them. Fulfilling this potential would require novel strategies to address the high levels of noise inherent in single-cell genomics – both technical, due to minute amounts of input material, and biological, e.g., due to short bursts of RNA transcription[30]. Future studies that couple technological advances in experimental preparation with novel computational approaches would enable analyses, based on hundreds or thousands of single cells, to reconstruct intracellular circuits, enumerate and redefine cell states and types, and transform our understanding of cellular decision-making on a genomic scale.

Methods Summary

BMDCs, prepared as previously described[12], were stimulated with LPS for 4h and then sorted as single cells or populations (10,000 cells) directly into TCL lysis buffer (Qiagen) supplemented with 1% v/v 2-mercaptoethanol. After performing an 2.2x clean up with Agencourt RNAClean XP Beads (Beckman Coulter), whole transcriptome-amplified cDNA products were generated using the SMARTer Ultra-low RNA Kit (Clontech), and conventional Illumina libraries were made and sequenced to an average depth of 27 million read pairs (HiSeq 2000, Illumina). Expression levels and splicing ratios were quantified using RSEM[14] and MISO[18], respectively. Additional experiments were performed using RNA-FISH (Panomics), Immunofluorescence, FACS, and single-cell qRT-PCR (Single Cell-to-CT (Invitrogen) and BioMark (Fludigm)). Full Methods and any associated references are provided in SI.