Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Gene expression variability across cells and species shapes innate immunity

Abstract

As the first line of defence against pathogens, cells mount an innate immune response, which varies widely from cell to cell. The response must be potent but carefully controlled to avoid self-damage. How these constraints have shaped the evolution of innate immunity remains poorly understood. Here we characterize the innate immune response’s transcriptional divergence between species and variability in expression among cells. Using bulk and single-cell transcriptomics in fibroblasts and mononuclear phagocytes from different species, challenged with immune stimuli, we map the architecture of the innate immune response. Transcriptionally diverging genes, including those that encode cytokines and chemokines, vary across cells and have distinct promoter structures. Conversely, genes that are involved in the regulation of this response, such as those that encode transcription factors and kinases, are conserved between species and display low cell-to-cell variability in expression. We suggest that this expression pattern, which is observed across species and conditions, has evolved as a mechanism for fine-tuned regulation to achieve an effective but balanced response.

This is a preview of subscription content, access via your institution

Access options

Rent or buy this article

Prices vary by article type

from$1.95

to$39.95

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: Response divergence across species in innate immune response.
Fig. 2: Transcriptionally divergent genes have unique functions and promoter architectures.
Fig. 3: Cell-to-cell variability in immune response corresponds to response divergence.
Fig. 4: Relationship of response divergence and other evolutionary modes.

Similar content being viewed by others

Data availability

Sequencing data have been deposited in ArrayExpress with the following accessions: E-MTAB-5918, E-MTAB-5919, E-MTAB-5920, E-MTAB-6754, E-MTAB-6773, E-MTAB-5988, E-MTAB-5989, E-MTAB-6831, E-MTAB-6066, E-MTAB-7032, E-MTAB-7037, E-MTAB-7051 and E-MTAB-7052.

References

  1. Borden, E. C. et al. Interferons at age 50: past, current and future impact on biomedicine. Nat. Rev. Drug Discov. 6, 975–990 (2007).

    Article  CAS  Google Scholar 

  2. Iwasaki, A. A virological view of innate immune recognition. Annu. Rev. Microbiol. 66, 177–196 (2012).

    Article  CAS  Google Scholar 

  3. Nielsen, R. et al. A scan for positively selected genes in the genomes of humans and chimpanzees. PLoS Biol. 3, e170 (2005).

    Article  Google Scholar 

  4. Haygood, R., Babbitt, C. C., Fedrigo, O. & Wray, G. A. Contrasts between adaptive coding and noncoding changes during human evolution. Proc. Natl Acad. Sci. USA 107, 7853–7857 (2010).

    Article  ADS  CAS  Google Scholar 

  5. Fumagalli, M. et al. Signatures of environmental genetic adaptation pinpoint pathogens as the main selective pressure through human evolution. PLoS Genet. 7, e1002355 (2011).

    Article  CAS  Google Scholar 

  6. Enard, D., Cai, L., Gwennap, C. & Petrov, D. A. Viruses are a dominant driver of protein adaptation in mammals. eLife 5, e12469 (2016).

    Article  Google Scholar 

  7. Barreiro, L. B. & Quintana-Murci, L. From evolutionary genetics to human immunology: how selection shapes host defence genes. Nat. Rev. Genet. 11, 17–30 (2010).

    Article  CAS  Google Scholar 

  8. Zhao, M., Zhang, J., Phatnani, H., Scheu, S. & Maniatis, T. Stochastic expression of the interferon-β gene. PLoS Biol. 10, e1001249 (2012).

    Article  CAS  Google Scholar 

  9. Avraham, R. et al. Pathogen cell-to-cell variability drives heterogeneity in host immune responses. Cell 162, 1309–1321 (2015).

    Article  CAS  Google Scholar 

  10. Shalek, A. K. et al. Single-cell RNA-seq reveals dynamic paracrine control of cellular variation. Nature 510, 363–369 (2014).

    Article  ADS  CAS  Google Scholar 

  11. Hwang, S. Y. et al. Biphasic RLR-IFN-β response controls the balance between antiviral immunity and cell damage. J. Immunol. 190, 1192–1200 (2013).

    Article  CAS  Google Scholar 

  12. Porritt, R. A. & Hertzog, P. J. Dynamic control of type I IFN signalling by an integrated network of negative regulators. Trends Immunol. 36, 150–160 (2015).

    Article  CAS  Google Scholar 

  13. Ivashkiv, L. B. & Donlin, L. T. Regulation of type I interferon responses. Nat. Rev. Immunol. 14, 36–49 (2014).

    Article  CAS  Google Scholar 

  14. Brinkworth, J. F. & Barreiro, L. B. The contribution of natural selection to present-day susceptibility to chronic inflammatory and autoimmune disease. Curr. Opin. Immunol. 31, 66–78 (2014).

    Article  CAS  Google Scholar 

  15. Kobayashi, K. S. & Flavell, R. A. Shielding the double-edged sword: negative regulation of the innate immune system. J. Leukoc. Biol. 75, 428–433 (2004).

    Article  CAS  Google Scholar 

  16. Kumar, H., Kawai, T. & Akira, S. Pathogen recognition by the innate immune system. Int. Rev. Immunol. 30, 16–34 (2011).

    Article  CAS  Google Scholar 

  17. Barreiro, L. B., Marioni, J. C., Blekhman, R., Stephens, M. & Gilad, Y. Functional comparison of innate immune signaling pathways in primates. PLoS Genet. 6, e1001249 (2010).

    Article  CAS  Google Scholar 

  18. Schroder, K. et al. Conservation and divergence in Toll-like receptor 4-regulated gene expression in primary human versus mouse macrophages. Proc. Natl Acad. Sci. USA 109, E944–E953 (2012).

    Article  CAS  Google Scholar 

  19. Shay, T. et al. Conservation and divergence in the transcriptional programs of the human and mouse immune systems. Proc. Natl Acad. Sci. USA 110, 2946–2951 (2013).

    Article  ADS  CAS  Google Scholar 

  20. Brawand, D. et al. The evolution of gene expression levels in mammalian organs. Nature 478, 343–348 (2011).

    Article  ADS  CAS  Google Scholar 

  21. Kalinka, A. T. et al. Gene expression divergence recapitulates the developmental hourglass model. Nature 468, 811–814 (2010).

    Article  ADS  CAS  Google Scholar 

  22. Khaitovich, P., Enard, W., Lachmann, M. & Pääbo, S. Evolution of primate gene expression. Nat. Rev. Genet. 7, 693–702 (2006).

    Article  CAS  Google Scholar 

  23. Levin, M. et al. The mid-developmental transition and the evolution of animal body plans. Nature 531, 637–641 (2016).

    Article  ADS  CAS  Google Scholar 

  24. Reilly, S. K. & Noonan, J. P. Evolution of gene regulation in humans. Annu. Rev. Genomics Hum. Genet. 17, 45–67 (2016).

    Article  CAS  Google Scholar 

  25. Tirosh, I., Weinberger, A., Carmi, M. & Barkai, N. A genetic signature of interspecies variations in gene expression. Nat. Genet. 38, 830–834 (2006).

    Article  CAS  Google Scholar 

  26. Haberle, V. & Lenhard, B. Promoter architectures and developmental gene regulation. Semin. Cell Dev. Biol. 57, 11–23 (2016).

    Article  CAS  Google Scholar 

  27. Lenhard, B., Sandelin, A. & Carninci, P. Metazoan promoters: emerging characteristics and insights into transcriptional regulation. Nat. Rev. Genet. 13, 233–245 (2012).

    Article  CAS  Google Scholar 

  28. Franz, K. M. & Kagan, J. C. Innate immune receptors as competitive determinants of cell fate. Mol. Cell 66, 750–760 (2017).

    Article  CAS  Google Scholar 

  29. Satija, R. & Shalek, A. K. Heterogeneity in immune responses: from populations to single cells. Trends Immunol. 35, 219–229 (2014).

    Article  CAS  Google Scholar 

  30. Newman, J. R. et al. Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature 441, 840–846 (2006).

    Article  ADS  CAS  Google Scholar 

  31. Faure, A. J., Schmiedel, J. M. & Lehner, B. Systematic analysis of the determinants of gene expression noise in embryonic stem cells. Cell Syst. 5, 471–484.e474 (2017).

    Article  CAS  Google Scholar 

  32. Rand, U. et al. Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Mol. Syst. Biol. 8, 584 (2012).

    Article  Google Scholar 

  33. Fumagalli, M. & Sironi, M. Human genome variability, natural selection and infectious diseases. Curr. Opin. Immunol. 30, 9–16 (2014).

    Article  CAS  Google Scholar 

  34. Johnson, W. E. & Sawyer, S. L. Molecular evolution of the antiretroviral TRIM5 gene. Immunogenetics 61, 163–176 (2009).

    Article  CAS  Google Scholar 

  35. Choo, S. W. et al. Pangolin genomes and the evolution of mammalian scales and immunity. Genome Res. 26, 1312–1322 (2016).

    Article  CAS  Google Scholar 

  36. Braun, B. A., Marcovitz, A., Camp, J. G., Jia, R. & Bejerano, G. Mx1 and Mx2 key antiviral proteins are surprisingly lost in toothed whales. Proc. Natl Acad. Sci. USA 112, 8036–8040 (2015).

    Article  ADS  CAS  Google Scholar 

  37. Xu, L. et al. Loss of RIG-I leads to a functional replacement with MDA5 in the Chinese tree shrew. Proc. Natl Acad. Sci. USA 113, 10950–10955 (2016).

    Article  CAS  Google Scholar 

  38. Sackton, T. B., Lazzaro, B. P. & Clark, A. G. Rapid expansion of immune-related gene families in the house fly, Musca domestica. Mol. Biol. Evol. 34, 857–872 (2017).

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Brunette, R. L. et al. Extensive evolutionary and functional diversity among mammalian AIM2-like receptors. J. Exp. Med. 209, 1969–1983 (2012).

    Article  CAS  Google Scholar 

  40. Malfavon-Borja, R., Wu, L. I., Emerman, M. & Malik, H. S. Birth, decay, and reconstruction of an ancient TRIMCyp gene fusion in primate genomes. Proc. Natl Acad. Sci. USA 110, E583–E592 (2013).

    Article  CAS  Google Scholar 

  41. Barber, M. F., Lee, E. M., Griffin, H. & Elde, N. C. Rapid evolution of primate type 2 immune response factors linked to asthma susceptibility. Genome Biol. Evol. 9, 1757–1765 (2017).

    Article  CAS  Google Scholar 

  42. Saeed, R. & Deane, C. M. Protein–protein interactions, evolutionary rate, abundance and age. BMC Bioinformatics 7, 128 (2006).

    Article  Google Scholar 

  43. Calderone, A., Licata, L. & Cesareni, G. VirusMentha: a new resource for virus-host protein interactions. Nucleic Acids Res. 43, D588–D592 (2015).

    Article  CAS  Google Scholar 

  44. Halehalli, R. R. & Nagarajaram, H. A. Molecular principles of human virus protein-protein interactions. Bioinformatics 31, 1025–1033 (2015).

    Article  CAS  Google Scholar 

  45. Pichlmair, A. et al. Viral immune modulators perturb the human molecular network by common and unique strategies. Nature 487, 486–490 (2012).

    Article  ADS  CAS  Google Scholar 

  46. Dyer, M. D., Murali, T. M. & Sobral, B. W. The landscape of human proteins interacting with viruses and other pathogens. PLoS Pathog. 4, e32 (2008).

    Article  Google Scholar 

  47. Tirosh, I. & Barkai, N. Two strategies for gene regulation by promoter nucleosomes. Genome Res. 18, 1084–1091 (2008).

    Article  CAS  Google Scholar 

  48. Crow, Y. J. & Manel, N. Aicardi-Goutières syndrome and the type I interferonopathies. Nat. Rev. Immunol. 15, 429–440 (2015).

    Article  CAS  Google Scholar 

  49. Hall, J. C. & Rosen, A. Type I interferons: crucial participants in disease amplification in autoimmunity. Nat. Rev. Rheumatol. 6, 40–49 (2010).

    Article  CAS  Google Scholar 

  50. Tisoncik, J. R. et al. Into the eye of the cytokine storm. Microbiol. Mol. Biol. Rev. 76, 16–32 (2012).

    Article  CAS  Google Scholar 

  51. Kilpinen, H. et al. Common genetic variation drives molecular heterogeneity in human iPSCs. Nature 546, 370–375 (2017).

    Article  ADS  CAS  Google Scholar 

  52. Schmidl, C., Rendeiro, A. F., Sheffield, N. C. & Bock, C. ChIPmentation: fast, robust, low-input ChIP-seq for histones and transcription factors. Nat. Methods 12, 963–965 (2015).

    Article  CAS  Google Scholar 

  53. Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protocols 9, 171–181 (2014).

    Article  CAS  Google Scholar 

  54. Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419 (2017).

    Article  CAS  Google Scholar 

  55. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).

    Article  CAS  Google Scholar 

  56. Nourmohammad, A. et al. Adaptive evolution of gene expression in Drosophila. Cell Reports 20, 1385–1395 (2017).

    Article  CAS  Google Scholar 

  57. Zhang, H. M. et al. AnimalTFDB: a comprehensive animal transcription factor database. Nucleic Acids Res. 40, D144–D149 (2012).

    Article  CAS  Google Scholar 

  58. Binns, D. et al. QuickGO: a web-based tool for Gene Ontology searching. Bioinformatics 25, 3045–3046 (2009).

    Article  CAS  Google Scholar 

  59. Kent, W. J. et al. The human genome browser at UCSC. Genome Res. 12, 996–1006 (2002).

    Article  CAS  Google Scholar 

  60. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).

    Article  CAS  Google Scholar 

  61. Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).

    Article  Google Scholar 

  62. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).

    Article  Google Scholar 

  63. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).

    Article  CAS  Google Scholar 

  64. Kuhn, R. M. et al. The UCSC genome browser database: update 2007. Nucleic Acids Res. 35, D668–D673 (2007).

    Article  CAS  Google Scholar 

  65. Mathelier, A. et al. JASPAR 2016: a major expansion and update of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 44, D110–D115 (2016).

    Article  CAS  Google Scholar 

  66. Grant, C. E., Bailey, T. L. & Noble, W. S. FIMO: scanning for occurrences of a given motif. Bioinformatics 27, 1017–1018 (2011).

    Article  CAS  Google Scholar 

  67. Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R. & Siepel, A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 20, 110–121 (2010).

    Article  CAS  Google Scholar 

  68. Kolodziejczyk, A. A. et al. Single cell RNA-sequencing of pluripotent states unlocks modular transcriptional variation. Cell Stem Cell 17, 471–485 (2015).

    Article  CAS  Google Scholar 

  69. Vallejos, C. A., Marioni, J. C. & Richardson, S. BASiCS: Bayesian analysis of single-cell sequencing data. PLOS Comput. Biol. 11, e1004333 (2015).

    Article  ADS  Google Scholar 

  70. Martinez-Jimenez, C. P. et al. Aging increases cell-to-cell transcriptional variability upon immune stimulation. Science 355, 1433–1436 (2017).

    Article  ADS  CAS  Google Scholar 

  71. Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P. L. & Ideker, T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27, 431–432 (2011).

    Article  CAS  Google Scholar 

  72. Lindblad-Toh, K. et al. A high-resolution map of human evolutionary constraint using 29 mammals. Nature 478, 476–482 (2011).

    Article  CAS  Google Scholar 

  73. Herrero, J. et al. Ensembl comparative genomics resources. Database 2016, bav096 (2016).

    Article  Google Scholar 

  74. De Bie, T., Cristianini, N., Demuth, J. P. & Hahn, M. W. CAFE: a computational tool for the study of gene family evolution. Bioinformatics 22, 1269–1271 (2006).

    Article  Google Scholar 

  75. Capra, J. A., Williams, A. G. & Pollard, K. S. ProteinHistorian: tools for the comparative analysis of eukaryote protein origin. PLOS Comput. Biol. 8, e1002567 (2012).

    Article  ADS  CAS  Google Scholar 

  76. Szklarczyk, D. et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452 (2015).

    Article  CAS  Google Scholar 

  77. Zheng, G. X. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).

    Article  ADS  CAS  Google Scholar 

  78. Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33, 495–502 (2015).

    Article  CAS  Google Scholar 

  79. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold-change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

    Article  Google Scholar 

Download references

Acknowledgements

We thank N. Eling, M. Fumagalli, Y. Gilad, O. Laufman, A. Marcovitz, J. Marioni, K. Meyer, M. Muffato, D. Odom, O. Stegle, A. Stern, M. Stubbington, V. Svensson and M. Ward for discussions; G. Emerton, A. Jinat, L. Mamanova, K. Polanski, A. Fullgrabe, N. George, S. Barnett, R. Boyd, S. Patel and C. Gomez for technical assistance; the Hipsci consortium for human fibroblast lines; and members of the Teichmann laboratory for support at various stages. This project was supported by ERC grants (ThDEFINE, ThSWITCH) and an EU FET-OPEN grant (MRG-GRAMMAR No 664918) and Wellcome Sanger core funding (Grant No WT206194). T.H. was supported by an HFSP Long-Term Fellowship and by EMBO Long-Term and Advanced fellowships. V.P. is funded by Fondazione Umberto Veronesi.

Reviewer information

Nature thanks L. Barreiro, I. Yanai and the other anonymous reviewer(s) for their contribution to the peer review of this work.

Author information

Authors and Affiliations

Authors

Contributions

T.H. and S.A.T. designed the project; T.H., X.C., R.J.M., R.R., N.K. and J.-E.P. performed experiments with help from V.P., G.D. and F.A.V.B.; T.H., X.C., R.J.M., R.R., T.G. and J.H. analysed the data with help from G.N., L.B.-C., G.T, A.N. and M.L.; J.F., E.S., P.V., I.K., M.D. and M.H. provided samples; S.A.T. supervised the project; T.H., R.R., N.K. and S.A.T. wrote the manuscript with input from all authors.

Corresponding authors

Correspondence to Tzachi Hagai or Sarah A. Teichmann.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Fibroblast response to dsRNA and IFNB across species.

To study the similarity in response to treatment across species, we plotted the fold-change values of all expressed genes (with one-to-one orthologues) between pairs of species (human–macaque, mouse–rat and human–mouse) in response to dsRNA (poly(I:C)) (ac). As a control, we performed the same procedure with IFNB stimulations (df). Fold-changes were inferred from differential expression analyses, determined by the exact test in the edgeR package6 and based on n = 6, 5, 3 and 3 individuals from human, macaque, rat and mouse, respectively. Spearman correlations between all expressed one-to-one orthologues are shown in grey, Spearman correlations between the subset of differentially expressed genes (FDR-corrected P < 0.01 in at least one species) appear in black. Number of genes shown is n = 11,035, 11,005, 11,137, 10,851, 10,826 and 10,957 in af, respectively. Genes are coloured blue if they were differentially expressed (FDR-corrected P < 0.01) in both species, purple if they were differentially expressed in only one species, or red if they were not differentially expressed. g, h, Density plots of ratio of fold-change in response to dsRNA or to IFNB. g, Comparison between human and macaque orthologues in dsRNA response. h, Comparison between human and mouse orthologues in IFNB response. i, Dendrogram based on the fold-change in response to dsRNA or to IFNB across 9,835 one-to-one orthologues in human, macaque, rat and mouse.

Extended Data Fig. 2 Correspondence of transcriptional divergence and divergence of active promoters and enhancers.

Comparison of divergence in transcriptional response to dsRNA with divergence of active chromatin marks in active promoters (a, profiled using H3K4me3 in proximity to gene’s TSS) and enhancers (b, H3K27ac without overlapping H3K4me3). Chromatin marks were linked to genes on the basis of their proximity to the gene’s TSS. Chromatin marks were obtained from n = 3 individuals in each of the four species, from fibroblasts stimulated with dsRNA or left untreated. The statistics are based on n = 855, 818 and 813 human genes that have a linked H3K4me3 mark with a syntenic region in macaque, rat and mouse, respectively (a); and on n = 326, 241 and 242 human genes that have a linked H3K27ac mark with a syntenic region in macaque, rat and mouse, respectively (b). Each panel shows the fraction of conserved marks between human and macaque, rat or mouse, in genes that have high, medium and low divergence in their transcriptional response. In each column, the histone mark’s signal was compared between human and the syntenic region in one of the three other species. If an active mark was found in the corresponding syntenic region, the linked gene was considered to have a conserved active mark (promoter or enhancer). The fractions of genes with conserved promoters (or enhancers) in each pair of species were compared between high- and low-divergence genes using a one-sided Fisher’s exact test. When comparing active promoter regions of high- versus low-divergence genes, we observe that low-divergence genes have a significantly higher fraction of conserved marks in rodents. This suggests an agreement between divergence at the transcriptional and chromatin levels in active promoter regions. In active enhancer regions, we do not observe these patterns, suggesting that the major contribution to divergence comes from promoters.

Extended Data Fig. 3 Comparison of response divergence of genes containing various promoter elements.

Comparison of response divergence between genes with and without a TATA-box and a CGI. Left, fibroblasts (n = 14, 14, 633 and 294 differentially expressed genes with only TATA-box element, with both CGI and TATA-box elements, with only CGI, and with neither element in their promoters, respectively); right, phagocytes (n = 13, 29, 1,718 and 576 differentially expressed genes with only a TATA-box element, with both CGI and TATA-box elements, with only a CGI, and with neither element in their promoters, respectively). Genes with a TATA-box without a CGI have higher response divergence than genes with both elements. Genes with a CGI but without a TATA-box diverge more slowly than genes with both elements. Genes with both elements do not differ significantly in their divergence from genes lacking both elements (one-sided Mann–Whitney test). Data in boxplots represent the median, first quartile and third quartile with lines extending to the furthest value within 1.5 of the IQR.

Extended Data Fig. 4 Response divergence of molecular processes upregulated in immune response.

Left, distributions of divergence values of n = 955 dsRNA-responsive genes in fibroblasts and subsets of this group belonging to different biological processes. For each functional subset, the distribution of divergence values is compared with the set of 955 dsRNA-responsive genes using a one-sided Mann–Whitney test. FDR-corrected P values are shown above each group and group size is shown inside each box. Right, distributions of divergence values of n = 2,336 LPS-responsive genes in mononuclear phagocytes and subsets of this group belonging to different biological processes. For each functional subset, the distribution of divergence values is compared with the set of 2,336 LPS-responsive genes. FDR-corrected P values (one-sided Mann–Whitney test) are shown above each group and group size is shown inside each box. Data in boxplots represent the median, first quartile and third quartile with lines extending to the furthest value within 1.5 of the IQR.

Extended Data Fig. 5 Cell-to-cell variability versus response divergence across species and conditions in fibroblasts after dsRNA stimulation.

Cell-to-cell variability values, as measured with DM across individual cells, compared with response divergence between species (grouped into low, medium and high divergence). Variability values are based on n = 29, 56, 55, 35 human cells, n = 20, 32, 29, 13 rhesus cells, n = 33, 70, 65, 40 rat cells, and n = 53, 81, 59, 30 mouse cells, stimulated with dsRNA for 0, 2, 4 and 8 h, respectively. Rows represent different dsRNA stimulation time points (0, 2, 4 and 8 h), and columns represent different species as shown. High-divergence genes were compared with low-divergence genes using a one-sided Mann–Whitney test. Data in boxplots represent the median, first quartile and third quartile with lines extending to the furthest value within 1.5 of the IQR.

Extended Data Fig. 6 Cell-to-cell variability versus response divergence across species and conditions in mononuclear phagocytes after LPS stimulation.

Cell-to-cell variability values, as measured with DM across cells, compared with response divergence between species (grouped into low, medium and high divergence). Variability values are based on n = 3,519, 4,321, 3,293, 2,126 mouse cells, n = 2,266, 2,839, 1,963, 1,607 rat cells, n = 3,275, 1,820, 1,522, 1,660 rabbit cells, and n = 1,748, 1,614, 1,899, 1,381 pig cells, stimulated with LPS for 0, 2, 4 and 6 h, respectively. Rows represent different LPS stimulation time points (0, 2, 4 and 6 h), and columns represent different species as shown. High-divergence genes were compared with low-divergence genes using a one-sided Mann–Whitney test. Data in boxplots represent the median, first quartile and third quartile with lines extending to the furthest value within 1.5 of the IQR.

Extended Data Fig. 7 Cell-to-cell variability of cytokine expression in single cell in situ RNA hybridization assay combined with flow cytometry (PrimeFlow).

PrimeFlow measurement of two cytokine genes (IFNB and CXCL10) that show high cell-to-cell variability in scRNA-seq. As controls, two genes matched on expression levels (ATXN2L and ADAM32) but that show low cell-to-cell variability in scRNA-seq data are shown. As the expression of cytokines is at the low end of the distribution, we also chose two genes with middle-range expression values (ADAMTSL3 and BRD2) as additional controls. The experiment was performed in n = 2 independent replicates, originating from the same individual. Both replicates are shown. a, Pseudocolour contour plot for RNA target expression in dsRNA-stimulated human fibroblasts. The x-axis shows area of side scatter (SSC-A) and the y-axis shows fluorescent signal for target RNA probes. RNA targets detected by the same fluorescent channel are displayed together. Top, IFNB and control genes BRD2 and ATXN2L, type 1 probe, Alexa FluorTM 647. Bottom, CXCL10 and control genes ADAMTSL3 and ADAM32, type 10 probe, Alexa FluorTM 568. The cytokine genes display a broader range of fluorescence signal than the controls. b, Histograms comparing fluorescence of cytokine and control pairs (IFNBBRD2 for type 1 probe and CXCL10ADAM32 for type 10 probe). The histograms show a bimodal distribution of expression signal for the two cytokine genes (IFNB and CXCL10, red), but not for controls (blue). This agrees with scRNA-seq data in which CXCL10 and IFNB display high levels of cell-to-cell variability.

Extended Data Fig. 8 Cell-to-cell variability levels and response divergence of cytokines, transcription factors and kinases in response to LPS stimulation of phagocytes.

A scatter plot showing divergence in response to LPS across species and transcriptional cell-to-cell variability in mouse mononuclear phagocytes following 4 h of LPS treatment, in n = 2,262 LPS-responsive genes. Purple, cytokines; green, transcription factors; beige, kinases. The distributions of divergence values and cell-to-cell variability values of each of the three functional groups are shown above and to the right of the scatter plot, respectively.

Extended Data Fig. 9 Cell-to-cell variability levels in cytokines, transcription factors and kinases across species and stimulation time points.

Violin plots showing the distribution of cell-to-cell variability values (DM) of cytokines, transcription factors and kinases during immune stimulation. Left, fibroblast dsRNA stimulation time course. Number of cells used in each species (at 2, 4, 8 h dsRNA, respectively): human, 56, 55, 35; macaque, 32, 29, 13; rat, 70, 65, 40; mouse, 81, 59, 30. Right, phagocyte LPS stimulation time course. Number of cells used in each species (at 2, 4, 6 h LPS, respectively): mouse, 4,321, 3,293, 2,126; rat, 2,839, 1,963, 1,607; rabbit, 1,820, 1,522, 1,660; pig, 1,614, 1,899, 1,381. For both panels, colours as in Fig. 3c. Comparisons between groups of genes were performed using one-sided Mann–Whitney tests. Violin plots show the kernel probability density of the data.

Extended Data Fig. 10 Percentage of cells expressing cytokines, transcription factors and kinases.

Histograms showing the percentage of fibroblasts expressing cytokines (top), transcription factors (middle) and kinases (bottom) following 4 h dsRNA stimulation, in human, macaque, rat and mouse cells (based on n = 55, 29, 65 and 59 cells, respectively). The percentage of expressing cells is divided into 13 bins (x-axis). The y-axis represents the fraction of genes from this gene class (for example, cytokines) that are expressed in each bin (for example, in human, nearly 30% of the cytokine genes (y-axis) are expressed in the first bin, corresponding to expression in fewer than 8% of cells).

Supplementary information

Supplementary Information

This file contains Supplementary Methods, Analyses and Discussion, Supplementary Tables 1-2, 5-7 and Supplementary Figures 1-17

Reporting Summary

Supplementary Table 3

GO-term enrichment analysis of genes in dsRNA treatment – see Supplementary Information for full description

Supplementary Table 4

A list of expressed genes in fibroblasts and phagocytes across species – see Supplementary Information for full description

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hagai, T., Chen, X., Miragaia, R.J. et al. Gene expression variability across cells and species shapes innate immunity. Nature 563, 197–202 (2018). https://doi.org/10.1038/s41586-018-0657-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41586-018-0657-2

Keywords

This article is cited by

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing