Research ArticleImmunology

Signals trigger state-specific transcriptional programs to support diversity and homeostasis in immune cells

See allHide authors and affiliations

Science Signaling  14 May 2019:
Vol. 12, Issue 581, eaao5820
DOI: 10.1126/scisignal.aao5820

Defining macrophage states

Macrophages are immune cells that play a major role in maintaining tissue homeostasis. Because of their high plasticity, macrophages can rapidly respond to various external and internal cues, such as exposure to infectious agents or metabolic stress. Through single-cell analyses of mRNA and proteins in macrophage populations and computational analysis, Fischer et al. characterized the differential transcriptional programs driven by the microbial product LPS and the free fatty acid palmitate, both of which engage the same cell surface receptor, TLR4. These data helped to characterize the pro- and anti-inflammatory macrophage states induced by both stimuli, determine how antagonistic genes interacted with each other, and measure transcriptional signaling efficiency in these cells.

Abstract

Macrophages play key roles in the immune systems of humans and other mammals. Here, we performed single-cell analyses of the mRNAs and proteins of human macrophages to compare their responses to the signaling molecules lipopolysaccharide (LPS), a component of Gram-negative bacteria, and palmitate (PAL), a free fatty acid. We found that, although both molecules signal through the cell surface protein Toll-like receptor 4 (TLR4), they stimulated the expression of different genes, resulting in specific pro- and anti-inflammatory cellular states for each signal. The effects of the glucocorticoid receptor, which antagonizes LPS signaling, and cyclic AMP–dependent transcription factor 3, which inhibits PAL-induced inflammation, on inflammatory response seemed largely determined by digital on-off events. Furthermore, the quantification of transcriptional variance and signaling entropy enabled the identification of cell state–specific deregulated molecular pathways. These data suggest that the preservation of signaling in distinct cells might confer diversity on macrophage populations essential to maintaining major cellular functions.

INTRODUCTION

Macrophages are crucial components of the innate and adaptive immune systems of humans and other animals, and they play lifelong roles in maintaining homeostasis in tissues and cell populations. An essential feature of macrophages is their high plasticity (their adaptability to changes in their environments). This is instrumental in their diverse functions in immunity and homeostasis and their ability to react to various external and intrinsic challenges. Humans are continually confronted with conditions for which they have not been primed by evolution because of the continual rise of new infectious agents, unpredictable environmental conditions (such as metabolic stress), and other events.

Previous studies have identified a number of mechanisms that underlie various functions fulfilled by mammalian macrophages (1). For example, CRISPR-knockout screen analyses of mouse macrophages have revealed that, during acute immune responses, macrophages are efficiently activated to enter a proinflammatory M1 state (2). Such effects can be triggered by signaling molecules such as lipopolysaccharide (LPS) presented by Gram-negative bacteria during infections (2). A different response is required when macrophages cope with excessive amounts of free fatty acids (FFAs), such as palmitate (PAL), which cause molecular stress. This process likely becomes more crucial later in life because PAL and other FFAs substantially increase the risk for metabolic and cardiovascular diseases over time through low-grade, inflammatory processes (3, 4).

Although both LPS and PAL induce inflammatory signaling through the cell surface protein Toll-like receptor 4 (TLR4), they rely on different mechanisms to do so. In the case of LPS, TLR4 activation is stimulated in concert with lymphocyte antigen 96 (Ly96, also referred to as MD-2) and CD14, leading to the intracellular recruitment of the adaptor protein myeloid differentiation primary response 88 (MYD88) and interleukin-1 (IL-1) receptor–associated kinases (IRAKs). This effect induces the downstream activation of pathways mediated by the transcriptional regulators nuclear factor κB (NF-κB), activator protein 1 (AP-1), and interferon regulatory factor 3 (IRF3), resulting in substantial expression of genes encoding proinflammatory factors, such as cytokines, chemokines, and other effectors of the innate immune response (5). In contrast, the activation of TLR4 by PAL through the ligand fetuin A leads to a low-grade inflammatory response and the induction of various metabolic pathways caused by FFA-induced stress (6).

In general, many studies of cellular responses have relied on the analysis of bulk cell populations, potentially losing crucial information regarding individual cells. However, the diversity and flexibility of cells are crucial in coping with environmental challenges and ensuring the continued homeostasis of physiology. It has been unclear, for example, whether populations of macrophages stimulated by LPS or PAL are all competent to respond in a similar manner to these different TLR4 effector molecules.

Miniaturization has tremendously advanced genome research over the past two decades (7). Emerging unbiased, single-cell methodologies are becoming increasingly popular as a means for detecting molecular heterogeneity within populations of tens of thousands of cells; microdroplet-based, single-cell RNA sequencing (scRNA-seq) (810) and other omics technologies (11) will enable a far more exact determination of (rare) cell types. Nevertheless, it seems that the detection of cell states featuring only a small number of specific genes remains difficult for these microdroplet-based profiling methods (12); in some cases, the gain of high transcriptome complexity per cell appears more relevant than analyzing a large number of cells (13). For our study, we used a highly sensitive, microfluidics-based, scRNA-seq approach, which, although it is limited in terms of the numbers of cells analyzed, allowed an in-depth characterization of transcriptional profiles. This was especially crucial for analyzing the relative mild cellular response induced by PAL.

Here, we aimed to benefit from single-cell approaches to address a number of questions. In terms of cellular plasticity, how do populations of immune cells, such as macrophages, adapt to different types of physiologically relevant stimuli, namely, LPS, which induces an acute inflammatory response, and PAL, which induces metabolic stress and a low-grade inflammatory response? With regard to gene regulation, can single-cell analyses shed new light on current gene interaction models that are almost exclusively based on conventional bulk experiments? Last, in terms of transcriptional signaling efficiency and coordination, are transcriptional adaptations to “bacterial activation” (LPS) and “metabolic stress” (PAL) equally efficient and can the degree of molecular coordination be quantified to describe and explain cellular deregulation in macrophages? A crucial methodological step toward answering these questions is to eliminate, as far as possible, any potentially uncontrollable or hidden physiological effects beyond the inherent stochastics of the cellular response to signal molecules. To achieve this goal, we focused on an isogenic, differentiated, and cell cycle–arrested human macrophage cell line (THP-1 cells in G0; fig. S1, A and B) and additionally analyzed human monocyte–derived macrophages from healthy donors. In our experiments, we further changed only one parameter to explore and quantify varying cellular responses by treating the cells for either 2 hours with LPS (as an established model for acute inflammatory responses) or 24 hours with PAL (a model for metabolic disease–associated chronic inflammation).

Using well-controlled human macrophage models and LPS and PAL stimulation followed by single-cell analyses, we (i) characterized pro- and anti-inflammatory macrophage states and transcriptional programs induced by both stimulations, (ii) found mutually exclusive, cell state–specific expression of major competing genes such as NR3C1 versus IL1B and ATF3 verus IL1B, and (iii) quantified the signaling efficiency of macrophage populations to signals.

RESULTS

Cellular plasticity: Macrophages form distinct transcriptional states after activation

To analyze the cellular plasticity of macrophages, we initially examined the transcriptomes of untreated, LPS-stimulated (100 ng/ml), and PAL-stimulated (250 μM) THP-1 cell–derived macrophages (hereafter referred to as THP-1 macrophages) using microfluidics-based, scRNA-seq. This analysis resulted in the sensitive detection of several thousands of genes per cell for each signal (see fig. S1C and Materials and Methods). Note that for LPS-treated THP-1 macrophages, we found that the inflammatory marker genes IL8 and IL1B were expressed in 84 and 63% of cells, respectively, whereas for PAL-treated THP-1 macrophage, 5 and 8% of cells showed expression of IL8 and IL1B, respectively (fig. S1D). Results derived from these analyses were corroborated by fully independent experiments using targeted fluorescence in situ hybridization (FISH) detection in hundreds of cells (fig. S1E). These data demonstrated the heterogeneity of the cellular response to external stimuli.

On the basis of single-cell transcriptome data, we next identified specific transcriptional states that represented snapshots of dynamic subpopulations of macrophages in an acute phase of inflammation (as induced by LPS) and metabolic stress including low-grade inflammation (as induced by PAL). Initially, we applied independent component analysis (ICA) based on the genes that were differentially expressed as found by RNA-seq of bulk RNA (fig. S2A). The gene sets for ICA dimension reduction were selected to focus on the effects of either LPS- and PAL-dependent stimulation of THP-1 macrophages. This method enabled a comparative analysis of signal-specific genes that would otherwise be masked by transcriptional variability related to housekeeping functions. The ICA showed that LPS- and PAL-treated macrophages both exhibited two major states after activation: an proinflammatory M1-like state and an anti-inflammatory M2-like state. Furthermore, clustering exposed an additional low-response cell state with a trend toward anti-inflammatory characteristics in case of LPS treatment [here defined as M2-like (low)] and a low-response state that trended toward proinflammatory cell profiles for PAL treatment [here defined as M1-like (low)] (Fig. 1, A and B). In summary, a maximum of three specific LPS-stimulated and three specific PAL-stimulated activation states were found (Fig. 1C).

Fig. 1 ICA-based separation of signal-induced cell transcriptional states.

(A) ICA of scRNA-seq gene expression data. Colors indicate identified cell states of LPS-treated [n = 83 cells; LPS, 22 hours phorbol 12-myristate 13-acetate (PMA) + 2 hours LPS + PMA], PAL-treated (n = 74 cells; PAL, 24 hours PAL + PMA), and untreated (UN) (n = 72 cells; UN, 24 hours PMA) THP-1 macrophages. (B) Average expression of treatment-induced differentially expressed genes, used as an ICA input, in identified cell states (n = 1001 genes for LPS treatment and n = 266 for PAL treatment). Statistical analysis was performed by one-way analysis of variance (ANOVA) with Tukey’s honestly significant difference (HSD) test. **P < 0.01; ***P < 0.001. n.s., not significant; FPKM, fragments per kilobase of per million reads. (C) Definition of macrophage states after activation with LPS or PAL.

Weighted gene coexpression network analysis (WGCNA) is a method to describe the correlation clusters of genes across gene expression datasets. WGCNA is particularly useful in finding modules of highly correlated genes (14). The expression of general inflammatory markers, such as IL1B, and anti-inflammatory markers, such as NR3C1 (LPS) and ATF3 (PAL), varied among the different cell states (Fig. 2, A and B). For example, IL1B was strongly expressed in LPS-treated cells of the M1-like state, whereas ATF3 was particularly expressed in PAL-treated cells of the M2-like state. WGCNA of LPS- and PAL-induced genes revealed that, in both LPS- and PAL-specific M1-like cells, compared to untreated cells, similar pathways were, in part, overrepresented (enriched), for example, pathways for Toll-like receptor (TLR) and peroxisome proliferator–activated receptor (PPAR) signaling. Nevertheless, distinctions could be made between the underlying hub genes and networks observed in the major M1- and M2-like cell states [Fig. 2C, figs. S2B and S3 (A and B), and table S1]. For example, stimulation by PAL induced lipid and lipoprotein metabolism pathways, including known transcriptional regulators of cholesterol homeostasis, such as ATF3 (15), which seemed less relevant during LPS-dependent activation (Fig. 2, C and D, and fig. S3B).

Fig. 2 ICA-based representation of individual genes and pathways in single cells and identified cell states.

(A) Expression of key pro- and anti-inflammatory genes in LPS- and PAL-induced cell states. Each dot represents the expression value in an individual cell from the corresponding cell state. (B) Expression of selected pro- and anti-inflammatory genes in individual cells in ICA-based two-dimensional space. Each point corresponds to an individual cell. (C) Differentially expressed key regulatory pathways in identified cell states. (D) ICA-based representation of key differentially expressed regulatory pathways. Each dot represents an average expression of the genes from the corresponding pathway in individual cells from the specified cell state. Statistical analysis for (A) and (C) was performed by one-way ANOVA with Tukey’s HSD test for LPS-treated (n = 83 cells), PAL-treated (n = 74 cells), and untreated THP-1 macrophages (n = 72 cells). *P < 0.05; **P < 0.01; ***P < 0.001; n.s., not significant.

Next, we aimed to complement our earlier analysis based on differentially expressed genes using the entire single-cell RNA expression datasets. To make use of all of the detected genes, we applied dimension reduction using self-organizing maps (SOMs) (16). The extent of expression of up to ~5000 expressed genes was thereby projected to 100 SOM components (fig. S4A). Single-cell SOMs (which provide a simplified, two-dimensional representation of the transcriptome) revealed that only a fraction of the ~5000 expressed genes appeared to be present in proximal cell groups, suggesting potential co-regulation of these genes (fig. S4B). Consistent with the results from our earlier single-cell analysis of differentially expressed genes, the hierarchical clustering of SOM component correlation matrices provided additional evidence that THP-1 macrophages could be clustered into two main states (fig. S4C): a proinflammatory M1-like state (characterized by the high expression of genes, such as IL1B, together with an enrichment of PPAR signaling and HIF-1α signaling pathways) and an anti-inflammatory M2-like state [characterized by the high expression of genes, such as the glucocorticoid receptor (GR) (NR3C1), and an enrichment of regulatory interleukin signaling] (fig. S5, A to C, and table S2). Note that, in particular, there were substantial differences between the hub genes inferred for the M2-like states of LPS- and PAL-treated cells (fig. S5D), indicating specific transcriptional network responses for each inflammatory stimulus.

Using the two computational approaches described here, we found LPS- and PAL-specific states (comprising M1- and M2-like states), which are consistent with those defined in other studies (17). A physiological situation that is less defined, in which multiple molecules act simultaneously, might reveal a number of substates (18). A deeper look at specific genes of interest in individual cells can reveal new insights. For example, it was shown that LPS can engage aerobic glycolysis through the expression IL10 to generate adenosine triphosphate independently of mitochondrial oxidative phosphorylation (OXPHOS) in bulk macrophages (19). Here, we detected IL10 almost exclusively in some cells of the M1-like state of LPS-activated cells, whereas we found expression of IL10-induced anti-inflammatory factors, such as STAT3 (signal transducer and activator of transcription 3), in almost all macrophage cells irrespective of their states (fig. S6A). Note that the IL10-expressing cells were mostly devoid of STAT3 mRNA, corroborating a study that showed a paracrine mechanism of gene induction to inhibit inflammatory mediators, such as IL-1β, in surrounding macrophages (fig. S6B) (19). The expression of genes associated with metabolic pathways (glycolysis, tricarboxylic acid cycle, and OXPHOS) was greater in the anti-inflammatory M2-like state compared to the proinflammatory M1-like state (fig. S6C). This suggests that a reanalysis of previous studies at the level of single cells might provide a deeper mechanistic understanding of their states.

A large number of genes that were strongly expressed in the proinflammatory cell state in response to LPS (fig. S7A) were repressed in THP-1 macrophages with stable knockdown of Myd88 (THP1-XBlue–defMyD88 cell line), as indicated by an average log 2 fold change of −2 for, on average, 70 genes per cell (fig. S7, A and B, and table S3). In contrast, the knockdown affected far fewer genes (on average, 25) that were specific to the M2-like state of macrophages (average log 2 fold change of −1.2). These data suggest that the loss of intracellular MYD88 signaling results in a substantial depletion of macrophages with proinflammatory M1-like profile (20). Furthermore, we observed that macrophages expressing proinflammatory M1-like genes exhibited a larger, flatter morphology than that of cells with anti-inflammatory M2-like profiles, which were smaller and rounder (Fig. 3, A and B). Treatment with dexamethasone, an anti-inflammatory GR activator, induced the transition of larger cells toward rounder phenotypes (Fig. 3, C and D), and cells exhibiting impairments in their proinflammatory capacity (Myd88-deficient macrophages) exhibited a rounder morphology. Note that even untreated cells showed both proinflammatory M1-like and anti-inflammatory M2-like morphologies (Fig. 3E). In addition, a microscopy analysis over more than 15 hours revealed fluctuations of M1- and M2-like (pre-) states of the same untreated cells over time (Fig. 3E). In summary, to answer our first question, deep scRNA-seq enabled us to define how populations of immune cells, such as macrophages, adapt to different stimuli, such as LPS and PAL, by generating distinct pro- and anti-inflammatory transcriptional (and corresponding morphological) states. Our next effort aimed at leveraging these rich data to obtain insights into the mechanistic cellular interplay of inflammatory mediators with major anti-inflammatory genes.

Fig. 3 Cell morphology can be linked to transcriptional states.

(A) Representative images from RNA FISH analysis of THP-1 macrophages treated with LPS (100 ng/ml) for IL1B (red) and NR3C1 (green) transcripts, as described in Materials and Methods. Bottom: Merged image of fluorescence channels and differential interference contrast images. (B) Quantification of cell size (arbitrary units) and eccentricity (0 = circle, 1 = ellipse) for cells with high expression of the indicated genes, as described in Materials and Methods. N indicates the number of cells analyzed. Data were acquired from two independent experiments. The red box plots represent data from IL1B-positive cells. White box plots represent data from cells with high expression (top 50%) of the genes indicated at the top of the panel. Statistical analysis was done by one-way ANOVA followed by Dunn’s multiple comparisons test. **P < 0.01; ***P < 0.001. a.u., arbitrary units; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; n.s., not significant. (C) Bright-field microscopy images were acquired at 1-hour intervals at constant positions. THP-1 macrophages were treated with LPS (100 ng/ml) alone or together with 1 μM dexamethasone (Dex). Cells in boxes at the 36-hour time point are shown in magnified view. Experiments were repeated three times. Representative images shown here were taken from one experiment. (D) Morphological properties of untreated THP-1 macrophages that had been differentiated for 72 hours and cells that were treated for 48 hours with LPS alone or in the presence of dexamethasone and differentiated THP1-XBlue–defMyD cells. Experiments were repeated three times. Representative images are shown. (E) Representative bright-field, live cell imaging results for cells with changing morphological properties. Image acquisition was performed every 2 min in Z-stack mode for a time period of 16 hours. Individual cells were tracked manually. Every line shows the same cells recorded at different times. Red color bars illustrate (nonquantitatively) big/flat M1-like cells, whereas blue color bars illustrate small/round M2-like cells. Experiments were repeated three times. Representative image regions are shown.

Gene regulation: Mutually exclusive expression of competing gene regulatory factors

In general, previous studies to elucidate fundamental biological mechanisms such as molecular antagonism (the partial or complete suppression of the effect or activity of the other) have been based on the assumption of largely homogenous populations of cells. Antagonistic interactions between molecules expressed at changing quantities within the same cells seemed to govern the switches. However, targeted analyses suggest that a mutually exclusive cellular expression of antagonistic genes or proteins might play a so-far underestimated role (21).

Thus, referring to our second main question, we aimed to investigate gene interactions by making use of single-cell instead of conventional bulk datasets. Note that our single-cell sequencing analyses indicated that genes that were previously thought to operate antagonistically within the same cell were expressed to a large degree in cells of different states. Consistent with the sequencing data, we confirmed the state-specific expression of major antagonistic transcriptional regulator genes, such as NR3C1 versus IL1B and ATF3 versus IL1B. For example, independent single-cell quantitative real-time polymerase chain reaction (sc-qPCR) experiments showed for LPS-treated macrophages a trend of mutually exclusive mRNA expression between the production of inflammatory IL1B in one subset of cells (M1-like) versus the anti-inflammatory genes NR3C1, JAK2, and IRAK3 in another [M2-like (low)] (Fig. 4A and fig. S8A). In contrast, the expression of IL1B was coordinated with that of other known inflammatory factors within the same cells, including HIF1A or IL8 (in the M1-like state) (fig. S8A). The anti-inflammatory genes NR3C1 and IRAK3 were also found to be coexpressed in respective cell states. Alternatively, using RNA FISH detection in independent experiments in 956 cells, we confirmed a substantial pattern of state-exclusive expression (fig. S8B), in which proinflammatory IL1B expression frequently occurred with that of proinflammatory IL8 but rarely with the expression of anti-inflammatory genes, such as NR3C1, IRAK3, JAK2, or FKBP5. These results indicate a mutually exclusive pattern of cellular expression of antagonistic pro- and anti-inflammatory genes. Furthermore, sc-qPCR data revealed that ATF3, an antagonist of metabolically induced chronic inflammation (15), was expressed in M2-like cells that were devoid of inflammatory mediators, such as IL-1β, suggesting that proinflammatory genes were efficiently repressed within these cells (Fig. 4A).

Fig. 4 The key antagonistic transcriptional regulator genes NR3C1/IL1B and ATF3/IL1B and their proteins are expressed in distinct transient cell states of LPS- and PAL-stimulated macrophages.

(A) scRNA-seq and sc-qPCR analyses of THP-1 macrophages indicate trends for mutually exclusive expression of the LPS-induced genes NR3C1 and IL1B and the PAL-induced genes ATF3 and IL1B. For the scRNA-seq experiments, n = 83 cells for LPS treatment and n = 74 cells for PAL treatment; for the sc-qPCR experiments, n = 88 cells for LPS treatment and n = 168 cells for PAL treatment. (B) Independent immunofluorescence detection of antagonistic proteins GR and IL-1β (for LPS treatment) and ATF3 and IL-1β (for PAL treatment) compared to the highly correlating protein pairs IL-8 and IL-1β (for LPS) and PPARG and IL-1β (for PAL) in primary macrophages. Each image represents an individual specimen for the detection of the indicated pair of proteins. Alexa Fluor 488– or Alexa Fluor 594–conjugated secondary antibodies were used for protein detection. Nuclei were detected with 4′,6-diamidino-2-phenylindole (DAPI) (blue staining, shown nuclei or nuclei borders). Cell borders defined by segmentation are shown in pink. Scale bar, 20 μm. (C) Quantification of immunofluorescence staining as shown in (B). N indicates the number of cells analyzed. R represents the Spearman’s rank correlation coefficient in (A) to (C).

Deeming it important to determine whether the trend for the mutually exclusive cellular gene expression was also reflected at the level of protein, we used immunofluorescence double staining to perform single-cell analyses of proteins in isogenic and primary human macrophages (Fig. 4, B and C, and fig. S8C). The results were largely consistent with the mRNA expression data we had captured from single cells: Our protein analyses revealed, to a large extent, cell state–specific, mutually exclusive expression patterns of proteins that belonged to central antagonistic pairs, such as anti-inflammatory GR and proinflammatory IL-1β, or anti-inflammatory activating transcription factor 3 (ATF3) and proinflammatory IL-1β.

If external signal molecules were acting as stimuli in prompting distinct bimodal patterns of responses on the part of individual cells, we wondered how increasing the concentration of these signal molecules would influence the percentage of cells in a particular state (indicated by state-specific genes) in terms of the expression of functionally relevant transcripts within those cells and what effects this would have on the population as a whole. Using sc-qPCR, we found that the number of macrophages expressing the proinflammatory genes IL1B and HIF1A was statistically significantly increased with increasing concentrations of LPS (0 to 1000 nM) (Fig. 5A, left), whereas the abundance of transcripts of these genes within each cell remained rather constant in response to concentrations of LPS between 10 and 1000 nM (Fig. 5A, right). These data indicate that stimulation with low (10 nM) concentrations of LPS might suffice to stimulate the maximal expression of key mRNAs within single cells. Note that as concentrations of LPS increased, we simultaneously detected an increase in the number of cells expressing the anti-inflammatory gene NR3C1 but without substantially changing the quantity of transcripts of the gene within individual cells (Fig. 5A). These effects seem counterintuitive under models based on bulk cell population experiments. Furthermore, the single-cell data obtained in this study suggest that a balance was maintained between the transcriptional responses of pro- and anti-inflammatory cells, which offers a means by which homeostasis can be maintained at the level of the entire cell population.

Fig. 5 Changes in cell state composition.

(A) Top: Increasing concentrations of LPS resulted in increases in the percentage of cells (THP-1 macrophages) that showed increased gene expression (right), whereas the increase in the number of transcripts per cell was largely constant (right). Gene expression was assessed using sc-qPCR for 172 cells treated with different concentrations of LPS. The percentage of (responsive) cells was calculated for every gene under investigation. All cells that showed expression of at least one gene of interest constituted the total quantity of cells (positive or responsive cells). For all positive cells, the numbers of transcripts per cell were calculated by absolute quantification relative to a DNA standard of known concentration and estimated molecules per volume. Data were derived from two independent experiments. **P < 0.01; ***P < 0.001; otherwise not significant. Statistical analysis was performed by two-way ANOVA with Bonferroni posttest; mean values (bars) and SDs (error bars) are indicated. Bottom: Density plots show the distribution of gene expression (transcripts per cell) as detected by sc-qPCR. Data were derived from two independent experiments. (B) Model for the inhibitory interaction of transcriptional regulators NR3C1 and ATF3. Multiple studies have shown the (competitive) inhibitory effects of GR and ATF3 on production of inflammatory mediators, such as IL-1β (51, 52). Our analyses of single cells indicated that macrophages expressing genes encoding GR or ATF3 were almost devoid of IL1B expression and of IL-1β protein, whereas cells that did not contain GR or ATF3 showed strong expression of IL1B.

In summary, our single-cell data indicate how antagonistic genes interact with each other. For example, the expression of transcriptional inhibitors such as NR3C1 (which encodes GR that combats acute inflammation) (22) or ATF3 (whose gene product acts against chronic metabolically induced inflammation) (15) can lead to an almost complete depletion of inflammatory mediators, such as the marker proteins IL-1β and IL-8, within a single macrophage (Fig. 4, A and B). These results suggest that single cells experience, to a large degree, digital “on” or “off” (0 or 1) switches of the expression of key inflammatory genes, committing macrophages to a transient pro- or rather anti-inflammatory state (Fig. 5B). One general advantage that may be achieved by this type of intercellular cooperation may involve an economization of intracellular resources (23). Thus, increasing the amounts of ligands (to greater than the minimal amount in our defined cell culture conditions) can enhance the transcriptional effects of nuclear receptors, such as GR.

Signaling efficiency and coordination: Quantification of the transcriptional cellular response to signal molecules

Last, single-cell sequencing data enabled us to answer our third main question, to quantify transcriptional signaling efficiency and coordination at the cell system level. This enabled us to gain insights into the adaptability of the macrophages to varying external stimuli (their plasticity, as defined earlier).

First, we explored the variance of gene expression between individual cells in response to changes in external stimuli by using two quantification methods: Shannon entropy evaluating uniformity of the expression profiles and the distance-to-median (DM) measure that was designed to be independent of gene or pathway expression abundance (see Materials and Methods). Expression variability was largely homogeneous among the states within each treatment (fig. S9). Furthermore, Shannon entropy analysis of scRNA-seq indicated high reproducibility, which was confirmed by additional independent sc-qPCR experiments (fig. S10A). Next, we estimated cell-to-cell variability using the DM measurement (see fig. S10, B and C, and Materials and Methods). In terms of the whole-transcriptome analysis of LPS- and PAL-specific states, the highest gene expression variance was observed in low-responsive, M2-like macrophages activated by LPS and M1-like macrophages activated by PAL (fig. S11A). In contrast, cells in M2-like states exhibited a lower amount of variance in gene expression. For signaling pathways induced by LPS or PAL, we also found the highest heterogeneity in the low-responsive, M2-like state after LPS treatment and in the M1-like state after PAL treatment, whereas the lowest heterogeneity was found in the M2-like states (fig. S11B). This was particularly true for pathways related to inflammation. To summarize, analysis of gene expression variance revealed cell state–specific heterogeneity of the expression of various sets of functionally related genes, indicating higher heterogeneity in the proinflammatory versus the anti-inflammatory states.

Next, by computing signaling entropy (SE) (2426), we quantified the transcriptional efficiency of macrophages to establish signal transduction networks. In other words, we analyzed the relative degree of molecular coordination of transcriptional responses of macrophages to the signals LPS and PAL. SE, a concept introduced by Teschendorff (25, 26), represents an overall measure of molecular pathway promiscuity by integrating single-cell gene expression data with protein-protein interaction (PPI) networks (Fig. 6A).

Fig. 6 Analysis of SE reveals the organization of intracellular transcriptional responses to the signal molecules LPS and PAL.

(A) SE network model. High SE is suggestive of high uncertainty, thus promiscuous, inefficient usage of network components, whereas lower SE indicates lower uncertainty and a more deterministic, coordinated usage of functional networks of cells. (B) Normalized global transcriptional SE for macrophage states under the indicated conditions for LPS-treated (n = 83 cells), PAL-treated (n = 74 cells), and untreated THP-1 macrophages (n = 72 cells). (C) SE analysis for the TNF-α and NF-κB signaling pathway (n = 519 genes). Statistical analysis for B and C was performed by one-way ANOVA with Tukey’s HSD test. *P < 0.05; **P < 0.01; ***P < 0.001; n.s., not significant. (D) Summary of the defined macrophage states and associated entropy states after treatment with PAL or LPS.

Analyzing the genes in terms of signals above background expression noise, the lowest SE was found for LPS treatment, a medium SE was found for PAL treatment, and the highest SE was found for unstimulated macrophages (Fig. 6B), indicating more conformity of cellular responses for genes that were highly expressed during activation. However, with regard to macrophage states, SE was rather high for the anti-inflammatory M2-like states in which gene expression variance was comparatively low, whereas the proinflammatory M1-like states showed lower SE but high gene expression variance (Fig. 6B and fig. S11A). These results indicate that cells of the M1-like state may establish a more coordinated transcriptional response of functionally linked genes than do cells of the M2-like state but that this is accompanied by a more heterogeneous enhancement of the expression of these genes.

Note that, at the level of key inflammatory pathways, such as those mediate by tumor necrosis factor–α (TNF-α) and NF-κB signaling, SE was decreased for the PAL-stimulated, M2-like state that also showed low gene expression variance (Fig. 6C). This result might suggest coordinated action across this metabolically induced, M2-like state in terms of usage and gene expression of functionally related genes, corroborating earlier mechanistic studies of low-grade inflammation processes in metabolic diseases, such as diabetes and atherosclerosis (27). In summary, the analysis of SE revealed varying degrees of cellular determinism to react to different stimuli, such as LPS or PAL (Fig. 6D). These results quantitatively indicate different usage of functionally relevant genes in the various transcriptional states. As shown here, this holistic cell system parameter analysis was only manageable through the use of unbiased single-cell transcriptome data, an emerging approach in the analysis of cell signaling.

DISCUSSION

In this study, we aimed to answer three main questions related to the adaption of macrophages to key immunological and metabolic triggers to gain insights on (i) how populations of immune cells, such as macrophages, adapt to different types of physiologically relevant stimuli (cellular plasticity), (ii) how antagonistic genes interact with each other (gene regulation), and (iii) transcriptional signaling efficiency and coordination. Through a single-cell approach, we analyzed gene expression in macrophage populations that were representative of either an acute inflammatory response (2 hours of LPS treatment) or metabolic stress (24 hours of PAL treatment). Answering our questions required well-controlled laboratory conditions and, particularly, resolution at the level of single cells, which could be achieved by applying particular single-cell genomics approaches. We endeavored to control conditions by changing only one parameter in vitro (LPS or PAL). However, local effects, such as proximity of cells to each other, can lead to interference, for example, through paracrine signaling (28). Nevertheless, in vivo applications would be much more difficult to monitor, leading to numerous additional partly hidden factors that might interfere with LPS or PAL. Note that macrophages derived from mice or other animal models are at least, in part, differentially regulated compared to those of humans. For example, human macrophages are endowed with the autocatalytic increase in the abundance of the nuclear receptor liver X receptor α (LXRα) to adapt to metabolic stress. This key pathway seems absent in mouse macrophages (2931).

It would be desirable to gain additional dynamic information in the future by analyzing single-cell transcriptomes for various time points for LPS and PAL treatment. High-resolution microscopy imaging of single cells including sophisticated computational tools (currently available to few laboratories worldwide) permits the tracking and analysis of few labeled biomolecules over time; this approach can only be carried out in vitro but not in vivo (32). Live cell imaging approaches and subsequent microfluidic-based single-cell sequencing technologies have been integrated in the detection of both the dynamics of LPS-dependent activation of the innate immune transcription factor NF-κB and the underlying global transcriptional response in the same single cell (33). This approach enabled distinct cytokine expression patterns to be correlated with the activation dynamics of NF-κB to explain (in part) cellular phenotypes.

Using the systems approaches described here, including independent scRNA-seq, sc-qPCR, FISH, and protein immunofluorescence detection experiments, we could decipher the adaptability of macrophages to LPS and PAL. The term “state” refers to analytical snapshots of transient macrophage subpopulations defined by specific transcriptional cellular profiles. Here, we first aimed to analyze the plasticity of macrophages in terms of their transcriptional adaptation to different types of physiologically relevant stimuli (LPS and PAL). After 2 hours of LPS treatment (adaptation to acute inflammation), we observed main M1- and M2-like states, as well as a minor M2-like low state (Fig. 1C). After 24 hours of PAL treatment (adaptation to metabolic stress), we found mainly M1- and M2-like states, as well as a minor M1-like low state (Fig. 1C). Compared to the LPS-induced states, the PAL-induced states showed a lower extent of gene expression, together with a reduced number of expressed genes and varying relevant gene sets (Fig. 2, A to C). Intracellular MYD88 signaling was found to be crucial in establishing the M1-like response (fig. S7B). Macrophages with M1-like gene expression profiles exhibited a larger, flatter morphology than the rounder and smaller M2-like cells (Fig. 3, A and B). Even untreated cells showed both proinflammatory M1-like and anti-inflammatory M2-like morphologies, which fluctuated transiently more than 15 hours between both states (Fig. 3E). Time-resolved detection using the earlier mentioned sophisticated imaging and cell tracking equipment might enable a more efficient analysis of the potential predetermination of cells to M1- or M2-like states at the point of time of treatment with LPS or PAL.

In particular, the rather weak transcriptional response of macrophages to PAL compared to the response to LPS rendered the analysis of cellular plasticity difficult. The detection of cell states relied on a comparably low number of treatment-induced genes with relatively small expression changes compared to those of unstimulated cells. Gaining high transcriptome complexity per cell became important; this could be achieved by applying microfluidics-based RNA-seq with high-sequencing depths compared to microdroplet-based procedures (12). In general, the results derived from our ICA approach and alternative machine learning approaches, such as SOMs, resulted in a largely congruent assignment of cellular states, defined by state-relevant genes and gene networks, thereby rendering the conclusions potentially more robust than relying on only one computational approach.

The conceptual approach of this study, namely, single-cell analyses of molecular perturbation of cells by varying compounds, may provide a way of refining the Connectivity Map Project (www.broadinstitute.org/connectivity-map-cmap). This and other major initiatives aim to generate a catalog of (bulk) transcriptional signatures representing systematic in vitro cellular perturbations with pharmacologic compounds to identify molecular connections. As shown here, such bulk analyses can blur underlying cellular heterogeneity, suggesting that single-cell transcriptional signatures might add valuable information to the cellular effects of compounds such as anti-inflammatory drugs. In summary, our results derived from single-cell analyses comprehensively describe cellular adaptation to external stimuli and provide a resource for developing hypotheses for mechanistic studies and for quantitative systems analyses.

Single-cell data enabled us to shed light on antagonistic gene-gene interactions of major gene regulatory factors of the inflammatory response that have so far gone largely unexplored. As exemplified by the major acute inflammation (LPS) inhibitor, the nuclear receptor GR, or by an inhibitor of metabolically induced chronic inflammation (PAL), the transcriptional regulator ATF3, we could observe almost digital on-off events. Our results, on the basis of scRNA-seq and independent targeted analyses of mRNAs and cellular proteins (Fig. 4B and fig. S8C), indicate that the molecular responses of macrophages seem, to a large degree, to be regulated in different transient transcriptional subpopulations (states). Our datasets show that single-cell analyses can be of high value for making valid hypotheses and performing proper mechanistic experiments in an adequate cellular context, corroborating other findings (26). In particular, unbiased, single-cell analyses can be used to avoid falsely linking genes that represent different cellular states and which may thus not directly interfere with each other (Fig. 5B). For example, ATF3 was identified as a metabolic target gene in macrophages that reduces the expression of TLR-induced proinflammatory cytokines. However, mechanistic insights were gained by bulk experiments, thereby suggesting competitive inhibition of ATF3 and inflammatory mediators, such as IL-1β. As shown here, single-cell analyses indicate strong inhibition of inflammatory response in ATF3-positive cells, whereas inflammatory markers, such as IL-1β, were rather strongly expressed in ATF3-negative cells. Similar trends for cellular on-off events were observed with GR for acute LPS-induced inflammation. Notably, the development of (quasi) cellular states may lead to transient incidents of intercellular cooperation, thereby balancing cellular homeostasis and economizing intracellular resources (Fig. 5, A and B) (20).

In general, cellular plasticity relies on the resilience of molecular networks rather than on a few single “targets.” Therefore, our third main question referred to analyzing transcriptional signaling efficiency and coordination of cells depending on varying stimuli. Concretely, we aimed to quantify the efficiencies of transcriptional adaptations of macrophages within an entire population to bacterial activation (LPS) and metabolic stress (PAL) and to relatively quantify the underlying degree of molecular coordination between cells. This can only be done by using whole transcriptome data from individual cells. As shown here (fig. S11), individual cells can, in part, react very differently to varying signal molecules, for example, with regard to the functions and overall number of genes involved. Calculation of the variance of gene expression for individual macrophages indicated a cell state–specific heterogeneity of the expression of genes and gene pathways. The highest heterogeneity was found in LPS-induced, proinflammatory M1-like states, suggesting a boost of gene induction compared to all other states defined in this study. The activation of macrophages seemed to be accompanied by higher “molecular disorder,” suggesting that genes with no relevant function exhibited (stochastic) transcriptional background noise.

Next, the determination of SE enabled the (relative) quantitative analyses of cell populations to make use of transcriptional signal transduction network components. In other words, this analysis sheds light on the coordinated usage of known functional networks of cellular components. SE relies on well-validated PPI networks. Using gene expression cutoff values maximizes the reliability of input gene expression data but includes smaller PPI networks. To focus on functionally relevant genes, we excluded background noise, resulting in our consideration of ~1700 expressed genes (with mean read counts per gene greater than 75). Note that for the M2-like states, SE analysis indicated a lower degree of functional connection of interacting molecules. Although potentially intriguing in terms of a higher disorder in M2-like compared to M1-like states, this might also be a result of interactions of genes or proteins relevant for so far less explored M2-like states that have not yet been explored and have not yet been documented in the PPI databases.

The quantitative concepts introduced here to describe transcriptional coordination in adaptations to changing signal molecules are still in their infancy. Single-cell and other approaches will benefit from the further development of validated resources, such as PPI databases. Moreover, using SE and single-cell analyses of many more cell types and perturbations will enable us to put the effects of observed statistically significant changes of SE in a broader biological context. Nevertheless, as was shown previously, similar changes to those shown in this study were observed for normal and tumor tissues (34).

In general, as shown here, using comprehensive transcriptional data derived from the single cell, the fundamental unit of biology, will enable (i) the description of cellular plasticity and the discovery of relevant cellular states, (ii) the identification of thus far hidden mechanisms of gene regulation, and (iii) the quantification of the degree of coordinated adaption of cells to signal molecules. The single-cell approaches applied here offer broader insights into fundamental principles of gene regulation that have been masked in the past by the limitations imposed by studying large cell populations. One eventual outcome should be to better pave the path toward our understanding of cellular signaling but at a cellular resolution. In general, single-cell analyses will shed new light on the paradigm that illnesses originate from the pathology of cells (35).

MATERIALS AND METHODS

Macrophage systems

To decipher the role of varying signal molecules on transcriptional responses, a well-controlled experimental human cellular model is needed to eliminate any other potential uncontrollable or hidden physiological effects beyond inherent stochastics of cellular response. For example, variation of cell cycle states should be kept to a minimum, and homogenous surrounding external conditions should be maximized as much as possible to address our research questions. Note that reprogramming of human cells, such as macrophages, is difficult to model with mouse cells. For example, in contrast to mouse macrophages, human macrophages appear, in part, to use different transcriptional programs (31, 36). Thus, we aimed to benefit from comparative single-cell analyses of individual isogenic, fully differentiated, cell cycle–arrested macrophages (G0) to investigate differential cellular responses to LPS (“acute inflammatory response”) and PAL (metabolic stress). The isogenic THP-1 cell model is considered a prime system for the analysis of human macrophages (37). We further aimed to corroborate potential effects using primary human macrophages.

Cell culture

Isogenic human THP-1 cells [a monocytic-like leukemia cell line; #TIB-202, American Type Culture Collection (ATCC)] and THP1-XBlue–defMyD cells (#thpx-dmyd, InvivoGen) were cultured in RPMI 1640 (#FG1215, Biochrom) supplemented with 10% fetal bovine serum (#S0615, Biochrom). THP-1 cells were seeded from frozen ATCC stocks (5 × 106 cells/ml). Cells were incubated at 37°C in an atmosphere of 5% CO2. Subculturing was done by adding fresh medium to reach a cell concentration of 2 × 105 cells/ml. Medium was completely renewed every week. Human primary macrophages were isolated from at least four individual buffy coats donated by healthy volunteers (provided by DRK Berlin with informed consent for research). Peripheral blood monocytes were isolated from buffy coats with Ficoll-Paque and MACS Monocyte Isolation Kit II with MACS LS columns (Miltenyi Biotec). Primary monocytes were cultured and differentiated to macrophages in RPMI 1640 (#FG1215, Biochrom) supplemented with 10% human serum (Human Serum Type AB, P30-2501, PAN-Biotech).

Cell treatments

THP-1 cells were differentiated to macrophages with 10 nM PMA (#P8139, Sigma-Aldrich) for 48 hours before being treated with either LPS or PAL. Cells to be stimulated with LPS were kept for an additional 21 or 22 hours in culturing medium with 10 nM PMA before being treated with LPS to be consistent with the conditions for PAL-stimulated and unstimulated cells. Primary macrophages were kept in culture for 7 days in the presence of 10% human serum before being treated with LPS or PAL. Cell cycle arrest was assessed by conventional flow cytometry analysis (fig. S1B). For LPS stimulations, primary macrophages and THP-1 cell–derived macrophages were treated with LPS (100 ng/ml; #L5293-2ML, Sigma-Aldrich) combined with 10 nM PMA for 2 hours for RNA analysis or for 3 hours for protein analysis. For PAL stimulation, 62.5 mM sodium palmitate (#P9767, Sigma Aldrich) powder was dissolved at 70°C for 30 min in 0.1 M NaOH. The PAL solution was then coupled with 25% bovine serum albumin (BSA; #A7030, Sigma-Aldrich) to achieve a molar ratio of PAL:BSA of ~7.5:1, which was incubated for 30 min at 37°C and diluted in culture medium with 10 nM PMA to achieve a final PAL concentration of 250 μM. PAL-containing medium was sterilized with a 0.22-μm vacuum filtering system (#SCGPT02RE, Merck Millipore), and the PAL concentration was controlled on the basis of assessment with a Free Fatty Acid Quantification kit (#K612, BioVision). For unstimulated (control) cells after 48 hours of differentiation, the culture medium was replaced with medium containing 10 nM PMA, and the cells were incubated for an additional 24 hours.

RNA isolation and quality control

RNA isolation was performed using the RNeasy Mini Kit (#74106, Qiagen) according to the manufacturer’s protocol with a 30-min on-column deoxyribonuclease digestion step. The RNA concentration was determined using a NanoDrop 2000 (Thermo Fisher Scientific), and additional quality control was performed using a Bioanalyzer (Agilent).

Single-cell mRNA sequencing

To explore the influence of two varying stimulatory signal molecules, LPS and PAL, we performed transcriptome-wide gene expression analysis of single cells. Differentiation of THP-1 cells into macrophages was accompanied by cell cycle arrest (fig. S1B). The cells were then left unstimulated, were treated with LPS (100 ng/ml; near-saturation) for 2 hours, or were treated with 250 μM PAL for 24 hours. Cells were stimulated with LPS (100 ng/ml) for 2 hours to induce an early inflammatory response so that we could detect transcriptional networks of primary rather than secondary response effects (2). We used the Fluidigm C1 Single-Cell Auto Prep System for single-cell preparation (#100-5761/#100-5760; #100-6201, Fluidigm) of 96 cells per treatment group followed by Nextera XT library preparation (#FC-131-1096/#FC-131-1002, Illumina) and Illumina-based RNA-seq (#PE-401-3001, #FC-401-3001, and #PE-121-1003, Illumina). In total, 234 cells were sequenced for this initial experiment [78 ± 10 cells/integrated fluidic circuit (IFC) average ± SEM] with the assumption of observing few major macrophage states. In addition, we generated bulk control samples from 1 × 106 cells and isolated total RNA from three biological replicates for each condition. Library preparation was performed with the conventional TruSeq stranded mRNA approach from Illumina according to the manufacturer’s protocol. All samples were sequenced in paired-end 100-nucleotide mode at the next-generation sequencing core facility of the Max Planck Institute for Molecular Genetics.

Primary analysis of single-cell and bulk sequencing data

The quality of the raw reads was assessed using sequence grooming tools from FastQC and RseQC (38). Reads were mapped to the human reference genome (GRCh37) using STAR version 2.4.0d (see fig. S12 for data quality representation). Each cell was sequenced to an average saturated depth of 5.8 ± 0.24 million mapped read pairs (table S4). A fairly uniform gene coverage was observed without significant 3′/5′-end bias (fig. S12A). Aggregated single-cell expression data agreed with the matching bulk controls with correlations plateauing once we had sampled ~30 to 40 cells (fig. S12, B and C). Gene expression estimation was determined using HTSeq-count (from HTSeq version 0.6.1) (39) for raw read count or Cufflinks (version 2.2.1) (40) for FPKM-normalized expression values, respectively. GENCODE annotation (release 19/GRCh37.p13) was used as a reference. Similar approaches for mapping and read counting were applied for the bulk RNA-seq data. On the basis of microscopic observations, 229 macrophage cells (LPS treatment, n = 83 cells; PAL treatment, n = 74 cells; unstimulated condition, n = 72 cells) were retained for further scRNA-seq analyses (pictures were taken from every microchamber of the Fluidigm C1 IFCs). Sequencing data derived from empty chambers or chambers with cell doublets or deformed or damaged cells were excluded (fig. S13). Data from microchambers showing a mean expression cutoff (derived from log10-transformed count data) less than 0.35, 0.25, and 0.15 for LPS-treated, PAL-treated, and untreated cells, respectively, were excluded. Differential expression analysis of bulk RNA-seq data was performed with DESeq2 software (41) with standard parameters.

ICA of single cells

Two-dimensional ICA was performed with the Sincell R package (42). Log-transformed [log10(FPKM + 1)] expression values in single cells were used as an input for the ICA. Only genes that were statistically significantly increased in expression in bulk RNA-seq (fold change, >2; P < 0.001) and with average expression in single cells greater than FPKM = 1 were used (1001 genes for LPS-treated cells and 266 genes for PAL-treated cells). Cell states were defined on the basis of cell density distribution along the IC 1 and IC 2 axes (fig. S2A).

WGCNA of single-cell data

WGCNA was performed with the WGCNA package in R (14). Specific LPS and PAL networks and topological overlap matrices were constructed on the basis of Pearson’s correlation using the same genes as selected for ICA. The soft threshold parameter was set to 10 for both the LPS and PAL networks on the basis of the criterion of scale-free topology. Hierarchical clustering with average agglomeration method on signed topological overlap matrices was performed, which was followed by cluster definition with the minimum gene cluster size set to 12. We identified 15 and 5 gene coexpression modules in the LPS and PAL networks, respectively. Analysis of overrepresented gene pathways was performed for defined gene modules using the DAVID functional annotation tool version 6.7 (43). The term “pathway overrepresentation” of the DAVID software that is mentioned in Results is synonymous to “enrichment.” DAVID analysis was done with default parameters. Briefly, gene lists (corresponding to modules that were inferred using WGCNA software and assigned to the defined cell subpopulations) were submitted for analysis by DAVID. The enriched (or overrepresented) pathways were defined on the basis of the standard analysis (P < 0.05), which calculates overrepresentation of the pathway gene sets in a given gene list relative to the human population background.

SOM-based dimension reduction

To separate cells based on their detectable transcriptome differences, we performed dimension reduction using Kohonen’s SOMs implemented in R (44) with a grid size of 10 by 10. The SOM algorithm computes two-dimensional representations of each single-cell transcriptome with cell-specific spot patterns whereby each spot includes several adjacent SOM components. Cell-specific overexpression spots are selected among all SOM components using a 98% quantile criterion. If SOM component expression values exceed the 98% quantile level of the expression range of all SOM components, then the individual genes associated with those SOM components define the most stringent genes that characterize a single-cell identity or profile. The dimension-reduced matrices were used to build correlation matrices with rcorr (R, Hmisc package). Statistically significant pairwise correlations were kept (P < 1 × 105) for hierarchical clustering analysis to determine basic cell states. SOM-reduced gene expression matrices were used for ICA using the fastICA R package. In ICA space, cells were annotated according to the results from hierarchical clustering.

SOM-based exploratory single-cell analysis

To identify genes specific for cell states, we performed differential gene expression analysis including computation of the fold change, weighted average difference score, and shrinkage t score according to previously published procedures (45). To quantify responsiveness of cells, we first computed cell-specific marker genes on the basis of SOM portrait analysis (45). Responsiveness was computed as treated versus untreated cell groups. Responsiveness to knockout of MYD88 was computed as treated control cells versus treated MYD88 knockdown cells. We then used differentially expressed genes determined from bulk mRNA-seq data using DESeq2 (41) and associated those to cell-specific genes (defined as cell-specific overexpression spots using a 98% quantile criterion). To identify state-specific major regulatory hubs genes, we defined 10 correlation clusters using treatment-specific average SOM portraits as described previously (45). Correlation cluster–specific genes served as input for FANTOM4 EdgeExpressDB analysis (46).

Quantitative real-time PCR

Single cells and 40-cell controls were sorted into 96-well plates filled with 5 μl of ribonuclease-free water using a BD FACSAria II (BD Biosciences). RNA from bulk samples was extracted and purified as described earlier. For individual cells and the 40-cell controls, the complementary DNA (cDNA) synthesis reaction was performed directly in 5 μl of the cell lysate using qScript cDNA SuperMix (#95048, Quanta BioSciences) in a standard thermal cycler (5, 30, and 5 min at 5°, 42°, and 85°C, respectively). High-performance liquid chromatography–purified custom-designed primers (table S4) were obtained from Sigma-Aldrich. Primer working solutions (2 μM) were freshly prepared. In addition, we used primers provided in RT2 Profiler Toll-Like Receptor Signaling Pathway and Human Glucocorticoid Signaling PCR Array (Qiagen). PerfeCTa SYBR Green SuperMix (#95054, Quanta BioSciences) was used for all qPCRs from bulk samples and single cells. Reactions were prepared in a 10-μl volume in 384-well plates (#PCR-384-LC480-W, Fisher Scientific GmbH). The qPCR run was performed using the LightCycler 480 II System (LC480, Roche). Melting curve analysis was performed to determine the specificity for each sample. To estimate absolute mRNA transcript numbers, an “Interplate Calibrator” standard (c = 106 copies/μl; #IPC250S, TATAA Biocenter) was used. The percentage of cells was calculated for every gene for the various LPS concentrations applied. All cells that showed expression of at least one gene constituted the total quantity of cells on which the percentage calculation was based. For correlation analysis, Cq values were transformed as follows: 2((Cmax− Cq−), with Cq denoting the quantitation cycle and Cmax denoting the number of performed PCR cycles. Afterward, relative expression (Er) values were scaled (by using R’s scale function) and zeroed (by subtracting the minimum value found in the data space) (47). Note that subsequent scaling of the data did not change the structure of the data nor the estimated correlation coefficients. Instead, scaling was used for better visualization (avoiding very small numbers).

RNA FISH

For RNA FISH preparations, 3.25 × 105 cells after treatment were fixed on poly-l-lysine–coated slides with fixation buffer [3.7% formaldehyde (#F8772-25ML, Sigma-Aldrich), 10× phosphate-buffered saline (#P5493-1L, Sigma-Aldrich), and nuclease-free water (#AM9937, Ambion)] for 10 min at room temperature. Fixed cells were washed and permeabilized in 75% ethanol at 4°C overnight. After washing with wash buffer [20× saline sodium citrate (SSC; #S6639-1L, Sigma-Aldrich), 10% formamide (#F9037-100ML, Sigma-Aldrich), and nuclease-free water] for 5 min at room temperature, cells were incubated at 37°C for 4.5 hours with 125 nM RNA FISH fluorescent probes in hybridization buffer [dextran sulfate (100 mg/ml; #D8906-10G, Sigma-Aldrich), 20× SSC, 10% formamide, and nuclease-free water]. The cells were then resuspended in wash buffer for 30 min at room temperature. Last, cells were washed in 2× SSC buffer for 5 min at room temperature. Microscopic slides were prepared with ProLong Gold antifade reagent with DAPI (#P36931, Thermo Fisher Scientific). Images were taken with a Zeiss Axio Observer Z1 wide-field fluorescence microscope equipped with a 60× oil immersion objective and a charge-coupled device camera. Five Z stacks were collected for all images with 0.3-μm spacing. Image analysis was performed using ZEN software (version 2.31). For each sample, 8-by-8 or 10-by-10 tiles were recorded. Binning was set to 3 × 3 pixel. Scaling per pixel was 0.330 μm by 0.330 μm. Bit depth was 14-bit. For background subtraction and normalization of tiles, Gaussian processing with high kernel density settings (kernel density value, 400) was performed to generate a background model for shading correction. Maximum intensity projection of deconvoluted images was used for fluorescence intensity quantification using EBImage (48). Briefly, images from different channels were normalized by subtracting the channel mean intensity from every pixel of the image, and nuclei were segmented and used for propagation to accomplish cell segmentation, as described by Pau et al. (48). After cell intensity–based background subtraction was performed, mean intensities per cell per channel were calculated and transformed [(SDcell − SDac) × Mcell]2, with SDcell denoting the SD of cell intensity, SDac denoting the mean SD of all cells per acquired picture, and Mcell denoting the mean intensity per cell. Cell size and eccentricity were calculated with EBImage as described previously (48). FISH analysis of LPS-treated cells confirmed the results from the scRNA-seq, sc-qPCR, and protein immunofluorescence analyses. Expression changes were less extensive in PAL-treated cells, which made this analysis more difficult. Because FISH detection turned out to be redundant, PAL-treated cells were not analyzed by FISH.

Immunofluorescence double staining

For immunofluorescence analysis, THP-1 cells or primary monocytes were differentiated on poly-l-lysine–coated glass slides in 24-well plates and subsequently treated as previously described with LPS or PAL or were left unstimulated. Fixation and permeabilization were performed on glass slides using the Transcription Factor Buffer Set (BD Pharmingen) according to the recommended manufacturer’s protocol. Blocking was performed with normal goat serum. Specimens were incubated with primary antibodies (mouse and rabbit) at 4°C overnight and subsequently incubated with goat anti-mouse and goat anti-rabbit immunoglobulin G secondary antibodies, conjugated with Alexa Fluor 488 or Alexa Fluor 594 at room temperature for 1 hour. ProLong Gold antifade reagent (#P36931, Thermo Fisher Scientific) with DAPI was used to prepare the specimens. Images were acquired with a Zeiss Axio Observer Z1 wide-field fluorescence microscope. Images consisted of 100 tiles (10 by 10) and were acquired using Z-stack mode (four Z layers) with 3 × 3 binning. Image processing was performed using ZEN 2 software (blue edition, Zeiss) followed by analysis with the EBImage package in R (48). In ZEN 2 software, we processed images using extended depth focus, stitching with fuse tiles and shading correction, and finally exported images in TIFF format for further processing with the EBImage (version 4.20.0) package in R. The EBImage workflow was based on a procedure described previously (48). First, the nucleus signals were processed using adaptive thresholding (function thresh), holes in the objects were filled (function fillHull), and binary segmentation was applied (function bwlabel). Next, cells that were partly captured on the image borders were removed. Then, the signals from two fluorophores were used for cell analysis including thresholding (function thresh), Voronoi-based segmentation to find cell borders given a labeling of the nuclei (propagate), encircling of the defined nuclei and cell borders (function paintObjects), and computing of the features of the defined cells (function computeFeatures). For each cell, the defined features included cell area, perimeter, mean radius, and radius SD, as well as the mean and total intensities of all three fluorophores (DAPI, Alexa Fluor 488, and Alexa Fluor 594). Cell area, perimeter, and radius SD features were scaled by the SD with centering (scale generic function in R) for THP-1 cells and primary macrophages separately. The cells with absolute values of these parameters that were greater than 2 SDs away from the mean were filtered out. Next, for each image, the DAPI mean intensity was scaled by SD with centering (scale function) and also filtered with a 2-SD threshold. Last, the cells with values in the top and bottom 5% (for THP-1 cells) or 10% (for primary macrophages) of the total Alexa Fluor 488 and Alexa Fluor 594 intensity range were removed, and the total intensity was further scaled with centering (using scale function). The cells with total intensity more than 1.5 (for THP-1 cells) or 2 (for primary macrophages) SDs away from the mean were removed. The remaining cells were used for the Spearman’s rank correlation calculation. In the experiment described earlier, we observed that relatively low abundant proteins, such as ATF3 in the ATF3/IL-1β pair, could be detected in a rather small proportion of cells, whereas for highly abundant protein pairs, such as PPARG/IL-1β, we observed a larger number of positive cells. The percentages of cells used for the correlation analysis were 75, 30, 82, and 9.6% for primary macrophages and 7, 10, 41, and 1.7% for THP-1 cell–derived macrophages for IL-1β/IL-8, IL-1β/GR, PPARG/IL-1β, and ATF3/IL-1β, respectively.

Variance estimation in single-cell data

The gene expression variance between cells was estimated similarly to the previously described DM measurement (49). The DM statistic is a measure of the relative variability of each gene, after accounting for the empirical mean-variance relationship. First, mean gene expression and CV2 were estimated for each gene across single cells, and the dependency of CV2 ~1/mean (FPKM) was computed using local polynomial regression (“loess function,” R Bioconductor). Second, the distance between the measured CV2 and the corresponding regression was calculated. The estimated difference represents a variance measure that describes cell-to-cell heterogeneity for each gene and does not depend on mean expression value, which enables direct comparison between different genes and gene groups.

Analysis of Shannon entropy and signaling entropy

Shannon entropy of genes based on gene expression only was calculated as described previously (50) using genes with mean expression >50 transcript counts. Given expression of a gene based on N cells, relative expression of a gene g for a cell c is given as pc|g = wg,c/∑1 ≤ cNwg,c, where wg,c is the expression level of the gene. The entropy of a gene’s expression distribution is Hg = ∑1 ≤ cNpc|g log 2(pc|g). Hg has units of bits and ranges from zero to log 2(N). The normalized entropy rate nH for a gene is calculated as Hc/log 2(N). Signaling entropy was determined according to previously published methods (25, 26) by integration of single-cell gene expression with a global PPI network [Human Protein Reference Database (HPRD) interaction network from Pathway Commons, 13 June 2012]. In the context of a single cell, SE reflects the amount of overall uncertainty in how information is passed on in the molecular interaction network. scRNA-seq read counts (mean read counts per gene, >75) of 1663, 1714, and 1916 expressed genes for untreated, LPS-treated, and PAL-treated macrophages, respectively, were log-transformed, quantile-normalized, and joined to the global HPRD PPI network. Using all genes without setting a detection threshold not only would increase network size but also would increase background noise, resulting in the highest SE for LPS-treated cells, a medium SE for PAL-treated cells, and the lowest SE for untreated cells. Edges of the network were weighted according to the gene expression of neighboring genes that constitute edges. The thereof constructed sample-specific stochastic matrix constituted average local network entropy for genes in the network, as Si = −∑pij log(pij), where pij is the normalized interaction probability of neighboring genes i and j. To facilitate comparison between sample SEs, the network entropy was divided by the maximum attainable value in the given network to derive the normalized local entropy rate.

SUPPLEMENTARY MATERIALS

stke.sciencemag.org/cgi/content/full/12/581/eaao5820/DC1

Fig. S1. Experimental approaches and technical background data.

Fig. S2. ICA of untreated, LPS-treated, and PAL-treated cells and WGCNA gene coexpression module construction.

Fig. S3. Characterization of WGCNA-defined gene coexpression modules and corresponding pathways.

Fig. S4. The SOM-based approach uses whole transcriptome data to determine cell transcriptional states in response to LPS and PAL.

Fig. S5. SOM-based analysis of whole transcriptome single-cell profiles.

Fig. S6. Single-cell gene expression analysis of IL10, STAT3, and IL1B.

Fig. S7. Differential expression analysis of THP-1 cells deficient in MyD88.

Fig. S8. Key antagonistic transcriptional regulators are expressed in distinct transient cell populations.

Fig. S9. Cumulative expression proportion counted for each cell in the different treatments and states.

Fig. S10. Variance estimation of gene expression data derived from single cells: Part I.

Fig. S11. Variance estimation of gene expression data derived from single cells: Part II.

Fig. S12. Quality control of single-cell sequencing data.

Fig. S13. Microfluidic IFC microchamber screening.

Table S1. WGCNA-based analysis of pathways and genes.

Table S2. Output list of SOM-based analysis of pathways and genes.

Table S3. A high number of genes specific for the LPS-induced, proinflammatory state that could be repressed by stable knockdown of MyD88.

Table S4. List of primers.

REFERENCES AND NOTES

Acknowledgments: We thank N. Rajewsky, M. Vingron, R. Zenobi, J. Vogel, and A. J. Ibanez for support and/or scientific discussion. We thank R. Hodge for reading the manuscript. Funding: This work was funded by the European Union [FP7, under grant agreement no. 262055 (ESGI)], the Max Planck Society, and the Helmholtz Association. Author contributions: C.F., M.M., and S.S. designed the experimental and analytical strategy. C.F., M.M., S.B., M.B., and P.G. performed the experiments. C.F., M.M., R.V., S.B., and S.S. analyzed the data. S.S. conceived and supervised the study. M.K. and S.S. managed the project. C.F., M.M., and S.S. wrote the manuscript with input from all other authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. R scripts to calculate the ICA, SOM, WGCNA, variance estimation, and signaling entropy, as well as the PPI network used, can be obtained here at https://github.com/SciGenMDC/Single-cell-macrophage. The sequencing data were deposited in the ArrayExpress database (www.ebi.ac.uk/arrayexpress) under the following accession numbers: E-MTAB-6064 (bulk RNA-seq) and E-MTAB-6075 (single-cell RNA-seq).
View Abstract

Navigate This Article