Quantification of dynamic protein interactions and phosphorylation in LPS signaling pathway by SWATH-MS

of the LPS signaling pathway. from macrophage cell line RAW 264.7 upon a serial of LPS treatments, and quantified proteins and phosphorylation events in these complexes using SWATH-MS. Almost all known players involved in LPS pathway have been identified and quantified in one experiment. Dynamic key protein profiles and their phosphorylation events in MyD88, TRAF6 and NEMO complexes under LPS treatment are shown, revealing the complex regulation network in LPS signaling pathway. In addition, we showed the stoichiometry of Myddosome, which can offer an insight into how TIRAP and MyD88 cooperate to ensure signal transduction and how IRAK complex are assembled.

Downstream of IRAK1 and IRAK4, TRAF6 is critical to the MyD88-dependent pathway(6). TRAF6 forms a complex with UBC13 and UEV1A, which activates TAK1(7). Subsequently TAK1 activates downstream NF-B and MAPK pathways. The MyD88-indepent signaling pathway is mediated by TRIF, an important TIRcontaining adaptor protein. TRIF recruits TRAF3 which then associates with TANK (TRAF family member-associated NF-B activator), TBK1 and IKKi to mediate IRF3 dimerization(8-10). IRF3, cooperating with NF-B, activates the transcription of target genes, such as Type I interferons and interferon-inducible genes(11). Although TLR4 signaling pathway has been intensively studied in the past few years, several important questions still remain unsolved. While LPS-induced NF-B activation is abolished in TAK1-deficient MEF cells, it is not affected in TAK1deficient macrophages (12,13). These reports indicate that other unknown players downstream of TRAF6, but not TAK1, are essential for LPS-induced NF-B activation in macrophages. Another intriguing question is how TIRAP-MyD88-IRAK complex (named Myddosome) is assembled and what the stoichiometry of this complex is. SWATH-MS (sequential window acquisition of all theoretical fragment ion spectra mass spectrometry) is a data-independent acquisition MS technique, which combines the advantages of shotgun MS and SRM (Selected Reaction Monitoring)(14, 15). SWATH-MS enables consistent and accurate protein quantification across multiple samples. Besides, SWATH-MS can offer relative protein quantification within one sample, allowing for calculation of protein complex stoichiometry (16). In this study, we used immunoprecipitation to purify time-resolved protein complexes from macrophage cell line RAW 264.7 upon a serial of LPS treatments, and quantified proteins and phosphorylation events in these complexes using SWATH-MS. Almost all known players involved in LPS pathway have been identified and quantified in one experiment. Dynamic key protein profiles and their phosphorylation events in MyD88, TRAF6 and NEMO complexes under LPS treatment are shown, revealing the complex regulation network in LPS signaling pathway. In addition, we showed the stoichiometry of Myddosome, which can offer an insight into how TIRAP and MyD88 cooperate to ensure signal transduction and how IRAK complex are assembled.

Methods Experimental Design and Statistical Rationale
3xFlag-KI TRAF6 RAW 264.7 cell line and 3xFlag-KI MyD88 RAW 264.7 cell line were treated with LPS for ten time points. 3xFlag-NEMO-reconstituted RAW 264.7 cell line was treated with LPS for seven time points. RAW 264.7 cell line expressing 3xFlag-vector was treated with LPS for ten time points. Protein complexes were purified with M2 beads. Four biological replicates were performed for 3xFlag-KI TRAF6 cells, and three biological replicates were performed for 3xFlag-KI MyD88 cells and 3xFlag-NEMO cells. Two biological replicates were performed for 3xFlagvector cells. One biological replicate is the repeat of the whole experiment from cell culture to MS analysis under the same conditions. No technical replicates (multiple injections for one peptide sample) was performed in the study. In the differential expression analysis, proteins with log2(fold change)>1 and -log10(P-value) >1.5 were supernatants were collected for immunoprecipitated overnight with anti-Flag M2 antibody-conjugated agarose at 4 °C. Beads containing protein complexes were washed four times with HBS lysis buffer. Bound Flag-immune complexes were eluted twice with 0.15 mg/mL of 3xFlag peptide with N-terminal biotin tag in HBS lysis buffer and then precipitated with 20% trichloroacetic acid (TCA). The protein pellets were washed three times with 1-mL cold acetone and dried in speedvac. TCA-precipitated proteins were re-suspended in 50 L8 M urea in 50 mM NH4HCO3, and 10 mM Tris(2-carboxyethyl) phosphine hydrochloride (TCEP) and 40 mM chloroacetamide (CAA) were added into reactions for 30 min at 37 °C for cysteine reduction and alkylation. Next, 8 M urea were diluted to 1.6 M urea with 50 mM NH4HCO3 and trypsin was added at the protein:trypin ratio of 50:1. Digestion was performed overnight at 37 °C. The biotin-3xFlag peptide was removed by the avidin beads. Peptides were acidified to a final concentration of 1% formic acid (FA), followed by desaltion using C18 StageTips. After desaltion, peptides were eluted with 70% acetonitrile/1% formic acid and dried. For phosphopeptide enrichment with IMAC, peptides were dissolved in 50 L 60%ACN/1%AA and incubated with 5 uL bead volume of IMAC beads. The peptides with IMAC beads was shaken for 30 min at room temperature. Non-phopshopeptides were washed with 25% ACN/0.1M NaCl/0.1%AA for three times followed by one-time water wash. Phosphopeptides were then eluted with 6% NH3 H2O and dried in speedvac.
Analysis of phosphoproteomics data SWATH-MS wiff files were first analyzed by Group-DIA for generation of mgf files as described above. These mgf files were converted to mzXML files, and DDA wiff files were converted to centroid mzXML files using qtofpeakpicker tool in TPP. These mzXML files were searched with Comet and X!Tandem (native and k-score plugin) altogether. The search parameters were set as followed: parent monoisotopic tolerance 50 ppm, modification 57.021464@C, potential modification mass 15.994915@M 79.966331@STY and maximum missed cleavage sites 2. The pep.xml files were scored using PeptideProphet with parameters -p0.05 -l7 -PPM -OAdPE -dDECOY and combined by iProphet with parameters DECOY=DECOY. The ipro.pep.xml files were subsequently analyzed by PTMProphet for phosphosites localization scoring. The phosphopeptides were first filtered at the iProphet values corresponding to 1% peptide FDR which was determined by Mayu (Version 1.07), and then filtered at localization score 0.7. After these filters, the phosphopeptide ions were for library building by SpectraST with CID-QTOF setting. The retention time of peptides in sptxt file were replaced with iRT time using spectrast2spectrast_irt.py script (downloaded from www.openswath.org), and iRT peptides used for retention time normalization were endogenous peptides. The sptxt file was made consensus nonabundant spectral library with the iRT retention time using spectraST.
Reprocessing of murine cell line L929 total spectral library As described previously, we used Mascot and X!tandem to process the L929 cell line DDA data. However, we found Comet performed better than Mascot. Therefore, we reprocessed 248 DDA runs with Comet and X!tandem. Database searches and spectral library building were performed using the same method as described above. Briefly, the DDA wiff files were first converted to centroid mzXML files using qtofpeakpicker tool in TPP followed by database searches with Comet and X!Tandem (native and k-score plugin). The pepXML files were validated with PeptideProphet followed by iProphet analysis. An iProphet probability cutoff (p=0.956653) was determined by MAYU (version 1.07), resulting in 110,211 peptides (FDR=0.22%) which correspond to 8,599 proteins (FDR=1.08%). The pepXML file was filtered to 1% FDR at protein level and converted to sptxt file using SpectraST with CID-QTOF setting. The retention time of peptides in sptxt file were replaced with iRT time using spectrast2spectrast_irt.py script, where ciRT were used for retention time normalization. After removing contaminant and decoy proteins hits, the sptxt file was made consensus nonabundant spectral library with the iRT retention time using spectraST.

Protein inference and quantification
The TRIC results were used for protein inference and quantification. Firstly, proteins with proteotypic peptides were considered as "uniquely identified", and the proteotypic peptides accounted for about 90% of all identified peptides (Table S1). Secondly, the peptides mapped to more than one protein entry were handled as followed. 1. The peptides shared with the proteins with proteotypic peptides were excluded for protein inference and quantification. 2. The peptides without evidence of unique protein mapping were considered as "from one protein representing the gene locus and expressed as the alphabetically first entry of the protein database (gene locus identification)" (30). To generate complete quantitative matrix of IP dataset, peptides only identified in all biological replicates of at least one timepoint were kept for extraction of quantitative information. Peptide intensities were directly from TRIC output results, where peptide intensities were calculated by summing the top five most intense fragment ion peak areas. In each dataset, all identified peptides from the specific protein were ranked by the average intensity in all runs. Subsequently, top 3 most intense peptides of the specific protein were selected and sum of these three peptides' intensities represented the protein intensity in each run. Where <3 peptides were detected, the available peak groups were summed. We found that a large number of missing values of protein quantitation were detected when the comprehensive library was used. This issue was not caused by stochastic precursor ion selection as in shotgun MS because SWATH-MS recorded all precursor ions, nor by false positive identification as 1% protein FDR at global level was applied. As shown in Figure 3C, the peptide peak group of TIRAP was unambiguously detected at the 30min run in TRAF6 dataset, but not at the 0min run. This result demonstrates that the peptide was present in the 30min run but not in the 0min run. Therefore, the peptide intensity in 0 min run was calculated as zero, and "zero" values were obviously inconvenient for downstream bioinformatics processing. Instead of missing value imputation, we used "background intensity strategy" to address this issue. The background intensity strategy was performed as followed. 1. If one peptide was detected in run A but not in run B, the RT (retention time) of the peptide in run A was used for location of the peptide in run B. The peptide RT in run A was transformed in iRT value, which was also considered as iRT value in run B. iRT in run B was then transformed into the actual RT in run B. Considering this RT maybe not precise, RT window of the peptide in B run was taken. The RT window of the peptide in B run was calculated by extending 10 min at the center of RT value. 2. In the RT window of the peptide in B run, the product ion mz intensities were extracted and summed at each cycle since no peaks were detected across the window. The summed mz intensity values were ranked and the median value was taken.
Because there was no peak in the retention time range, the summed intensity of base line was significantly lower than other peptide peak intensities. Thus we compared intensities of 100 peaks with the summed intensities of their base line in a 20 min retention time window. We found the ratios of "peak intensity" to "the summed intensity of base line" was about 5.5-6.0, and we took the average number "5.7". The median value was multiplied with 5.7 and considered as background intensity.

Targeted analysis of phopshoproteomics SWATH-MS data using Peakview software
The spectray library was generated from DDA and pseudo-spectra files, and converted to Peakview-compatible file. The parameters of Peakveiw software were set as followed. "Number of transitions per peptide" was "6". "False Discovery Rate" was "1%". "XIC Extraction window" was "10min", and "XIC width (Da)" was "0.05".

The heatmap generation and differential expression analysis
The protein intensities were input into Perseus software (31). The Hierarchical clustering and differential expression analysis were conducted with default settings.

Manual inspection of peptide XICs
The XICs of peptides of significantly changed proteins were manually inspected, and several rules were used to remove falsely identified peptides by the software. Firstly, we removed the peptides whose XICs had poor quality. Secondly, we removed the proteins whose peptide XICs showed different quantitative trends in biological replicates over LPS treatment. Thirdly, the peptides with no obvious peak groups were removed (Supplementary file 1).

MS Data availability
All SWATH MS raw data, all spectral libraries and XICs figures have been made available through the PeptideAtlas database (www.peptideatlas.org/PASS/PASS01283).

Western Blot
Proteins of the cell lysates or immunoprecipitants were separated by SDS-PAGE and then transferred onto polyvinylidene difluoride membranes. The membranes were blocked with 5% bovine serum albumin at room temperature for 1 hour and then incubated with primary antibodies overnight at 4 °C, followed by incubation with secondary antibodies for 1 hour at room temperature. The luminescent signals of immunoblotting were analyzed using an ImageQuant LAS 4000 Scanner (GE Healthcare). Antibodies for Phospho-p65(3033S), Phospho-p38(9216S), Phospho- were purchased from Cell Signaling Technology; Antibodies for NEMO(18474-1-AP) and IRAK1(10478-2-AP) were purchased from Proteintech Group Inc. Antibody for TRAF6(sc-7221) was purchased from Santa Cruz Biotechnology, Inc.
Luciferase reporter assay HEK293T cells were transiently transfected with the NF-B firefly luciferase reporter plasmid, pRL-TK-Renilla-luciferase plasmid and the indicated expression constructs. After 24hr, the luciferase activity in the total cell lysate was measured with the Dualluciferase reporter assay system (Promega). Firefly luciferase activity was normalized to Renilla luciferase activity.

RNA extraction and quantitative RT-PCR
Total RNA was isolated from cells using RNAiso Plus (Takara) according to the manufacturer's instructions. cDNA was prepared with M-MLV reverse transcriptase and oligo-dT primers. Quantitative RT-PCR was performed in a CFX96 RealTime System (Bio-Rad) by using SYBR Green reagent along with gene-specific primers. All the results were analyzed by relative quantification by normalizing to the internal control GAPDH RNA level. Primer sequences are available on request.

Results
An AP-SWATH-MS workflow using combination of the pre-built spectral library and the internal library generated directly from SWATH-MS data Affinity purification (AP)-SWATH workflow has been utilized to investigate dynamic protein interactions in signaling pathways (32-34). A spectral library pre-built from existing shotgun MS data was commonly employed for interpretation of SWATH-MS data in those studies. Meanwhile, an internal library directly generated from SWATH-MS data can also be used. The former often identifies more peptides but the latter can reveal novel peptides not detected by pre-built spectral library(35). This is most likely due to that pre-built library could be more complete but the internal library contains more specific peptide information derived directly from given SWATH-MS analysis. Therefore, more complete interpretation of SWATH-MS data can be achieved by combining the comprehensive pre-built library and the internal library built from SWATH-MS data. As TRAF6 is an effector essential for LPS signal transduction, we used measurement of dynamic interactions between TRAF6 and other proteins in LPS-stimulated macrophages as an example to describe the AP-SWATH workflow. To validate the crucial role of TRAF6 in LPS signaling pathway, we first knocked out Traf6 in macrophage cell line RAW 264.7. As anticipated, deletion of TRAF6 impaired LPS responses as indicated by poor activation of MAP kinases (p38, JNK and ERK) and NF-B ( Figure S1A). To purify endogenous TRAF6 complex, we generated Flagknock-in TRAF6 RAW 264.7 cell line. Because of the similar LPS responsiveness between wildtype (WT) RAW264.7 and Flag-knock-in TRAF6 cells ( Figure S1B and S1F), we used SWATH-MS to quantify proteins in anti-Flag immunoprecipitation samples from the Flag-KI TRAF6 cell line. The cells were treated with LPS for 10 different time periods in biological quadruplicates and TRAF6 complexes were then purified by anti-Flag immunoprecipitation and subjected to SWATH-MS analysis ( Figure 1A).
We have previously built a spectral library (pre-built library) containing 109,323 peptide sequences which were assigned to 8,599 mouse proteins by extensive fractionation at protein and peptide level of murine cell line L929 (36). Untargeted analysis software, Group-DIA (36), was used to generate internal library directly from TRAF6 SWATH-MS data ( Figure 1B). The internal library or the pre-built library was input into OpenSWATH (13) for detection of peptides. At 1% protein global FDR level, 16,418 peptide sequences were detected with pre-built library, whereas 9,423 peptides were identified with internal library. Among the identified peptides, 6,720 peptides were detected by both libraries (Figure 1C and Table S1). Indeed, pre-built librarybased strategy detected more peptides than internal library-based approach, and the internal library-based strategy did detect peptides that were not available in the prebuilt library. Importantly, TRAF1, TIRAP and TANK, which are known to play a role in LPS signaling pathway, were identified by the internal library, but not by pre-built library (Table S1). To comprehensively analyze TRAF6 SWATH-MS data, we combined these two libraries and detected about double the amount of peptides and 50% more proteins in comparison with those detected by the internal library alone at 1% protein FDR ( Figure 1C and Table S1). With the combined library, 16,519 unique proteotypic peptide sequences assigned to 2,921 proteins were identified. To evaluate whether we achieved an in-depth exploration of TRAF6 SWATH-MS data, we searched literatures and found 18 of the identified proteins were previously reported to be involved in LPS-induced cell activation (1). Some of them have not been reported to be associated with TRAF6. Of these 18 proteins, 15 proteins were identified with more than one peptide (Table S1). Thus, using combination of comprehensive pre-built library and internal library in AP-SWATH workflow provides in-depth exploration into protein interactions. In addition to analyzing the TRAF6 complex, we also analyzed MyD88 and NEMO complexes using AP-SWATH workflow as described in Figure 1. Flag-knock-in MyD88 and 3xFlag-NEMO-reconstituted Nemo -/-RAW 264.7 cell lines behaved similarly to WT RAW 264.7 cells in response to LPS stimulation ( Figure S1C, D (Figure  2A, S2A and S3A). To further determine quantitative reproducibility in biological replicates, we computed the coefficient of variation (CV) at each time point. For all LPS treatment time points in three datasets, median CVs of log2-transformed protein abundance were below 10% (Figure 2B, S2B and S3B). In TRAF6 dataset, the minimal CV was 5.16 ± 3.66 % (expressed as median ± standard deviation) at the first time point, and the maximal CV was 7.768 ± 4.082 % at the fourth time point. Taken together, this demonstrates excellent reproducibility of the entire experiment.
To evaluate the dynamic range of the MyD88, TRAF6 and NEMO interactome, we elected to estimate the range of their abundances via the 'TOP3 peptides' approach. Top3 label-free quantification method shows high sensitivity and accuracy, which also is the preferred method in the absence of reference protein measurements (38, 39). We plotted the estimated log10 protein abundances for the interacting proteins ordered from high to low abundance. The interacting protein abundances in TRAF6 dataset were spread over about 4.5-5 orders of magnitude in each time point (Figure 2C, S2C and  S3C). The wide dynamic range of protein abundances showed the effective sensitivity of our AP-SWATH workflow.

Identification of high-confidence interactors of MyD88, TRAF6 and NEMO
To obtain an overview of interacting proteins, we performed hierarchical clustering of MyD88, TRAF6 and NEMO interactions across different timepoints of LPS treatment ( Figure 3A, S4A and S5A). Majority of detected proteins in TRAF6 immunoprecipitates did not change in abundance across 10 time points of LPS treatment, implying that these proteins either consistently bound to TRAF6 or were nonspecifically pulled down ( Figure 3A). Some proteins in TRAF6 dataset showed dynamic changes and were clustered tightly into two clusters ( Figure 3A). Cluster 1 included the proteins recruited to TRAF6 complex at 30-60 min after LPS treatment. TNIP1, TNIP2, TIRAP, IRAK1, IRAK2, IRAK3, IRAK4, TRAF2 and MyD88 were in this cluster. Cluster 2 contained the proteins whose intensities peaked at 120-360 min after LPS treatment. Those were MRP, OASL1, RSAD2, TNAP3, CMPK2, IRG1, TRAF1, IL1B, and CCL5. Similarly, clusters were identified from MyD88 and NEMO datasets (Figure S4A and S5A). To further analyze the interactors of the bait proteins, we conducted a differential expression analysis individually at each time point ( Figure  3B, S4B and S5B). Based on fold of changes (|log2(fold change)|>1) and P value (-log10(P-value)>1.5) at all time points, we deduced 91, 102 and 54 proteins from 2,641, 2,991 and 2,657 quantified proteins in MyD88, TRAF6 and NEMO datasets respectively as the dynamic interactors of baits (Table S2). To further confirm these proteins were really increased or decreased in the given IPs at different time points of LPS treatment, we manually checked the XICs of all peptides of the proteins ( Figure  3C and Supplementary file 1), and confirmed 24, 25 and 12 protein abundances increased in MyD88, TRAF6 and NEMO IP complexes respectively (Table S2).
Western blot results were also highly consistent with quantitative SWATH-MS data ( Figure 3D, S4C and S5C). Of these proteins, some were exclusively detected in the bait protein complexes at LPS treatment for a long time, such as at 240-360 min. To determine whether the detection of these proteins were associated with induction of their expression, we measured the mRNA expression levels of some of the proteins and found they were markedly upregulated by LPS treatment (Figure S6A), implying that detection of these proteins in the bait protein complexes at 240-360 min of LPS treatment could be resulted from the increase of these protein levels rather than specific recruitment to bait proteins. To obtain a protein list of background binders of M2 beads, we generated a RAW 264.7 cell line expressing 3xFlag-vector and conducted immunoprecipitations using M2 beads in this cell line. The data of ten timepoints of LPS treatment in biological duplicates were collected. LPS-dependent background binders of M2 beads were determined by differential expression analysis (Figure S6B), and were removed from the interactors of the bait protein IPs. Ultimately, we obtained 20, 24 and 12 high-confidence MyD88, TRAF6 and NEMO interactors, respectively ( Table S2). As expected, majority of these proteins were well-established to be involved in LPS signaling pathway, but we identified an unknown protein, WRNIP1, which was recruited to MyD88 in an LPS-dependent way. WRNIP1 was recently reported to be involved in RIG-I-mediated antiviral signaling pathway (40), but the role of WRNIP1 in LPS signaling pathway remains undefined.

Dynamics of the recruitment of different IRAKs in the initial LPS signaling cascade
Dynamic profiles of high-confidence interactors were shown in the Figure 3E. In MyD88 and TRAF6 complexes, the highest abundance protein was IRAK2. The relative ratio of IRAK4:IRAK2 in MyD88 complex were consistent with that in TRAF6 complex (Table S3). However, the abundance of IRAK3 was 100 time lower than that of IRAK2 in MyD88 complex at 45-60 min, while its abundance was about one tenth of IRAK2 abundance at 30-60 min in TRAF6 complex.  Figure 3F). IRAK2-NEMO association has not been documented, and biological implication of the interaction requires more research.

The Stoichiometry of MyD88 and TRAF6 complexes revealed a signal amplification mechanism for MyD88-dependent TLR4 signaling
Besides providing relative quantification of proteins across different samples, label-free SWATH-MS intensity can be used to estimate subunit stoichiometry in one sample (16). It should be noted that "Top3 peptide approach" label-free absolute quantification was employed in this study, which was not as accurate as stable-isotope labeled peptide methods.
We first determined the stoichiometry of TIRAP and MyD88 in TRAF6 dataset. MyD88 and TIRAP both harbor TIR (Toll-Interleukin-I receptor) domains, which can associate directly with TIR domains of TLR (Toll-Like Receptors). Indeed, MyD88 was reported to directly interact with the TIR domains of TLR3, 7 and 9, but the binding of MyD88 to TLR4 requires TIRAP (46, 47). TRAF6 was recruited to MyD88 complex at 30min of LPS treatment and dissociated at 60min (Figure 3E), therefore we analyzed the TIRAP:MyD88 stoichiometry in TRAF6 dataset during 30-60min after LPS treatment( Table S3). As shown in Figure 4A, the stoichiometry of TIRAP: MyD88 was about 1:5. This data suggested that TIRAP binds to TLR4, thereby allowing multiple MyD88 recruitment to TIRAP, which subsequently leads to signal cascade amplification. Next, we attempted to investigate the stoichiometry of IRAK family proteins. MyD88:IRAK4 in TRAF6 dataset was about 3:1. The stoichiometry of IRAK1-IRAK4 and IRAK2-IRAK4 complexes was calculated in MyD88 dataset where IRAK3 abundance was about 100-time lower than that of the other three proteins ( Table S3). Figure 4B, the ratio of IRAK1:IRAK4 changed significantly during 15-90 min after LPS stimulation. In contrast, the stoichiometry of IRAK2-IRAK4 complex showed more stability and remained about "3:1". In TRAF6 complex, the IRAK2-IRAK4 stoichiometry was about "5:1", slightly larger than "3:1" in MyD88 complex ( Figure 4A). Taking into account that MyD88 acts upstream of IRAK proteins and TRAF6 functions downstream, IRAK protein complex might have disassembled in TRAF6 complex. Therefore, the 3:1 stoichiometry of IRAK2-IRAK4 should be more reasonable. A previous study described a "1:1" of IRAK2:IRAK4 based on the structure of death domains of IRAK2 and IRAK4 (48), but our data demonstrated that it is "3:1" in MyD88 complex. These observations suggested IRAK2 was tightly associated with IRAK4 during recruitment to and dissociation from MyD88 protein.

As shown in
On the basis of these results, we proposed a signal amplification model in MyD88mediated LPS pathway ( Figure 4C). Although "Top3" quantification method has been reported to be about two-fold error (38), some conclusions are not affected by the error. Specifically, ratio of "MyD88:TIRAP" was more than "1", suggesting the signal was amplified in TIRAP-MyD88 complex. Additionally, ratio of "IRAK2-IRAK4" in TRAF6 and MyD88 complex should be greater than "1", which also indicates potential signal amplification in this complex.

LPS-induced dynamic phosphorylation of high-confidence interactors of TRAF6, NEMO and MyD88
In addition to identification of key components in MyD88, TRAF6 and NEMO complexes, we also identified and quantified the phosphorylation events in these complexes. Phosphopeptides in IP samples were enriched using IMAC, followed by shotgun MS and SWATH-MS analysis separately. The phosphopeptides identified by these two methods were combined. Filtered at 1% peptide level and localization score of 0.7, 1,893, 1,269 and 1,122 distinct phosphopeptides were identified in MyD88, TRAF6 and NEMO complexes, respectively (Table S4). Among these phosphopeptides, 1,006, 1,022 and 935 phosphopeptides can be quantified across the timepoints of LPS treatment in MyD88, TRAF6 and NEMO complexes, respectively (Table S4). We then focused on quantified phosphopeptides of aforementioned high-confidence interactors.
To remove weak intensity phosphopeptides, we manually checked and extracted their XICs and excluded ambiguously localized phosphopeptides (Table S4). We identified 103 phosphosites in these high-confidence interactors and 41 of them have not been reported (Table S4). We first compared the phosphopetide profiles with their corresponding non-phophopeptide profiles. To obtain the stoichiometry of the phosphorylation sites, temporal phosphorylation levels of each site were subsequently compared with their protein levels ( Figure 5A). S389 and S412 on TAK1, S136, T140, S160, S171, S376, S436 and S616 on IRAK2, and S186 on IRAK4 were quantified across ten LPS timepoints in MyD88 dataset. S412 on TAK1 and S186 on IRAK4 showed a similar pattern with their proteins, suggesting the changed phosphorylation probably were attributed to varied protein levels. S389 on TAK1 exhibited distinct profile compared with the protein profile, suggesting S389 was phosphorylated in an LPS-dependent way. Of seven quantified phosphosites on IRAK2, 5 phosphosites showed a similar pattern with IRAK2 protein. S160 and S171 phosphorylation in IRAK2 apparently showed different profiles from that of their protein, suggesting the phosphorylation on these sites occurred after IRAK2 recruitment to MyD88 ( Figure  5B). In TRAF6 complex, phosphorylation of S525 on TRAF6 was LPS-induced, while S291 phosphorylation was LPS-independent. Notably, S185 and S188 phosphorylation on IRAK1 in TRAF6 complex peaked at 15min and decreased sharply, whereas IRAK1 protein abundance in TRAF 6 complex increased from 15min and stayed at the plateau from 30-60min. This discrepancy between IRAK1 protein and S185 and S188 phosphorylation of IRAK1 suggested the two Serine amino acids on IRAK1 in TRAF6 complex were not phosphorylated. Given that IRAK1 was first recruited to MyD88 and then shifted to TRAF6, S185 and S188 phosphorylation on IRAK1 in MyD88 complex were probably dephosphorylated followed by translocation to TRAF6 (Figure 5C). In NEMO complex, phosphorylation of S148, but not S380, was regulated in an LPSdependent manner. Similarly, S672 and S751 on IKKB, S389 on TAK1, and S381 and S533 on TNAP3 (also named as A20) were phosphorylated in an LPS-dependent manner. S381 on A20 was reported to be phosphorylated by IKK , which enhanced the function of A20 to down-regulate pro-inflammatory cytokines (49). Since TRAF6, IRAK1, IKK , NEMO and TAK1 were involved in NF-B activation, we sought to examine whether the identified phosphosites on them affected the ability to induce NF-B activation. To this end, we transfected wildtype and phosphorylation site mutant constructs into 293T cells and determined their ability to induce NF-B luciferase reporter activity. Disruption of most phosphorylation sites did not affect these protein function to induce NF-B activity ( Figure S7). However, we did find that substitution of S148, but not 380, to alanine impaired NEMO overexpression induced NF-B reporter gene expression ( Figure 5D). This result showed that S148 phosphorylation on NEMO most likely played a role in NF-B activation.

Overview of LPS-induced signaling events of protein interactions and phosphorylation
To have an overview of LPS-induced signaling events, we overlaid the data of temporal interaction of signaling molecules and the phosphoproteomics data onto a literaturederived LPS signaling pathway and summarized them in Figure 6. Almost all known key players in LPS signaling pathway are in the identified proteins using our AP-SWATH workflow. We also identified a number of high-confidence LPS-induced phosphorylation sites, in which about half were not reported before (phosphosites in red in Figure 6). This study may serve as a resource for in-depth understanding of LPS signaling pathway.

Discussion
In this study, we quantitatively measured time-resolved protein interactions and phosphorylation in LPS signaling pathway using SWATH-MS. This approach is critically dependent on the peptide assay library that is used to identify and quantify them. The peptide libraries are usually constructed by DDA datasets. In this study, we used a comprehensive spectral library containing 8,599 murine proteins to analyze the IP SWATH-MS data. Even with this library, not all peptides in SWATH-MS data from IP samples were identified. Targeted analysis using an internal library directly from SWATH-MS led to the identification of extra key proteins that are involved in LPS signaling pathway. This suggested that in-depth exploration of IP SWATH-MS necessitates a combination of the external library and the internal library. Whether this combination is required for the full interpretation of SWATH-MS data from complex samples such as cell lysates or tissues needs further research. We quantified about 2,500-3,000 proteins in TRAF6, MyD88 and NEMO IP datasets across about 21-40 IP samples. Of these proteins, the majority were probably background binders of M2 beads as their abundances did not change during LPS treatment. To identify LPS-dependent background binders, we performed immunoprecipitations in 3xFlag-vector-expressing cells at a serial of time points of LPS stimulation. We detected 943 proteins in IP samples using 3xFlag-vector-expressing cells (Table S2). These proteins contained "normal background binders" and "LPSdependent background binders", which were both excluded from interactors of the bait proteins. Importantly, the number of identified proteins in negative-control experiments were apparently less than those identified in MyD88, TRAF6 and NEMO datasets. This suggested that pull-downs in negative-control cells actually represent only part of the background proteins in the IP using the bait protein.
We selected the proteins whose abundance in the IPs significantly changed upon LPS treatment. Selection based on this criterion has successfully retrieved high-confidence interactors of the bait proteins, and the majority of the interactors have been proved in other studies. Nevertheless, there are other reasons except for LPS-induced recruitment that can lead to our detection of changes in protein abundances in the complexes. One is that the amounts of a number of proteins were increased after LPS stimulation, which could then increase the background binding of these proteins. Another reason is that proteins would non-specifically bind to M2 beads if they formed aggregates after LPS stimulation. Although the formation of large protein aggregates was not reported in LPS-treated RAW 264.7 cells, we had employed high speed centrifugation to remove the insoluble large aggregates.
It should be pointed out that the criterion will miss the interactors that are always binding to bait proteins. Indeed, IKKA and IKKB that are well-established binders of NEMO are not selected in this study. However, the percent of these proteins of all interactors (2/34=5.88%) is relatively low. In addition, we showed that automated analysis by OpenSWATH requires manual inspection. Peak groups may be wrongly assigned or the XICs of correctly assigned peptides may be of poor quality, which are the reasons for us to remove the peptide results by OpenSWATH analysis. Furthermore, we removed the proteins whose peptide XICs showed different quantitative trends in biological replicates over LPS treatment (Supplementary file 1). Based on the temporal interactor profiles, we demonstrated differential recruiting pattern of IRAK family proteins in MyD88, TRAF6 and NEMO complexes, and proposed a model of assembling and de-assembling of IRAKs in the signaling complexes of LPS (Figure 3F and 4C). Our data suggest IRAK1 recruitment and phosphorylation was not necessarily related at least in RAW 264.7 cells (Figure 5C). It is unexpected that IRAK2 associated with NEMO at 15min-120min after LPS treatment. Since it is well known that IRAK2 and IRAK4 formed a tight complex in Myddosome, the IRAK2 in NEMO complex may perform unknown functions. Stoichiometry of a given complex was often deduced from the structure data. However, the stoichiometry in biological reactions could be different or dynamically changeable. SWATH-MS enabled the calculation of relative protein quantitation in one sample, and thus allowed us to determine the ratio of MyD88, TIRAP, IRAK2 and IRAK4 in Myddosome.
We also identified a large number of high-confidence phosphosites in the interactors in TRAF6, MyD88 and NEMO complexes. Many of them were not reported before.
Especially, 12 phosphosites on IRAK2 were identified, only one of which is included in the Phosphosite.org database. We also showed that one of the identified phosphosites on NEMO affected the ability of NEMO to induce NF-B activation. Finally, we believe the data of our AP-SWATH-MS analysis will offer an important resource for further LPS pathway research. Group-DIA (data-independent acquisition) was used to construct pseudo-spectra from SWATH-MS data, and an internal library was made by these pseudo-spectra. The internal library and a pre-made external library of murine cell line were used either independently or in combination (combined library) in OpenSWATH analysis of SWATH-MS data. The final results were filtered at 1% global protein level. C. Comparison of the numbers of peptides and proteins detected in TRAF6 complexes (TRAF6 dataset) by using internal, external or combined library.