Discovery of Mouse Spleen Signaling Responses to Anthrax using Label-Free Quantitative Phosphoproteomics via Mass Spectrometry*

Inhalational anthrax is caused by spores of the bacterium Bacillus anthracis (B. anthracis), and is an extremely dangerous disease that can kill unvaccinated victims within 2 weeks. Modern antibiotic-based therapy can increase the survival rate to ∼50%, but only if administered presymptomatically (within 24–48 h of exposure). To discover host signaling responses to presymptomatic anthrax, label-free quantitative phosphoproteomics via liquid chromatography coupled to mass spectrometry was used to compare spleens from uninfected and spore-challenged mice over a 72 h time-course. Spleen proteins were denatured using urea, reduced using dithiothreitol, alkylated using iodoacetamide, and digested into peptides using trypsin, and the resulting phosphopeptides were enriched using titanium dioxide solid-phase extraction and analyzed by nano-liquid chromatography-Linear Trap Quadrupole-Orbitrap-MS(/MS). The fragment ion spectra were processed using DeconMSn and searched using both Mascot and SEQUEST resulting in 252,626 confident identifications of 6248 phosphopeptides (corresponding to 5782 phosphorylation sites). The precursor ion spectra were deisotoped using Decon2LS and aligned using MultiAlign resulting in the confident quantitation of 3265 of the identified phosphopeptides. ANOVAs were used to produce a q-value ranked list of host signaling responses. Late-stage (48–72 h postchallenge) Sterne strain (lethal) infections resulted in global alterations to the spleen phosphoproteome. In contrast, ΔSterne strain (asymptomatic; missing the anthrax toxin) infections resulted in 188 (5.8%) significantly altered (q<0.05) phosphopeptides. Twenty-six highly tentative phosphorylation responses to early-stage (24 h postchallenge) anthrax were discovered (q<0.5), and ten of these originated from eight proteins that have known roles in the host immune response. These tentative early-anthrax host response signaling events within mouse spleens may translate into presymptomatic diagnostic biomarkers of human anthrax detectable within circulating immune cells, and could aid in the identification of pathogenic mechanisms and therapeutic targets.

Anthrax is a life-threatening infectious disease that affects humans and many animal species, and Bacillus anthracis (B. anthracis), the causative agent of anthrax, is a Gram-positive, rod-shaped, nonmotile, endospore-forming bacterium (1,2). During the 20th century, effective human and animal vaccines almost completely eradicated anthrax as a naturally occurring disease from the industrialized world. However, anthrax remains a major concern because several nations stockpiled weaponized B. anthracis spores, and because B. anthracis has been used as a bioterror weapon (e.g. the 2001 anthrax letters attack in the United States) (3). As a result, B. anthracis is classified as a Category A Select Agent by the Centers for Disease Control and Prevention (some attenuated strains are excluded).
The major routes of infection occur via inhalation, cutaneous abrasion, and ingestion of B. anthracis spores, which can survive for years in the environment. Inhalational anthrax is the most severe form of the disease, and in the past it was almost invariably lethal within ϳ2 weeks. Fortunately, modern antibiotic-based therapies have significantly reduced anthrax mortality. As a result, 55% of the inhalational anthrax victims of the 2001 U. S. attack survived. However, antibiotics are effective only against proliferating B. anthracis, and therefore improve patient outcome only if administered presymptomatically (within 24 -48 h of exposure), when the level of anthrax toxemia is low.
The major toxic factors of B. anthracis are the lethal and edema toxins encoded by the pXO1 plasmid, although other bacterial factors might additionally play important pathogenic roles. Lethal toxin is a proteolytic inhibitor of mitogen-activated protein kinase kinase, and edema toxin is a bacterial adenylate cyclase capable of increasing the levels of intracellular cyclic AMP. Both of these enzymatic activities have a profound effect on cell function and viability. The antiphagocytic capsule produced by proteins encoded by the pXO2 plasmid is also an important virulence factor, but it has no inherent toxicity. Elimination of these two plasmids greatly reduces the virulence of B. anthracis. For example, the Sterne strain (pXO1 ϩ , pXO2 Ϫ ) is an important vaccine strain that causes a mild infection in most animals, including humans. In DBA/2J mice, however, Sterne B. anthracis causes a highly lethal infection. Consequently, Sterne-challenged DBA/2J mice have become a valuable model of human anthrax. In contrast, the plasmidless ⌬Sterne strain (pXO1 Ϫ , pXO2 Ϫ ) causes an asymptomatic, abortive infection in DBA/2J mice.
The molecular events that result in the high mortality of anthrax victims is still an active area of research, and there is no known cell model or biomarker that is predictive of patient outcome. After the 2001 attacks, police in the USA have investigated roughly 500 "white powder" hoaxes per year (4), and prophylactically prescribing powerful antibiotics can result in uninfected patients developing serious adverse effects. In a large-scale anthrax hoax or in an actual attack, first responders and hospitals would have great difficulty distinguishing uninfected and presymptomatic patients. Therefore, efforts in our laboratory are focused on identifying diagnostic host response biomarkers of presymptomatic anthrax, and we recently identified candidate anthrax biomarkers in the low-molecular-mass serum proteome of Sterne and ⌬Sterne spore-challenged DBA/2J mice using liquid chromatography coupled to (tandem) mass spectrometry (LC-MS[/MS]) 1 (5). Proteomics (6 -12), including phosphoproteomics (13)(14)(15)(16), has greatly benefited from recent advances in LC-MS(/MS), though tentative biomarkers still require validation by orthogonal experimentation.
In this investigation, host signaling responses to anthrax were discovered using label-free quantitative phosphoproteomics to compare spleens from uninfected and spore-challenged mice. The spleen was studied because of its role in the immune system, because it is among the first of the organs to display the presence of disseminating B. anthracis in the murine model, and because of the difficulty of reproducibly preparing sufficient masses of circulating immune cells directly from the blood of spore-challenged mice. Synergistically with its role as a hematopoietic organ (i.e. as a site of immune cell maturation, activation, and proliferation), the spleen also functions as a secondary lymphoid organ both by directly filtering the blood to capture pathogens and antigens, and by transiently concentrating phagocytes and antigen presenting cells (e.g. dendritic cells, granulocytes, macrophages, and monocytes) and antigen-specific B-and T-lymphocytes, resulting in an increase in the efficiency of the host immune response (17)(18)(19)(20). The splenic marginated pool of immune cells is in dynamic equilibrium with the circulating immune cells, and includes 10 -15% of the body's B lymphocytes (20,21), 25% of the body's T lymphocytes (20,21), 21-34% of the body's granulocytes (22,23), and 50% of the body's monocytes (24) (this reserve of monocytes can exit the spleen and differentiate to become dendritic cells and macrophages). The strategy of using comparative proteomics to study altered tissues to generate lists of candidate biomarkers has been employed previously, primarily because of the difficulty of discovering biomarkers directly from the blood (25)(26)(27)(28)(29)(30)(31)(32). Host response signaling events within mouse spleens may translate into presymptomatic diagnostic biomarkers of human anthrax detectable within circulating immune cells, and could aid in the identification of pathogenic mechanisms and therapeutic targets. In this study, thousands of mouse spleen phosphopeptides were confidently identified and quantitated, demonstrating that label-free quantitative phosphoproteomics is a viable alternative to stable isotope labeling (e.g. 18 O, isobaric tag for relative and absolute quantification, Tandem Mass Tags, and stable isotope labeling with amino acids in cell culture).

EXPERIMENTAL PROCEDURES
Mouse Infection and Tissue Homogenization-The George Mason University Institutional Animal Care and Use Committee and the Biocon Animal Care and Use Committee/Institutional Review Board (Biocon Inc., Rockville, MD) approved all of the animal experimentation protocols. Using BSL-2 safety protocols at Biocon Inc., mice (Mus musculus strain DBA/2J, female, 8 weeks old, The Jackson Laboratory, Bar Harbor, ME) were either not injected at all (five mice) or were intraperitoneally injected with 100 l of sterile H 2 O (15 mice), or 100 l ⌬Sterne (pXO1 Ϫ , pXO2 Ϫ ) B. anthracis spores in sterile H 2 O (5 ϫ 10 6 spores/100 l; 15 mice), or 100 l 34F2 Sterne (pXO1 ϩ , pXO2 Ϫ ) B. anthracis spores in sterile H 2 O (5 ϫ 10 6 spores/100 l; 25 mice). The B. anthracis Sterne 34F2 and ⌬Sterne strains (33) and the spore preparation protocol (34) have been described previously. At t ϭ 0 h, the five uninjected mice were prepared by terminal bleed from the orbital sinus using a glass Pasteur pipette, followed by cervical dislocation and dissection of the spleen, which was then placed into a 2-ml cryogenic vial and immediately snap frozen in liquid nitrogen. At 24, 48, and 72 h postchallenge, five mice from each of the remaining experimental conditions (H 2 O, ⌬Sterne, and Sterne) were prepared using the same dissection protocol. Fourteen of the Sterne-challenged mice were found dead in their cages prior to 72 h postchallenge (Sterne-infected mice have a high mortality rate at this time point) and were discarded, so only a single 72 h post-Sternechallenge mouse spleen was isolated. The spleen samples were organized into five sample blocks by placing one biological replicate from each experimental condition into each block, and then randomizing the order of the samples within each block. The samples were processed in this (random) order within these blocks to reduce any possible systematic biases (35). Each sample was homogenized by manual disruption into 600 l of freshly prepared Lysis Buffer (50 mM tris-HCl pH 8, 8 M urea), and BCA protein concentration assays (Thermo Fisher Scientific Inc., Waltham, MA) were performed.
Sample Digestion and Phosphopeptide Enrichment-Two milligrams (protein mass) of each spleen homogenate (ϳ30 -50% of a mouse spleen) was diluted with Lysis Buffer to a final volume of 370 l, and 200 ng of bovine ␤-casein was added to each sample to serve as a phosphoprotein standard. Dithiothreitol was added to a final concentration of 10 mM, and the samples were incubated at 60°C for 30 min to reduce disulfide bridges. The sample proteins were alkylated by adding iodoacetamide to a final concentration of 50 mM and incubating the samples at room temperature for 20 min in darkness. The samples were diluted with H 2 O and 500 mM NH 4 HCO 3 pH 9 such that the final urea concentration was 1 M and the final NH 4 HCO 3 concentration was 100 mM. The sample proteins were digested into peptides by adding sequencing grade modified trypsin (1:100 [protein w:w] trypsin:sample; Promega Corp., Madison, WI) and incubating the samples at 37°C for 21.5 h.
The samples were acidified in preparation for C-18 solid-phase extraction (SPE) by adding acetic acid to a final concentration of 2% v/v, and 8 pmol of phosphorylated angiotensin II was added to each sample to serve as a phosphopeptide standard. Each sample was centrifuged 40 min at 2000 ϫ g at room temperature to pellet particulates (to avoid clogging the SPE column), and then the supernatants underwent C-18 solid phase extraction using a Sep-Pak C-18 SPE 1 cc 100 mg column (Waters Corp., Milford, MA). Centrifugation in a swinging-bucket rotor was used to equilibrate the Sep-Pak columns and to apply the samples and the wash buffer (0.1% v/v formic acid). The elution buffer (0.1% v/v formic acid, 80% v/v acetonitrile [ACN]) was applied slowly using a rubber bulb to produce positive air pressure.
The Sep-Pak eluates were concentrated in a SpeedVac (no heating; Thermo Fisher Scientific Inc.) to 150 l. Next, 150 l of freshly prepared TiO 2 Loading Buffer (5% v/v trifluoroacetic acid [TFA], 80% v/v ACN, 225 mg/ml 2,5-dihydroxybenzoic acid [DHB]) was added to each sample, and the samples were microcentrifuged 10 min at 16,000 ϫ g at room temperature to pellet any particulates that might clog the TiO 2 SPE tip. Phosphopeptides from each supernatant were enriched using a Titansphere PHOS-TiO 200 l centrifuge tip containing 3 mg of TiO 2 media (GL Sciences Inc., Torrance, CA) using a protocol based on (36). Centrifugation was used to equilibrate the tips, apply the samples, apply the first wash buffer (freshly prepared 2% v/v TFA, 80% v/v ACN, 40 mg/ml DHB), the second wash buffer (2% v/v TFA, 80% v/v ACN), the first elution buffer (5% w/v [ϳ3M] NH 4 OH, to elute the phosphopeptides from the TiO 2 resin), and the second elution buffer (2% v/v TFA, 80% v/v ACN, to elute the phosphopeptides from the Empore C-8 bonded silica membrane below the TiO 2 resin inside the tips). The eluates were immediately concentrated in a SpeedVac (no heating) to 20 l to evaporate the ACN and NH 4 OH, and the samples were acidified in preparation for C-18 SPE by adding 60 l of 1% v/v acetic acid.
The samples underwent C-18 SPE using a ZipTip 10 l pipette tip (Millipore Corp., Billerica, MA) with the Sep-Pak C-18 SPE wash and elution buffers described above. During the development of this protocol, it was found that skipping the ZipTip SPE resulted in clogging of the nano-LC-electrospray ionization (ESI) column/tip, and that performing this step solved this problem. Forty l of H 2 O was added to each ZipTip eluate, and then the samples were concentrated in a SpeedVac (no heating) to 10 l to evaporate the ACN. One pmol of angiotensin I was added to each sample to serve as a peptide standard, and 7 l of 1% v/v acetic acid was added to each sample to acidify it in preparation for reversed-phase LC-MS(/MS). The samples were stored at Ϫ80°C until LC-MS(/MS).
Mass Spectrometry-Nano-LC-ESI-MS(/MS) analyses were performed using a Surveyor LC system coupled to an LTQ-Orbitrap mass spectrometer (both Thermo Fisher Scientific Inc.; Buffer A: 0.1% v/v formic acid; Buffer B: 0.1% v/v formic acid, 80% v/v ACN) (37)(38). Each nano-LC-ESI column/tip was prepared by laser-pulling a tip (Laser Based Micropipette Puller Model P-2000, Sutter Instrument Co., Novato, CA) onto a 30 cm length of 100 m I.D. coated silica capillary tubing (Polymicro Technologies LLC, Phoenix, AZ), which was then packed with 15 cm of Magic C18AQ C-18 media (5 m diameter, 200 A pores, Michrom Bioresources Inc., Auburn, CA) (note that this was just a single piece of capillary; there was no connector or frit between the column and tip). One nano-LC-ESI column/tip was used for sample blocks 1-2, and a second column/tip was used for sample blocks 3-5. Each sample was manually loaded using a pressure cell (Brechbuhler Inc., Houston, TX), the nano-LC-ESI column/tip was connected to the LC system, and the postsplit flow rate was manually calibrated to be 350 nL/min (ϳ70 bar, ϳ70 l/min presplit) using a calibrated micropipette. Note that because the sample was manually loaded in this manner, the phosphopeptides were never exposed to metal or other surfaces that might have caused sample loss (39). A 2 kV ESI voltage was applied and the ESI spray was observed. LC-MS(/MS) was performed using an instrument method that included a 120 min linear gradient (10 -40% Buffer B) and a 40 min 100% Buffer B column regeneration step. The LTQ-Orbitrap was set to perform simultaneous Orbitrap-MS (400 -1600 m/z, resolution ϭ 100,000) and shotgun CID LTQ-MS 2 (no Electron Transfer Dissociation, Higher Energy Collisional Dissociation, or Pulsed Q Collision Induced Dissociation was used) against the top eight most intense ions (top six for sample block 5). Dynamic exclusion was enabled to avoid repeatedly selecting intense ions for fragmentation (sample blocks 1 and 4 excluded from Ϫ0.1 to 1.6 m/z units about the precursor ion for 25 s; block 2 excluded from Ϫ0.1 to 1.6 m/z for 15 s; blocks 3 and 5 excluded from Ϫ1.1 to 1.6 m/z for 60 s). Multistage activation was performed during the analyses of sample blocks 4 and 5 (Neutral Loss Mass List ϭ Ϫ97.98, Ϫ48.99, Ϫ32.66 m/z units). These spectra are available for download as .RAW files at the Proteome Commons Tranche data repository at https://proteomecommons.org/(Tranche Hash: j7rN4iodsyEQbWjTaWAqMthHEdlfs CqJc/WeLuA5xsrD0Tiu7NOEAwQWpWxhpAwlϩhRPiBRzscϩK0c TFg2KipfpϩLDEAAAAAAAA9rwϭϭ).
A FASTA text file consisting of five concatenated sets of protein sequences was used to search the spectra: (i) protein and peptide standards, (ii) common contaminants (e.g. porcine trypsin, human keratin), (iii) Bacillus anthracis Sterne proteins (this dataset did not include any plasmid proteins, retrieved July  In addition to searching the normal protein sequence file, Mascot and SEQUEST were both used to search the two decoy databases. The decoy database searches were performed separately from the normal searches, and the estimated false discovery rate (FDR) was calculated using: FDR ϭ 0.5 ϫ (the number of reversed identifications ϩ the number of scrambled identifications)/(the number of normal identifications). It has been argued that this FDR estimate overestimates the true FDR, and an alternative method using concatenated (normal ϩ decoy) databases has been proposed (46). Empirically, a highly confident peptide identification will result from a highly confident first hit (the "top-ranked hit"; i.e. the highest-confidence peptide sequence that resulted from a FASTA database search against a single experimental fragment ion spectrum) and a much less confident second hit (the second-place hit). Although the concatenated database method may result in a more accurate FDR estimate, it unfortunately overestimates the confidence of the second hits (46). Therefore, the decoy database searches were performed separately, which resulted in relatively accurate second hit scores at the cost of possibly overestimating the FDR estimates.
The resulting Mascot (.dat) and SEQUEST (.out) output files were converted to tab-delimited text files using Mascot Output Parser (v2003-11-13, generously provided by Dr. Matthew Monroe, Pacific Northwest National Laboratory, Richland, WA) and Peptide File Extractor (console version, v1.1.3519.25650, http://omics.pnl.gov/software/). Microsoft Access 2003 was used to import each of the resulting text files and to perform queries in an automated fashion using macros. Systematic precursor ion mass measurement errors were determined and the search-space was reduced from 20 ppm (Ϯ10 ppm) to 7 ppm (Ϯ3.5 ppm about the systematic error) to reduce the FDR (supplemental Fig. 1) (a much more advanced algorithm to remove systematic precursor ion mass measurement errors has recently been published (47,48)). Phosphopeptide FDR Estimator (v2009 -10 Ϫ 26) (49) was used to calculate the Ambiguity Score (50) of each phosphorylation site identified by Mascot and SEQUEST, and also to perform a discriminant analysis using the SEQUEST data (supplemental Fig. 2; the resulting q-values are included as supplemental data but otherwise were not used to filter the data or to estimate FDRs of the filtered data). All of these steps were automated using Microsoft DOS command line scripting and Microsoft Access 2003 macros.
For each MS 2 spectrum produced by DeconMSn, only the topranked hit was retained (one from each of the six searches: Mascot/ SEQUEST, normal/reversed/scrambled). A two-step process was then used to filter the resulting data to reduce the FDR. First, if Mascot and SEQUEST agreed on the phosphopeptide identification (ignoring post-translational modification [PTM] localization), relatively mild data filtration criteria were employed (see supplemental Table 2 for a detailed description of all of the data filtration criteria). In the rare case that a single spectrum resulted in filter-passing Mascot and SEQUEST identifications of two different peptides, these identifications were discarded (only 515 spectra resulted in this situation). In the second step, strict data filtration criteria were employed to filter identifications made by only one of the two search engines (this resulted in a relatively small number of additional filter-passing identifications). FDRs of the filtered data were estimated using the decoy database method described above, and the data filtration criteria were designed so that Յ2.5% of the phosphopeptides were wrongly identified within each of 15 data categories (precursor ion charge states ϩ1 -ϩ5; identifications made by Mascot only, SEQUEST only, or both). Overall, this resulted in 252,626 filter-passing phosphopep-tide identifications (FDR ϭ 0.36%; an MS 2 spectrum confidently identified by both Mascot and SEQUEST still counted as a single identification) of 6248 different phosphopeptides (FDR ϭ 3.4%). This corresponded to 5782 phosphorylation sites. If phosphopeptides identified by only a single spectrum were excluded ("one-hit-wonders" are generally considered to be low-confidence), 5172 different phosphopeptides were identified (FDR ϭ 1.9%), and if Ն3 identifications were required per phosphopeptide, 4652 different phosphopeptides were identified (FDR ϭ 1.3%). Compared with using either Mascot or SEQUEST alone, using both search engines (i.e. "consensus scoring") resulted in roughly triple the number of confidently identified phosphopeptides, in agreement with previous studies (51)(52)(53)(54)(55).
Phosphopeptide Quantitation and ANOVAs-The LTQ-Orbitrap output files were deisotoped using Decon2LS (v1.0.2964.22547, http://omics.pnl.gov/software/) (56), an implementation of the THRASH algorithm (41). The resulting data were analyzed using Mul-tiAlign (v1. MultiAlign clusters that peak-matched to multiple phosphopeptides [ignoring PTM localization] were discarded unless the vast majority of the matching phosphopeptide identifications were of only one phosphopeptide (i.e. Նfourfold the number of the second-most-numerous match). The vast majority of MultiAlign clusters only matched to a single phosphopeptide because of the enormous reduction in sample complexity that resulted from the TiO 2 SPE. MultiAlign clusters with a mass range Ͼ16 ppm were discarded, "split" MultiAlign clusters (clusters confidently peak-matched to the same phosphopeptide [ignoring PTM localization]) were joined, and phosphopeptides that contained a missed trypsin cleavage site were discarded. MultiAlign was used to calculate the peak area (elution time versus intensity) of each feature, and each peak area was used as a relative measure of the abundance of the analyte. This resulted in 3265 confidently identified and quantitated phosphopeptides (FDR ϭ 1.81% estimated by peak-matching the decoy database identifications). If Ն2 identifications were required per phosphopeptide, 3006 different phosphopeptides were confidently identified and quantitated (FDR ϭ 1.3%), and if Ն3 identifications were required, 2837 different phosphopeptides were confidently identified and quantitated (FDR ϭ 0.97%).
To eliminate systematic abundance value biases between the samples (due to, for example, slightly different sample phosphopeptide masses), a central tendency normalization (58) was performed (supplemental Fig. 3). Specifically, one LC-MS(/MS) dataset was selected to be the baseline, and for each "alignee" dataset, and for each confidently identified and quantitated phosphopeptide, the alignee to baseline abundance value ratio was calculated. The median of these ratios was calculated for each alignee and used as a normalization factor. Normalization was performed by dividing all of the alignee abundance values by this normalization factor.
The abundance values were log 10 -transformed because this increased the normality of the distribution of their z-scores (supplemental Fig. 4). Then, to determine which phosphopeptide abundance values were significantly different between the experimental groups, two different types of analysis of variance (ANOVA) tests were performed for each phosphopeptide using a MATLAB script (supplemental material) (35,59). The first ANOVA was simply a parametric one-way ANOVA (missing values were excluded). The second ANOVA was a nonparametric Kruskal-Wallis one-way ANOVA by ranks (missing values were tied for the lowest rank). Each ANOVA resulted in a p value that is the probability of obtaining by random chance data at least as disproportionate as that observed. Both types of ANOVAs were employed because whereas the parametric ANOVA is a much more powerful test in general, it is unable to test for phosphopeptides that are present in some experimental conditions and absent from others. Other researchers have addressed the problem of missing values from proteomics datasets using data imputation (60), but this strategy is controversial so it was not used in this study. The resulting parametric p values were used to estimate parametric q-values using QVALUE (v1.1, http://genomics.princeton.edu/ storeylab/qvalue/index.html) (61), and the nonparametric p values were used to estimate nonparametric q-values in the same way. A q-value is defined as the proportion of false positives incurred (i.e. the FDR) when a particular p value is considered significant. A single, robust q-value for each phosphopeptide could not be calculated because the parametric and nonparametric p values were not independent. Consequently, for each phosphopeptide the minimum of the parametric and nonparametric q-values was calculated and was used for data filtration and for ranking the confidence in the phosphorylation responses, and was considered an approximate estimate of the q-value of the combined ANOVAs. Unless specified otherwise, "q" or "q-value" refers to the minimum of the parametric and nonparametric q-values. To test the efficacy of the analysis strategy of the quantitative data, selected ion chromatograms of one of the phosphorylation responses were manually inspected (supplemental Fig. 5). Hierarchical cluster analyses of the significant data were performed using Genesis (v1.7.2, http://genome.tugraz.at/) (62) with average linkage correlations determined by Euclidean distance.

RESULTS
In this investigation, spleens from uninfected and sporechallenged mice were analyzed using label-free quantitative phosphoproteomics to discover host signaling responses to anthrax. A variety of organs and biofluids were retained from each mouse, and we ultimately focused on the spleen because of its role in the immune system, and because repro-ducibly isolating sufficient masses of circulating immune cells from whole blood using BSL2 protocols was found to be very challenging. Additionally, hemolysis was caused by the advanced Sterne infections, and the resulting highly variable volumes of blood recoverable from each mouse were often insufficient for our experiments.
Initially, method development experiments were performed to design and optimize sample preparation and mass spectrometry. Notably, several different phosphopeptide enrichment procedures were tested. Microcolumns prepared using a pressure cell, frit, and capillary tubing, and packed with Titansphere 5 m ) became available recently, and they performed just as well as the Titansphere microcolumns while enabling high sample throughput (ϳ20 samples/5 h), which in turn reduced sample processing time and eliminated the reproducibility concern. Three alterations to the protocol were investigated. First, using lactic acid instead of DHB (as in (64)) didn't significantly change the number of identified phosphopeptides per LC-LTQ-Orbitrap-MS(/MS) analysis. Second, including an insoluble component of the spleen homogenate in the tryptic digestion resulted in clogging of the Sep-Pak C-18 SPE, and therefore was excluded from consideration. Lastly, simplifying the TiO 2 SPE procedure was found to have no adverse effect (described in detail in the Supplemental Protocol).
Similarly, a data analysis strategy needed to be developed to confidently discover the phosphoprotein signaling responses. To confidently identify the phosphopeptides, both Mascot and SEQUEST were used to perform database searching. This greatly mitigated the effect of the often poor fragmentation common to CID spectra of phosphopeptides, because both algorithms greatly complement one another. In fact, when the same FDR was required (of unique phosphopeptides), each algorithm alone resulted in roughly one-third of the identified phosphopeptides compared with using both search engines and requiring that both identified the same phosphopeptide (ignoring PTM localization). This is not surprising, as the use of multiple search engines (i.e. "consensus scoring") has been shown to significantly improve peptide identification confidence (51)(52)(53)(54)(55). Additionally, a short Perl script was written (supplemental material) to remove tentative neutral-loss precursor ion peaks from the fragment ion spectra .dta files prior to the Mascot and SEQUEST analyses, and this significantly increased the number of identified phospho-peptides at the same FDR. Also, the use of DeconMSn instead of extract_msn enabled operation of the LTQ-Orbitrap with monoisotopic precursor ion selection turned off, and this also increased the number of identified phosphopeptides at the same FDR.
In a preliminary comparative phosphoproteomics experiment using six mice (two uninfected, two 24 h Sterneinfected, and two 48 h Sterne-infected), spectrum counting (simply counting confident phosphopeptide identifications as a rough measure of abundance) was found to be insufficient to identify tentative phosphoprotein signaling responses. However, reanalyzing the data using MultiAlign resulted in the discovery of 97 phosphopeptides that were altered significantly within the 48 h Sterne samples compared with the other four samples (44 were more abundant and 53 were less abundant). The two phosphoproteins affected most dramatically (both were ϳ30-fold more abundant) have known roles in the immune system: interferon-␥-inducible p47 GTPase and Z-DNA binding protein (a cytosolic DNA sensor).
Based on the results of the preliminary experiments, a scaled-up discovery experiment using 60 mice was designed and performed (Fig. 1). At t ϭ 0 h, five uninjected negative control (nc) mice were dissected, 15 nc mice were injected with H 2 O, 15 mice were injected with ⌬Sterne spores, and 25 mice were injected with Sterne spores. At t ϭ 24, 48, and 72 h, five H 2 O, five ⌬Sterne, and five Sterne mice were dissected (14 of the Sterne mice died in their cages prior to t ϭ 72 h and were discarded, so the 72 h Sterne experimental condition contained only one mouse spleen). Both the lethal Sterne strain and the asymptomatic ⌬Sterne strain were studied so that the virulence-specificity of the signaling responses could be determined. Signaling responses to Sterne but not ⌬Sterne might be useful markers of lethal anthrax or of virulent bacterial infection in general. The 46 resulting spleens were homogenized, and the resulting homogenate proteins were reduced using dithiothreitol, alkylated using iodoacetamide, and digested into peptides using trypsin. Phosphopeptides were enriched using TiO 2 SPE and analyzed by nano-LC-LTQ-Orbitrap-MS(/MS). These spectra are available for download as .RAW files at the Proteome Commons Tranche data repository at https://proteomecommons.org/(Tranche Hash: j7rN4iodsyEQbWjTaWAqMthHEdlfsCqJc/WeLuA5xsrD0Tiu 7NOEAwQWpWxhpAwlϩhRPiBRzscϩK0cTFg2KipfpϩLDE AAAAAAAA9rwϭϭ).
The resulting 46 LTQ-Orbitrap datasets described above were also analyzed using Decon2LS and MultiAlign. Pairwise Pearson correlations of the log 10 -transformed MultiAlign abundance values displayed as a heatmap (Fig. 2) clearly show that the phosphoproteome profiles of the very sick (Ն48 h Sterne) mice were globally different from the other profiles but not from each other. This was not surprising as these mice were moribund, and their spleens were pale pink (the normal spleens were blood red) and smaller than normal (roughly half the volume) based on visual observation.
The MultiAlign data were peak-matched to the filter-passing phosphopeptide identifications, and additional data filtration criteria were imposed. Of the 6248 confidently identified phosphopeptides, 3265 were confidently quantitated (FDR ϭ 1.81% estimated by peak-matching the decoy database identifications). If Ն2 identifications were required per phosphopeptide, 3006 different phosphopeptides were confidently identified and quantitated (FDR ϭ 1.3%), and if Ն3 identifications were required, 2837 different phosphopeptides were confidently identified and quantitated (FDR ϭ 0.97%). These data are included as supplemental Table 4.
A global normalization of the abundance values was performed to eliminate systematic biases between the samples because of, for example, slightly different sample phosphopeptide masses (supplemental Fig. 3). Then log 10 -transformation was performed because it increased the normality of the z-score distribution of the abundance values (supplemental Fig. 4) indicating that the unadjusted values were approximately log-normal. A plot of the geometric mean ( G ) versus geometric standard deviation ( G ) of the globally normalized but otherwise unadjusted abundance values indicated that G Ͻ 2 for most of the phosphopeptides (Fig. 3). This indicates that phosphorylation abundance changes of twofold upon infection might be statistically discernable. Also, it was noted that G decreased as G increased, likely because the Or-bitrap intensity measurement is less accurate at lower precursor ion intensities.
Following globally normalizing and log 10 -transforming the abundance values, ANOVAs were used to identify phosphopeptides that were disproportionately abundant in one or more of the experimental groups (35,59). Of the 3265 confidently identified and quantitated phosphopeptides, 1173 (36%) were found to be significantly altered (q Ͻ 0.05). Note that by requiring q Ͻ 0.05, 5% of the data (59 phosphopeptides) are estimated to have been wrongly considered significant by the ANOVAs. The significantly altered phosphopeptide data were analyzed using a hierarchical cluster analysis and displayed as a heatmap (Fig. 4 Top; the corresponding data is included as supplemental Table 5). The Ն48 h Sterne samples were very different from the other samples, in agreement with Fig. 2. It should be noted that to justify performing the global normalization, it was assumed that any abundance differences between the samples would average out. This assumption was true even for the Ն48 h Sterne samples, as the global normalization successfully removed systematic biases from even these datasets (supplemental Fig. 3; note that the MultiAlign and global normalization baseline was the same sample, and that it was one of the 24 h ⌬Sterne samples). MultiAlign phosphopeptide abundance values were globally normalized, and then the geometric mean ( G ) and geometric standard deviation ( G ) were calculated for each analyte across the 46 spleen samples. Confidence intervals of log-normal data take the form G G Ϫn -G G n (e.g. the 95% confidence interval is from G G Ϫ2 to G G 2 , or from 0.5 ϫ G to 2 ϫ G if G ϭ 2 0.5 ). Phosphopeptide abundances altered twofold are potentially statistically discernable.  3265) were globally normalized, log 10 -transformed, and analyzed by ANOVAs across the seven experimental conditions. The significant data (qϽ0.05, n ϭ 1173) were z-score transformed across the 46 samples, and analyzed using a hierarchical cluster analysis. Each row represents a confidently identified phosphopeptide, each column represents a mouse spleen, and relative abundance values are indicated by color (green ϭ relatively low, red ϭ relatively high, gray ϭ missing data). Most of the phosphopeptides were altered in the Ն48 h Sterne samples, correlating strongly with anthrax pathogenesis as the Sterne-infected mice began to die at 48 h postchallenge and the mice infected with the ⌬Sterne strain (missing the anthrax toxin) were asymptomatic. Bottom Heatmap: All of the data except for the Ն48 h Sterne data were reanalyzed by ANOVAs across the five remaining experimental conditions, and the significant data (qϽ0.05, n ϭ 188) were analyzed using a hierarchical cluster analysis. The resulting heatmap depicts the mouse spleen signal transduction cascade responding to the asymptomatic ⌬Sterne infection.
Ideally, a post hoc analysis such as the Tukey-Kramer test would have been used to perform pairwise comparisons to determine which phosphopeptides were significantly different between specific pairs of experimental groups. Unfortunately, an algorithm to calculate the FDR of such analyses has yet to be discovered. Therefore, to discover which phosphopeptides were altered between the other experimental conditions (i.e. other than Ն48 h Sterne), the ANOVAs were recalculated excluding the Ն48 h Sterne data. This time, of the 3265 confidently identified and quantitated phosphopeptides, 188 (5.8%) were found to be significantly altered (q Ͻ 0.05). Note that by requiring q Ͻ 0.05, 5% of the data (nine phosphopeptides) are estimated to have been wrongly considered significant by the ANOVAs. The significantly altered phosphopeptide data were analyzed using a hierarchical cluster analysis and displayed as a heatmap (Fig. 4 Bottom; the corresponding data is included as supplemental Table 6). Almost all of these 188 phosphopeptides were disproportionately abundant within the Ն48 h ⌬Sterne samples. Clearly these phosphoproteins were a major component of the mouse spleen immune response to the asymptomatic ⌬Sterne challenge. Notably, this response was mostly absent at 24 h postchallenge, and it was decreased in intensity by 72 h postchallenge, in agreement with the abortive nature of the ⌬Sterne infection.
To discover phosphopeptides that were altered at 24 h postinfection, the ANOVAs were recalculated yet again, this time excluding the Ն48 h ⌬Sterne and Sterne data. This time, of the 3265 confidently identified and quantitated phosphopeptides, 26 (0.8%) were tentatively found to be altered (qϽ0.5). Note that by requiring q Ͻ 0.5, 50% of the data (13 phosphopeptides) are estimated to have been wrongly considered significant by the ANOVAs. The tentatively altered phosphopeptide data were analyzed using a hierarchical cluster analysis and displayed as a heatmap ( Fig. 5; the corresponding data is included as supplemental Table 7) . Almost all of the 26 phosphopeptides were affected the same way by the ⌬Sterne and Sterne strains, indicating similarity between the early host responses to the toxigenic and nontoxigenic infections. It's reasonable to question whether or not these early signaling events are specific to B. anthracis infection, or if they are not anthrax-specific but rather are common to bacterial infection in general.
Ten of the 26 phosphopeptides were derived from eight proteins that have known roles in the immune system. Two of these phosphoproteins, Pram1 and Hemgn, each had two phosphorylation sites (on two separate phosphopeptides) that were affected similarly during early-anthrax (Fig. 5), indicating that these proteins were probably genuinely altered at 24 h postchallenge. Because the phosphopeptides from each of these two proteins clustered together, it's reasonable to speculate that these proteins were multiply phosphorylated/dephosphorylated by the same kinase/phosphatase. Of the 26 phosphopeptides, the one with the lowest q-value was de-rived from interferon-␥-inducible p47 GTPase, a regulator of immunity to intracellular pathogens (Gm12250; p[parametric] ϭ 0.0000345; q[parametric] ϭ 0.11). Overall, these 26 phosphopeptides constitute a q-value-ranked list of candidate biomarkers of early-stage anthrax available for targeted validation experiments using circulating immune cells. Ideally, a set of highly-confident (q Ͻ 0.001) early-anthrax candidate biomarkers would have been discovered. However, comparative proteome profiling of highly similar samples often has the problem of a small number of slightly altered proteins, making it challenging to confidently discover the earliest biomarkers of diseases in general.
In an attempt to further refine the list of 26 phosphopeptides, the results of the three sets of ANOVAs were compared, and it was found that 15 of the phosphopeptides were discovered by all three analyses (Fig. 6). Thus, these 15 phosphopeptides were mildly altered at 24 h postinfection, and then more significantly altered at Ն48 h postinfection. Six of the 15 phosphopeptides were derived from five proteins that have known roles in the host immune system. These include both Hemgn phosphopeptides, one of the two Pram1 phosphopeptides, and the interferon-␥-inducible p47 GTPase (Gm12250) phosphopeptide mentioned above. The other two phosphopeptides originated from Dok1 and Samsn1. Because these five phosphoproteins have known immune system functions, were mildly altered in the 24 h anthrax samples, and were significantly altered in the Ն48 h anthrax samples, these constitute the highest-confidence candidate phosphoprotein biomarkers of early-stage anthrax discovered by this study. FIG. 5. Mouse spleen signaling responses to early-stage anthrax. To tentatively discover early-stage signaling responses, the analysis described in Fig. 4 was repeated except that the Ն48 h ⌬Sterne and Sterne data were excluded. Twenty-six phosphopeptides were tentatively found to be altered (qϽ0.5). Note that by requiring qϽ0.5, 50% of the data (13 phosphopeptides) were estimated to have been wrongly considered significant by the ANOVAs. Of these 26 phosphopeptides, ten were derived from eight proteins that have known roles in the immune system. Two of these eight proteins, Pram1 and Hemgn, each had two phosphopeptides that responded similarly (i.e. clustered together) at 24 h postchallenge, indicating that both proteins were probably genuinely affected during early-stage anthrax.

DISCUSSION
Recently developed sample preparation, mass spectrometry, and data analysis technologies have greatly advanced quantitative phosphoproteomics. In this study, the spleens of B. anthracis spore-challenged mice were analyzed using label-free quantitative phosphoproteomics, thousands of phosphopeptides were confidently identified and quantitated, and many candidate biomarkers of anthrax were discovered. Sterne infections resulted in global alterations to the spleen phosphoproteome starting at 48 h postchallenge, correlating well with the morphological changes to the spleens and the severity of the symptoms at this stage of the disease. Asymptomatic ⌬Sterne infections resulted in 188 significantly altered phosphopeptides, and the corresponding phosphoproteins likely have important host response functions in the spleen.
Twenty-six phosphopeptides were tentatively discovered to be altered at 24 h postchallenge. Eight of the corresponding proteins have known roles in the immune system: (i) Dido1 (death-inducer obliterator), an interleukin-3 activated apoptotic factor, (ii) Dok1 (docking protein), a negative regulator downstream of Toll-like receptors, (iii) Gm12250 (interferon-␥-inducible p47 GTPase), a regulator of immunity to intracellular pathogens, (iv) Hemgn (hemogen), a hematopoietic cellspecific protein that regulates proliferation and differentiation and resists cell death through activation of NFB, (v) Map3k5 (mitogen-activated protein kinase kinase kinase), which is activated by the T-cell costimulatory receptor and regulates apoptosis upon infection and B-cell activation and proliferation, (vi) Mllt4 (afadin), a regulator of hematopoietic cell-cell adhesion, (vii) Pram1 (PML-RARA-regulated adapter), a regulator of T-cell receptor mediated activation of, for example, NFB, and (viii) Samsn1 (SH3 SAM containing hematopoietic adaptor), a B-cell receptor immunoinhibitor. These represent excellent candidate phosphoprotein biomarkers for validation studies of circulating immune cells using, for example, triple quadrupole mass spectrometry or reverse phase protein mi-croarrays. A validated mouse spleen signaling response of early-stage anthrax might translate into a presymptomatic biomarker of human anthrax detectable within circulating immune cells. Alternatively, additional discovery experiments would likely result in a set of higher-confidence candidate biomarkers. For example, label-free quantitative phosphoproteomics could be used to study circulating immune cells isolated from subhuman primates given an intratracheal administration of a B. anthracis spore aerosol of a highly pathogenic BSL-3 strain such as the Ames strain. Additionally, phosphoproteome profiling of human circulating immune cells incubated with B. anthracis bacteria or from recipients of anthrax vaccination might be informative. Used as diagnostic markers, phosphoprotein biomarkers of presymptomatic anthrax could potentially provide physicians with a key tool to fight anthrax in the event of a bioterror attack. More broadly, they may be useful biomarkers of immune response in general.
The recent development of phosphopeptide enrichment strategies, high-resolution nano-LC-MS(/MS), software designed to analyze mass spectrometry data, and countless other technologies have enabled confident, high-throughput, shotgun phosphopeptide identification and simultaneous comparative quantitation. Besides innumerable applications to life science research, future clinical applications of quantitative phosphoproteomics might include rapid phosphoproteome profiling of tumors to differentially diagnose cancer and of circulating immune cells to test for infection, inflammation, and immune disorders. Combined with other omics technologies, quantitative phosphoproteomics might help foster the emerging science of systems biology, and together these synergistic analytical methods will hopefully lead to more comprehensive models of biological processes and the development of personalized medicine.