Mice

Female C57BL/6J mice (wild-type, stock 000664), B6.129S7-Ifngtm1Ts/J (Ifng KO, stock 002287), B6.129P2-Il18tm1Aki/J (Il18 KO, stock 004130), C57BL/6J-Il1bem2Lutzy/Mmjax (Il1b KO, stock 068082-JAX) and B6.129S-Tnftm1Gkl/J (Tnf KO, stock 005540) were obtained from the Jackson Laboratories. Animals were housed in specific pathogen-free and BSL2 conditions at The University of Chicago, and all experiments were performed in accordance with the US National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by The University of Chicago Institutional Animal Care and Use Committee.

Mouse models of endotoxemia and sepsis

For LPS endotoxemia, mice were injected intraperitoneally with either lethal (10–15 mg per kg body weight) or sublethal (3–5 mg per kg body weight) doses of LPS derived from Escherichia coli O55:B5 (Sigma-Aldrich) diluted in PBS. Dosing was established for each lot of LPS by in vivo titration. CLP was performed as described by others59,60. Briefly, mice were anesthetized with isoflurane. A 1- to 2-cm midline laparotomy was performed and the cecum was exposed. The cecum was ligated with 6-0 silk sutures (Ethicon) and perforated as follows to vary disease severity: (1) mild sepsis: ligate at distal 33% position and perforate once with a 21-gauge needle; (2) moderate sepsis: ligate at distal 40% position and perforate twice with a 19-gauge needle; (3) severe sepsis: ligate immediately below the ileocecal valve and perforate twice with a 19-gauge needle. The cecum was tucked back into the peritoneum and gently squeezed to extrude a small amount of fecal content. The peritoneal wall was closed using absorbable suture. The skin was closed with surgical staples. To resuscitate animals, 1 ml of saline was injected subcutaneously. Mice were temporarily placed on a heating pad for recovery. Sham-operated mice underwent the same procedure except that the cecum was neither tied nor perforated.

Recombinant cytokine injections

C57BL/6J mice were injected intravenously with 2.5 µg of recombinant TNF, IL-1β, IL-6, IL-10, IL-18 or IFN-γ used alone (6 singles) or in pairwise combinations (15 pairs).

Neutralizing antibody and drug treatments

For neutralizing antibodies, C57BL/6J or indicated knockout mice were injected intraperitoneally with 50 µg of TNF (clone BE0058, BioXCell), IL-18 (clone BE0237, BioXCell), IFN-γ (clone BE0055, BioXCell) or IL-1β (clone BE0246, BioXCell) neutralizing antibodies in 100 µl of PBS 1 h before LPS injection.

Blood analysis

Mouse whole blood was harvested by cardiac puncture and plasma and serum were isolated using lithium heparin-coated Microtainer blood collection tubes (BD, 365965) and Microtainer blood collection tubes (BD, 365978), respectively. For flow cytometric, bead-based immunoassays, plasma was diluted and processed using the LEGENDplex Mouse Inflammation Panel (BioLegend, 740446) and Mouse Macrophage/Microglia Panel (BioLegend, 740846) kits. Data were acquired on the NovoCyte flow cytometer (Acea Biosciences/Agilent) and analyzed using the LEGNEDplex software v8 (BioLegend). To measure tissue injury marker levels in sera, samples were processed with the following kits for BUN (BioAssay Systems DIUR-100), ALT (Cayman Chemical, 700260) and troponin-I (Life Diagnostics, CTNI-1-HS) levels according to the manufacturer’s instructions.

Tissue harvest

Tissues were harvested, frozen and stored as previously described16,17. Mice were anesthetized with 2,2,2-tribromoethanol (250–500 mg per kg body weight) and perfused transcardially with PBS containing 10 mM EDTA (to avoid signal contamination from blood in tissues). Before perfusion, blood was collected by cardiac puncture and stored on ice and, immediately after perfusion, tissues were placed in RNA-preserving solution (5.3 M ammonium sulfate, 25 mM sodium citrate, 20 mM EDTA) and kept at 4 °C overnight before transfer at −80 °C for storage. For each mouse, we harvested up to 13 tissues in total: iLNs, flank skin, thymus, heart, lung, spleen, kidney, small intestine, colon, liver, brain, bone marrow and PBMCs. Small intestine and colon were cut longitudinally and washed extensively in PBS to completely remove feces contamination. Bone marrow cells were collected from femora and tibiae, stored overnight in RNA-preserving solution at 4 °C, centrifuged at 5,000g for 5 min at 4 °C and cell pellets were stored at −80 °C.

Whole-tissue RNA extraction

Whole-tissue RNA extraction was performed as described previously17. Briefly, tissues stored in RNA-preserving solution were thawed and transferred to 2-ml tubes containing 700–1,500 µl (depending on tissue) of PureZOL (Bio-Rad, 7326890) or homemade TRIzol-like solution (38% phenol, 0.8 M guanidine thiocyanate, 0.4 M ammonium thiocyanate, 0.1 M sodium acetate, 5% glycerol). Tissues were lysed by adding 2.8-mm ceramic beads (OMNI International, 19–646) and running 1–3 cycles of 5–45 s at 3,500 r.p.m. on the PowerLyzer 24 (QIAGEN). For liver, brain and small intestine samples, tissues were lysed with 3–5 ml using M tubes (Miltenyi Biotec, 130-096-335) and running 1–4 cycles of the RNA_02.01 program on the gentleMACS Octo Dissociator (Miltenyi Biotec). Next, lysates were processed in deep 96-well plates (USA Scientific, 1896–2000) by adding chloroform for phase separation by centrifugation, followed by precipitation of total RNA in the aqueous phase using magnetic beads coated with silane (Dynabeads MyOne Silane; Thermo Fisher Scientific, 37002D), buffer RLT (QIAGEN, 79216) and ethanol. Genomic DNA contamination was removed by on-bead DNase I (Thermo Fisher Scientific, AM2239) treatment at 37 °C for 20 min. After washing steps with 80% ethanol, RNA was eluted from beads. This RNA extraction protocol was performed on the Bravo Automated Liquid Handling Platform (Agilent)17. Sample concentrations were measured using a Nanodrop One (Thermo Scientific). RNA quality was confirmed using a Tapestation 4200 (Agilent Technologies). The samples with low RNA quality were excluded from the subsequent experiments.

RNA sequencing

For each tissue sample, full-length cDNA was synthesized in 20 µl final reaction volume containing the following: (1) 10 µl of 10 ng µl–1 RNA; (2) 1 µl containing 2 pmol of a custom RT primer biotinylated in 5′ and containing sequences from 5′ to 3′ for the Illumina read 1 primer, a 6-bp sample barcode (up to 384), a 10-bp unique molecular identifier (UMI) and an anchored oligo(dT)30 for priming61; and (3) 9 µl of RT mix containing 4 µl of 5× RT buffer, 1 µl of 10 mM dNTPs, 2 pmol of template switching oligo and 0.25 µl of Maxima H Minus Reverse Transcriptase (Thermo Scientific, EP0753). First, barcoded RT primers were added to RNA, which were then denatured at 72 °C for 1 min followed by snap cooling on ice. Second, the RT mix was added and plates were incubated at 42 °C for 120 min. For each library, double-stranded cDNA from up to 384 samples were pooled using DNA Clean & Concentrator-5 columns (Zymo Research, D4013) and residual RT primers were removed using exonuclease I (New England Biolabs, M0293). Full-length cDNAs were amplified with 5 to 8 cycles of single-primer PCR using the Advantage 2 PCR Kit (Clontech, 639206) and cleaned up using SPRIselect magnetic beads (Beckman Coulter, B23318). cDNA was quantified with a Qubit dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific, 32851) and 50 ng of cDNA per pool of samples was tagmented using the Tagment DNA Enzyme I (Illumina, 20034197) and amplified using the NEBNext Ultra II Q5 Master Mix (New England BioLabs, M0544L). Libraries were gel purified using 2% E-Gel EX Agarose Gels (Thermo Fisher Scientific, G402002), quantified with a Qubit dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific, Q32851) and a Tapestation 4200 (Agilent Technologies) and sequenced on the NextSeq 550 platform (Illumina).

Custom, whole-mouse spatial transcriptomics using Array-seq

Mice injected with LPS or left untreated as control were euthanized with CO2, frozen in a dry ice–hexane bath after removing all body hair and teeth and stored at −80 °C until use. Frozen section preparation and section transfer were carried out by modifying Kawamoto’s film method62,63,64. Frozen mice were embedded in a cryo-embedding medium and sectioned (10-µm thickness) using a Leica CM3600-XP cryomacrotome. Resulting whole-mouse sections were transferred onto custom, large-format microarrays (30-µm spot diameter with 36.65 µm center-to-center distance between spots), which were repurposed for spatial transcriptomics measurements using the Array-seq method. After transfer, sections were fixed in methanol, stained with H&E and imaged on an Olympus VS2000 slide scanner (×20 magnification). Sections were permeabilized (1% pepsin), incubated for in-tissue reverse transcription and treated with proteinase K for tissue removal. Resulting full-length, single-stranded cDNAs were denatured and retrieved from the array using potassium hydroxide and purified by column clean up (Zymo Research). cDNA was processed for single-primer PCR amplication followed by sequencing library construction using tagmentation (Nextera DNA Library Prep Kit) and final PCR amplification. Resulting libraries were sequenced on the NovaSeq 6000 (Illumina) and sequencing data was preprocessed using STAR/STARsolo (version 2.7.10a)79 (https://github.com/alexdobin/STAR/blob/master/docs/STARsolo.md/) for read alignment using the GRCm39 mouse reference genome, spatial barcode demultiplexing and UMI counting. Resulting spatial transcriptomics data was normalized, processed for differential expression analysis, and visualized using custom Python 3.8.5 (https://www.python.org/) scripts and existing packages, including Scanpy (version 1.9.1)65, scikit-image (version 1.1.3)66 and Seaborn (version 0.11.2)67,68 and scikit-learn (version 0.24.2). Cell-type deconvolution for each spatial transcriptomics spot was done using the CARD package (version 1.0.0)69.

Commercial, kidney spatial transcriptomics using Visium

Mouse kidneys were dissected from LPS-injected or control mice without transcardial perfusion and frozen in optimal cutting temperature media. In total, 10-µm frozen tissue sections were cut with a CM1850 Cryostat (Leica) and mounted onto a Visium Spatial Gene Expression library preparation slide (10x Genomics). Samples were processed to generate spatial transcriptomics sequencing libraries according to the manufacturer’s instructions. In brief, sections were fixed in 100% methanol and stained with H&E reagents. H&E-stained sections were imaged using a CRi Panoramic MIDI Whole Slide Scanner with ×20 magnification. Sections were then permeabilized with 0.1% pepsin in 0.01 M HCl for 14 min at 37 °C and processed for in-tissue reverse transcription followed by on-slide second-strand synthesis. Resulting cDNA was used to construct sequencing libraries that were sequenced on the NextSeq 550 platform (Illumina), with 28 bases for read 1 and 56 for read 2 and at a depth of 78–114 million total reads per sample. The output data of each sequencing run (Illumina BCL files) were converted into FASTQ files using Bcl2Fastq v.2.19.1. The Space Ranger software (v.1.2.0, 10x Genomics) was used to process, align and summarize the FASTQ files against a GRCm39 mouse reference genome. Raw UMI count spot matrices, spot coordinates and images were imported into Python using Scanpy (v.1.9.1)65. Raw UMI counts were log10 normalized and clustered using a Louvain algorithm (resolution of 0.35). Differential expression between control and LPS-treated samples was performed using Scanpy’s rank_genes_groups function using a Wilcoxon rank-sum test. Spatially resolved counts of differentially expressed genes were overlaid with corresponding grayscale H&E images and visualized using Seaborn v.0.11.2 (https://github.com/mwaskom/seaborn/).

Histology

Tissue processing, embedding, sectioning, immunohistochemistry using purified mouse Ly6G (clone 1A8, BioLegend) and F4/80 (clone BM8, BioLegend) antibodies, or TUNEL (terminal deoxynucleotidyl transferase dUTP nick end labeling) staining was performed by the Human Tissue Resource Center at the University of Chicago. Section images were obtained using the Slideview VS200 Research Slide Scanner (Olympus). Image analysis and quantification (Ly6G+, TUNEL+ and F4/80+ areas) were performed using ImageJ (https://imagej.nih.gov/ij/).

Flow cytometry

To analyze splenic B cells, total splenocytes were obtained by mashing spleens on 70-µm filters followed by red blood cell lysis (Lonza). To analysis red blood cell content in the bone marrow, total bone marrow cells were flushed out of femora and tibiae using PBS. Single-cell suspensions were stained in the presence of Fc receptor-blocking antibodies (mouse CD16/32, clone 93) using the following antibodies (BioLegend): CD19-FITC (clone 1D3/CD19, 152403), B220-PerCP (clone RA3-6B2, 103233), CD93-PE (clone AA4.1, 136503), CD23-APC (clone B3B4, 101619), CD21-Pacific Blue (clone 7E9, 123413), Ter119-FITC (clone TER-119, 116205) and CD45-APC-Cy7 (clone 30-F11, 103115). Cell viability was measured using Zombie Yellow Fixable Viability kit (423103) or DAPI. Flow cytometry data were acquired on the NovoCyte flow cytometer (Acea Biosciences/Agilent Technologies) using NovoExpress (version 1.3.0) and analyzed using FlowJo (BD).

RNA-seq data analysis

Sequencing read files were processed to generate UMI70 count matrices using the Python toolkit from the bcbio-nextgen project version 1.1.5 (https://bcbio-nextgen.readthedocs.io/en/latest/). In brief, reads were aligned to the mouse mm10 transcriptome with RapMap71. Quality-control metrics were compiled with a combination of FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), Qualimap and MultiQC (https://github.com/ewels/MultiQC/)72,73. Samples were demultiplexed using barcodes stored in read 1 (first 6 bases), and raw UMI count matrices were computed using UMIs stored in read 1 (bases 7 to 16; https://github.com/vals/umis/).

Differential expression analysis was done using custom scripts in R version 4.2.0 (https://www.r-project.org/). Raw count matrices were filtered to keep genes with at least 20 counts per million or five UMIs in two samples and normalized across samples using the calcNormFactor function in edgeR74. We identified genes with at least a twofold expression difference and indicated Benjamini and Hochberg FDR-adjusted P value and fold expression difference by comparing treated tissues and matching control tissues using limma. Data analysis was also performed with existing packages, including tidyverse (version 2.0.0), data.table (version 1.14.8), cmapR (version 1.8.0), RColorBrewer (version 1.1–3), enrichR (version 3.1), ggrepel (version 0.9.3), patchwork (version 1.1.2), cowplot (version 1.1.0), glue (version 1.4.2), fs (version 1.3.2) and Matrix (version 1.2–18).

To assess the expression profiles of known sepsis biomarkers, we used a set of 258 genes reported as sepsis biomarkers by others24. We defined the absolute average log2 fold change of these 258 genes within each RNA-seq profile as the sepsis biomarker score.

Heat maps for RNA-seq data display the indicated numbers of transcripts, and color intensities are determined by log2 fold-change value for each heat map. The rows of each heat map were ordered by k-means clustering of log2 fold-change values in R or Morpheus (https://software.broadinstitute.org/morpheus/). All heat maps were generated using ComplexHeatmap (version 2.12.1; https://github.com/jokergoo/ComplexHeatmap/) and circlize (version 0.4.15; https://github.com/jokergoo/circlize/) packages in R75,76.

Statistical modeling of cytokine pairwise effects on tissue gene expression

To assess the extent to which pairwise administration of cytokines (that is, TNF plus IL-18, IFN-γ or IL-1β) resulted in nonadditive changes (that is, synergistic or antagonistic interactions) in gene expression levels across tissues in mice, we developed an interaction scoring method based on a linear modeling method adapted from previous work77, using custom scripts in R (https://www.r-project.org/).

First, raw, tissue RNA-seq count matrices were normalized across samples using the calcNormFactor function in edgeR74 and subsequently filtered to keep genes with at least 15 counts per million in two samples. Data were log-transformed and a linear model was fit using the limma package78. We then computed the following contrasts for each pair (AB) of interest and its component singles (A, B) and unstimulated control mice:

$$\begin}\,1\,}=}-},\\ }\,2\,}=}-},\\ }\,}=}-},\\ }\,}=[}-}]+[}-}],}\\ }\,}=(}-})-[(}-})+(}-})].\end$$

Where, ‘pair effect’ is equivalent to the observed gene expression value for a given pair, while ‘additive effect’ is equivalent to the predicted gene expression value for that pair if it is assumed to be equal to the sum of the component singles. The ‘interaction effect’ is equal to the difference between these two values and is used as the score for assessing nonadditive interactions.

We identified genes with significantly different expression within each contrast and across all contrasts using a Benjamini and Hochberg correction for multiple-hypothesis testing and an FDR of 0.1. We next classified each gene, for each organ and pair treatment, as ‘synergistic’, ‘antagonistic’ or ‘additive’, depending on its score and the gene expression values of the pair and its component singles. Genes without >0.5 absolute difference in log2 fold change in at least two of the three experimentally measured treatment conditions for a given pair and organ (single 1 effect, single 2 effect, pair effect) were considered to have roughly the same expression across all samples and were excluded from further classification to avoid classifying genes with very high or low baseline expression in the singles (and, therefore, very high or low predicted additive effects but no additional increase or decrease in gene expression at the pair level) as synergistic or antagonistic.

Using the standard deviation for each contrast determined via limma, we calculated an error value, E, for each gene, as the average of the standard deviations for all experimentally measured samples for that gene. Where the score was >2 × E, and the score was significant (FDR < 0.1) OR score > 1 (>2-fold difference between predicted and observed gene expression values), the gene was classified as synergistic. The gene was only classified as significantly synergistic if the score was determined to be significant at the chosen FDR (0.1). Following the same logic, if the score < −2 × E and score significant (FDR < 0.1) OR score < −1, the gene was classified as antagonistic in a particular pair and organ. Again, only if the score was determined to be significant at the chosen FDR (0.1), was the gene classified as significantly antagonistic.

The total number of DEGs was calculated by totaling any gene that showed significant differential expression (FDR < 0.1) in single 1, single 2 or pair, compared to control. The percentage of all DEGs for a given pair and organ that were classified as synergistic, additive or antagonistic was then calculated.

Public RNA-seq data

To compare the expression profile of bacterial sepsis (this study) with that of viral sepsis induced using tissues from mice infected with a lethal dose of vaccinia virus strain Western Reserve16, we used our previously published bulk RNA-seq data (GSE87633).

Statistics

Statistical analyses were performed by R, using limma, or one-way ANOVA with Tukey–Kramer test. Data collection and analysis were not performed blind to the conditions of the experiments. No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those reported in previous publications16. Data distribution was assumed to be normal, but this was not formally tested. For experiments that require treatments, age-matched and sex-matched animals were randomly assigned into each group.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.