An Integrative Framework Reveals Signaling-to-Transcription Events in Toll-like Receptor Signaling

Abstract

Building an integrated view of cellular responses to environmental cues remains a fundamental challenge due to the complexity of intracellular networks in mammalian cells. Here, we introduce an integrative biochemical and genetic framework to dissect signal transduction events using multiple data types and, in particular, to unify signaling and transcriptional networks. Using the Toll-like receptor (TLR) system as a model cellular response, we generate multifaceted datasets on physical, enzymatic, and functional interactions and integrate these data to reveal biochemical paths that connect TLR4 signaling to transcription. We define the roles of proximal TLR4 kinases, identify and functionally test two dozen candidate regulators, and demonstrate a role for Ap1ar (encoding the Gadkin protein) and its binding partner, Picalm, potentially linking vesicle transport with pro-inflammatory responses. Our study thus demonstrates how deciphering dynamic cellular responses by integrating datasets on various regulatory layers defines key components and higher-order logic underlying signaling-to-transcription pathways.

Keywords: TLRs; Toll-like receptors; large-scale in vitro kinase assay; pathogen-sensing pathways; phosphoproteomics; protein-protein interactions; signaling; transcriptional network analysis.

Figure 1. TLR4 stimulation with LPS leads to global and dynamic changes in the phosphoproteome of dendritic cells (DCs)

(A–B) Diagram highlighting general principles of cellular signaling-to-transcription events (A) and their transposition to the TLR4 pathway (B). (C) Temporal changes in the phosphoproteome of LPS-stimulated DCs. Shown are the distributions of log2 fold changes of phosphosites (X axis) between LPS-treated and untreated cells at indicated times after LPS stimulation, as density (top of each panel) and dot plots (bottom of each panel, with MS2 spectra count in Y axis and showing phosphosites measured in all 8 time points). (D) Comparison between the phosphoproteome and total proteome of LPS-stimulated DCs. Shown are distributions of log2 fold changes of phosphosites (X axis) and proteins (Y axis) between LPS-treated and untreated cells at 120 (top) and 360 (bottom) min post-stimulation. See also Figure S1 and Table S1.

Figure 2. Temporal analysis of the LPS-induced phosphoproteome reveals known and candidate regulators of TLR4 signaling

(A) Temporal phosphorylation profiles during LPS stimulation in DCs. Log2 fold changes between LPS-treated and untreated cells for 3557 phosphosites (rows) detected in at least 6 out of 8 time points (columns). Phosphosites are partitioned into 10 clusters using k-means (color bars, right). White indicates missing values. (B) Median log2 fold changes between LPS-treated and untreated cells (y axis) and median absolute deviation (MAD, colored error bar) at each time point (x axis) for phosphosites in all 10 k-means clusters from A. Known TLR pathway proteins detected in each cluster are indicated on the right. Parentheses indicate the number of phosphosites per proteins (when > 1). See also Figure S2 and Table S2.

Figure 3. Genetic perturbations of phosphorylated proteins identify putative regulators of TLR4 signaling

(A) Overview of phosphoprotein candidate selection for functional analysis. (B) Perturbation profiles of the 27 phosphoproteins that significantly impacted TLR4 outputs. Shown are the perturbed candidates and control phosphoproteins (columns) and the log2 fold changes for each target gene (rows) between gene-specific and control shRNAs. The right-most column categorizes target genes into antiviral (light green) and inflammatory (light orange) programs. See also Figure S3 and Table S3.

Figure 4. Identification of candidate regulators in the MYD88-dependent inflammatory pathway

(A) Perturbation profiles of genes affecting the MYD88 pathway. Shown are 4 perturbed candidate genes and MYD88 (columns) and the log2 fold changes between gene-specific and control shRNAs (rows) of 10 target genes. The right-most column categorizes target genes into antiviral (light green) and inflammatory (light orange) programs. (B) Expression levels (qPCR) relative to control shRNAs (left bars, dark grey) for two antiviral cytokines (Ifit1 and Cxcl10) and for three inflammatory cytokines (Il6, Cxcl1, and Tnf), following LPS stimulation in DCs using two independent shRNAs. Bottom tick marks separate shRNAs controls and each gene (‘Average’ indicate the mean value for all 8 control shRNAs). Two to three replicates for each experiment; error bars are the standard deviations. (C) Inhibition of transcription of inflammation cytokines in Ap1ar^−/− DCs. mRNA levels (qPCR; relative to Gapdh) for indicated inflammatory (light orange) and antiviral (light green) cytokines in three replicates per time point. Error bars represent the standard deviation. (D) Interaction proteomics identified putative binders for AP1AR in DCs. Log2 fold change (X axis) of proteins identified between DCs expressing V5-tagged-AP1AR and -GFP (control bait) plotted against the number of peptides identified per protein (Y axis). (E) Perturbation profiles of indicated genes (columns) and the log2 fold changes between gene-specific and control shRNAs (rows) of 150 target genes. The right-most column categorizes target genes into antiviral (light green) and inflammatory (light orange) programs. (F) Impact of SAMHD1 mutations on human fibroblast cell response to LPS. Human fibroblasts from healthy (H) or mutant-carrying patients (M; with homozygous c.445C>T p.Gln149* for M1 and c.1609-1G>C for M2) were stimulated with LPS or left untreated as control, and indicated inflammatory (light orange) and antiviral (light green) cytokine levels were measured by qPCR (relative to GAPDH). Error bars represent the standard deviation. See also Figure S4 and Table S4.

Figure 5. Similarities in perturbation profiles and overlap with transcription factor target genes suggest three functional modules for the 29 candidate phosphoproteins

(A) Functional classification based on similarity of perturbation profiles. Shown is a correlation matrix (Pearson correlation coefficient) of the perturbation profiles from Figures 3B and 4E combined. (B) Intersection between genes affected by a phosphoprotein perturbation and genes whose promoters are bound by transcription factors (TFs). Shown are the overlaps between genes affected by 29 candidate signaling regulators knockdowns (columns, including positive control genes) and genes whose promoters are bound by 20 TFs (rows). P values, hypergeometric test (purple: significant correlation; white: no correlation). See also Figure S5 and Table S5.

Figure 6. Physical and functional proteomics assays pinpoint binding and phosphorylation events downstream of the Myd88 adaptor and associated kinases

(A–B) Affinity purification followed by proteomics. Shown are dot plots of SILAC ratios for proteins identified in DCs overexpressing V5-tagged MYD88 (A) or IRAK2 (B). Cells were stimulated with LPS for 30 min, and protein complexes purified using anti-V5 antibodies coupled to magnetic beads. Each axis represents an independent experiment. (C) Diagram depicting our experimental approach for measuring the impact of gene knockout on the TLR4-regulated phosphoproteome of mouse BMDCs. (D) Phosphoproteomics in KO cells. Left, shown is a heatmap for SILAC ratios of phosphosites (rows) in 4 KO models (columns) at 30 min after LPS stimulation compared to control wild-type cells, as indicated (grey, missing values). Middle, shown in light brown are phosphosites with significant up- or down-regulation in KO vs WT. Right, shown in black are the phosphosites belonging to known TLR proteins. (E) Diagram depicting our experimental approach for large-scale in vitro kinase assays using native protein lysates from BMDCs and phosphoproteomics. (F–G) In vitro kinase (IVK) assay followed by phosphoproteomics. Shown are scatter plots of SILAC ratios of phosphosites identified using purified kinases: IRAK4 (D) and TBK1 (E). Light grey, all data points; dark grey, phosphosites with FDR < 0.1 in IVK; red, phosphosites with FDR < 0.1 in both IVK and in cells stimulated with LPS, which highlights the overlap between IVK and phosphoproteome measurements on stimulated cells (denoted as IVK + cells). Gene names at the bottom right of each plot indicate known TLR components with the number of phosphosites in parenthesis. See also Figure S6 and Table S6.

Figure 7. An integrative analysis reveals known and candidate signaling-to-transcription paths and helps parse the effects of Myd88 and associated kinases in the TLR4 system

(A) A computational framework for integrative analysis of the functional and physical proteomics datasets collected in this study (from left to right): a background interaction network is assembled using database and local data, nodes and edges are scored based on experimental evidence from this work, and statistically significant relationships determined by bootstrap analysis. (B) Cumulative number of significant relationships (boostrap p-value < 0.0005, FDR < 0.05) identified between ‘seed nodes’ (29) and any of the transcriptional regulators detected in BMDCs (782 possible ‘target nodes’ in total) using ‘background network’ (dark grey) and ‘weighted network’ (light grey) methods. (C) Total number of relationships linking seeds (29) and known TLR transcription regulators (14) for ‘background network’ (dark grey) and ‘weighted network’ (light grey) methods. (D) Significant relationships (420 pairs) found between 29 seeds (columns) and 95 transcriptional regulators (rows). Modules from Figures 3B and 5A are shown (columns) in light green (I), purple (II), and orange (III). Transcriptional regulators with phosphosites with significant up- or down-regulation in Myd88^−/−/Ticam1^−/− vs WT and in time series are indicated on the right (light brown). P-values, bootstrap (purple). (E) An interaction network connects 27 seeds (blue) to 95 transcriptional regulators (red) through the top 60 intermediate (yellow) nodes that were ranked based on centrality measure (see Experimental Procedures). (F) Centrality score of the top 60 intermediate nodes across the three modules from D. (G) t-distributed stochastic neighbor embedding (t-SNE) analysis of the effects of gene knockout (data from Figure 6D) on the phosphorylation levels of nodes present in the paths mediating the seed-transcriptional regulator relationships identified in panel D. Shown are all of the 391 out of 420 relationships affected by Myd88^−/−/Ticam1^−/− (grey dots). The effects of Irak4, Myd88 and Irak2 on these paths are overlaid in orange, blue and yellow, respectively. See also Figure S7 and Table S7.