Cluster analysis of host cytokine responses to biodefense pathogens in a whole blood ex vivo exposure model (WEEM).

BACKGROUND: Rapid detection and therapeutic intervention for infectious and emerging diseases is a major scientific goal in biodefense and public health. Toward this end, cytokine profiles in human blood were investigated using a human whole blood ex vivo exposure model, called WEEM.
RESULTS: Samples of whole blood from healthy volunteers were incubated with seven pathogens including Yersinia pseudotuberculosis, Yersinia enterocolitica, Bacillus anthracis, and multiple strains of Yersinia pestis, and multiplexed protein expression profiling was conducted on supernatants of these cultures with an antibody array to detect 30 cytokines simultaneously. Levels of 8 cytokines, IL-1α, IL-1β, IL-6, IL-8, IL-10, IP-10, MCP-1 and TNFα, were significantly up-regulated in plasma after bacterial exposures of 4 hours. Statistical clustering was applied to group the pathogens based on the host response protein expression profiles. The nearest phylogenetic neighbors clustered more closely than the more distant pathogens, and all seven pathogens were clearly differentiated from the unexposed control. In addition, the Y. pestis and Yersinia near neighbors were differentiated from the B. anthracis strains.
CONCLUSIONS: Cluster analysis, based on host response cytokine profiles, indicates that distinct patterns of immunomodulatory proteins are induced by the different pathogen exposures and these patterns may enable further development into biomarkers for diagnosing pathogen exposure.

Background

Yersinia pestis and Bacillus anthracis are two pathogens of significant concern to public health from a biodefense perspective [1,2]. Y. pestis, the causative agent of plague, is a Gram-negative, highly communicable coccobacillus that has been responsible for three historic pandemics with high mortality rates [3-5]. The microorganism possesses a Type III secretion mechanism common to several human, animal and plant pathogens, whereby a series of pathogen-specific structural proteins form a syringe-like structure capable of injecting virulence factors into the mammalian host cell. These virulence factors then facilitate pathogen use of host nutrients and thwart the host immune response, ultimately causing cell and host death [6,7]. Naturally occurring plague can be transmitted from infected fleas and rodents to humans, and although the pathogen can be phagocytosed, it can also resist destruction by manipulating the host defense mechanism(s), potentially through antigenic mimicry [8]. Y. pestis then multiplies rapidly leading to necrosis of lymph nodes, a condition known as bubonic plague, which can result in death if untreated [2]. In some cases the infection can spread through the blood stream resulting in systemic plague (septicemia) or to the lungs resulting in the highly contagious and deadly form of the disease known as pneumonic plague. There are currently no rapid, widely available diagnostic tests for plague, and the most common treatment is streptomycin [2,3], an antibiotic with adverse effects.

Two other species from the genus Yersinia are also human pathogens: Y. pseudotuberculosis and Y. enterocolitica[9,10]. Despite their high degree of sequence similarity to Y. pestis, these two near neighbors of Y. pestis manifest in very different symptoms, ranging from abdominal pain to septicemia in humans, usually caused by infection through contaminated food. Infections caused by Y. pseudotuberculosis or Y. enterocolitica can be effectively treated with antibiotics and in most cases are self-limiting. Notably, Y. pestis is reported to have evolved from Y. pseudotuberculosis within the past 10,000 years [11].

B. anthracis is a Gram-positive, rod-shaped spore-forming bacterial pathogen and the causative agent of anthrax [12,13]. Human, livestock, and wildlife mortalities attributable to anthrax occur in numerous regions of the world, although the majority of cases are found in less industrialized nations [14]. Three forms of the disease have been described: cutaneous, intestinal and inhalational. While cutaneous and intestinal forms may be less severe, inhalational anthrax is often fatal without prompt antibiotic treatment [13]. The primary mechanisms of virulence employed by B. anthracis are associated with two virulence plasmids designated pXO1 and pXO2 [15]. The net effect of these plasmids is virtually unhindered proliferation of B. anthracis within the host, hemorrhaging, cardio-pulmonary collapse, and death.

The regulation of production of host cytokines by both Yersinia and B. anthracis has been described previously. Pickering A. K. et. al. measured cytokine levels in human dendritic cell supernatant and in mouse peritoneal macrophages exposed to B. anthracis spores [16]. They observed significant increase in TNF-α, IL-6, IL-1β, IL-8, and IL-12 in human dendritic cell supernatants by 5 hours post-exposure. High levels of IL-6, and TNF-α were observed in the supernatant from B. anthracis infected mouse peritoneal macrophages [16]. In a mouse model, 6 cytokines, namely IL-12p70, TNF, IFN-γ, MCP-1, IL-10, and IL-6, were increased significantly in mouse lung at 48 hours of Y. pestis infection [17]. In previous work comparing exposures to different bacterial pathogens, distinct patterns of cytokine expression levels were found that could discriminate the particular host response [18], including while using pathogen-specific LPS in whole blood [19].

The hypothesis for the present study is that exposure to diverse bacterial pathogen strains would result in distinct cytokine profiles in the host, with strains from the same species exhibiting more similar profiles than strains from phylogenetically distant species. A multiplex cytokine protein chip was used, and a multivariate approach was taken that combined expression data on multiple cytokines. Multivariate clustering techniques were used to establish cytokine expression profiles after ex vivo exposure of whole blood to seven pathogens.

Methods

Bacterial strains and culture conditions

The bacterial strains used in this study include: B. anthracis Ames (virulent), B. anthracis Sterne (vaccine strain), Y. pestis KIM5 D27 (attenuated, pgm-). Y. pestis India/P (attenuated, pgm-), and Y. pestis NYC (virulent), Y. pseudotuberculosis serotype 1 PB1, and Y. enterocolitica WA serovar 0:8. Bacteria were grown on tryptose blood agar slants at 26°C for 1-2 days and subsequently collected using 2 ml of 0.033M potassium-phosphate, pH 7.0;.bacterial densities were measured at OD620 (1 OD620 = 1.2 x 109 colony forming units/ml).

Whole blood ex vivo exposure model (WEEM)

Human blood was collected from a healthy donor by venipuncture using CPT Vacutainer tubes (Becton Dickinson) containing citrate. Informed consent was obtained and our blood collection protocol was approved by the LLNL IRB committee. Separate CPT tubes were used for the unexposed control and 7 different bacterial exposures (B. anthracis Ames, B. anthracis Sterne, Y. pestis NYC, Y. pestis India/P, Y. pestis KIM5 D27, Y. pseudotuberculosis, and Y. enterocolitica). Bacteria were added to blood within 15 minutes of collection at a multiplicity of infection ratio of 5:1. This ratio was determined against white blood cells in whole blood:∼7 x 106 cells/ml. Each whole blood sample was incubated with bacteria for 4 hours at 37°C in 5% CO2 Following incubation, plasma was collected by centrifugation at 2000 x g for 10 min at 4°C. The control plasma was obtained in the same way and treated with 0.033 M potassium-phosphate as a mock exposure. These plasma samples were used for cytokine measurements.

Cytokine immunoassays with protein arrays

The measurements of cytokines were performed using Zyomyx Protein Profiling Biochips (Hayward, CA). These protein arrays allow the simultaneous quantification of 30 biologically relevant cytokines, as determined by Zyomyx, Inc: IL-1α, IL-1β, IL-2, IL-3, IL-4, IL-5, IL-6, IL-7, IL-8, IL-10, IL-12(p40), IL-12(p40/p70), IL-12(p70), IL-13, IL-15, TNFα, TNFβ, Eotaxin, MCP-1, MCP-3, TRAIL, CD95(sFas), MIG, sICAM-1, IP-10, CD23, TGF-β, GM-CSF, GCSF, IFN-γ. Each cytokine assay was optimized for the Zyomyx Protein Profiling Biochip based on many factors including the availability of antibodies and the sensitivity and specificity of antibody-cytokine interactions. Each protein array chip is designed with 6 independent microfluidic channels that allow up to 6 samples to be loaded into isolated regions of an array. Antibodies specific for 30 analytes were arrayed in each channel, and each antibody was arrayed in redundancy on 5 pillars within the channel. Accordingly, a cytokine measurement represents the average of 5 measurements. All immunoassay steps, including sample loading, washing, and detection, were performed with a fully automated biochip processing station (Zyomyx Assay 1200 workstation). Eight protein array chips were used in these experiments. Two chips were used for generating calibration curves with a calibration standard kit containing 30 analytes (Zyomyx, Inc.). Sample (40 μl) was injected into each channel of the protein array chips. Standard solutions were applied to two channels of each chip for chip-to-chip normalization. Triplicates of control and pathogen-exposed plasmas were applied randomly to four channels of 6 protein array chips. Protein arrays were scanned at 532 nm with Zyomyx Scanner 100 after immunoassays. Zyomyx Data Reduction software was used for normalization, calculation of calibration curves. Dixon’s test was used to remove outliers, and the median feature intensity was background subtracted. Concentrations of cytokines in plasma samples were determined by a four parameter logistic model.

Cluster analysis of cytokine data

Multiple hierarchical clustering methods were used to group the pathogen exposures based on the multivariate cytokine expression profiles induced in a host infection model system. First, hierarchical agglomerative clustering [20] was applied to group the control and the seven pathogen-exposed samples based on their cytokine concentration profiles. Each of the eight samples was characterized by the multivariate vector of its average log10 concentration over the cytokines. The samples were then clustered based on the following distance measures between the samples and between the clusters.

Distance between two samples was defined using two distance metrics:

· Euclidean distance

· Correlation distance: (1 – Spearman correlation coefficient between the samples)

Distance between two clusters was defined using three methods:

· Complete linkage (furthest neighbor): the largest distance between members of the clusters

· Single linkage (nearest neighbor): the smallest distance between members of the clusters

· Average linkage (group average): the average distance between members of the clusters

Given a pair of distance metrics between samples and clusters, the algorithm was initialized with the eight samples forming eight different clusters and then processed iteratively by joining the two most similar clusters. The tree was built starting from the individual samples, using an agglomerative (bottom up) approach. The resulting hierarchy of clusters was displayed as a dendrogram. These traditional clustering methods provide a quick, exploratory overview of the data. However, these methods do not estimate the optimal number of clusters in the data; rather, the clustering is performed exhaustively from the lowest possible level of the hierarchy where each sample forms its own cluster, to the highest level where all samples are grouped into one cluster.

In addition to the traditional hierarchical agglomerative clustering method, the hierarchical ordered partitioning and collapsing hybrid (HOPACH) algorithm was also applied to the cytokine measurements [21]. In contrast with the previous approaches where the tree was built starting from the individual samples as the leaf nodes, HOPACH used a hybrid divisive-agglomerative approach: it started from the root cluster containing all the samples (divisive, top down approach), then divided the root down to leaf nodes, with an extra collapsing (agglomerative) step after each iteration that combined similar clusters. Based on the correlation distance between samples, HOPACH determined the split that minimized a measure of cluster homogeneity called the median split silhouette. While computationally more expensive than the previous methods, HOPACH was expected to perform better because of its dynamic approach to update and potentially revise the clusters at every step of the iteration. Furthermore, HOPACH also estimated the optimal number of clusters from the data, and thus offered another advantage over the previous methods.

Computations were performed in the R computing environment (http://www.r-project.org/) and the HOPACH package [21].

Results

Cytokine levels were examined using an ex vivo model, termed WEEM for whole blood x vivoexposure model. Individual samples of anti-coagulated human blood were incubated with B. anthracis Ames, B. anthracis Sterne, Y. pestis KIM5 D27, Y. pestis NYC, Y. pestis India/P, Y. enterocolitica, and Y. pseudotuberculosis. After 4 hours exposure, blood cells were removed by low-speed centrifugation and concentrations of 30 cytokines in the plasma were measured with protein arrays. Concentrations of fourteen cytokines, GCSF, IFNγ, GM-CSF, IL-7, IL-12(p70), IL-12(p40/p70), IL-13, IL-2, IL-3, IL-4, IL-5, MCP-3, TGFβ, and TNFβ were below the limit of detection in this study. The following 16 cytokines were detected: Eotaxin, IL-10, IL-12(p40), IL-15, IL-1α, IL-1β, IL-6, IL-8, IP-10, MCP-1, MIG, TNFα, TRAIL, sCD23, sCD95, and sICAM-1 (Figure 1). To determine if there were significant differences among the levels of cytokines in the control and pathogen exposed plasma samples, F-tests were performed. For thirteen of these 16 cytokines, all three replicates were detected and these cytokines were subjected to F-tests. Statistical analysis indicated that 8 cytokines (IL-1α, IL-1β, IL-6, IL-8, IL-10, IP-10, MCP-1, and TNFα) had differentially elevated expression profiles following different bacterial exposures. Figure 2 shows the concentrations (pg/ml) of these cytokines in the control and bacteria exposed plasma samples. The F-tests revealed that the other five cytokines containing complete datasets, TRAIL, sCD23, sCD95, MIG, and sICAM-1, had no significant difference between bacterial exposures and the mock-exposed control. Moreover, there was a great variation in absolute concentrations between cytokines. For example, the concentrations of TNFα, sCD23, and sICAM-1 were as high as 1 x 104 -105 pg/ml, whereas IL-10 was much lower, about 16 pg/ml.

Figure 1

Scatter plots of 16 cytokine concentrations detected in human blood followingbacterial exposures. Cytokine concentrations were displayed on a logarithmic scale. The cytokines shown here were detected out of the 30 cytokines in the arrays. The 8 cytokines that were found to be statistically differentially expressed among these samples are highlighted with rectangular boxes. Each mark delineates the average of triplicate exposure samples. Each exposure sample is loaded onto a protein array chip that contains 5 independent measurements per cytokine meaning that fifteen measurements are used to obtain these data.

Figure 2

Concentrations of 8 cytokines in human whole blood afterexposure to pathogens. The control was a mock-exposed sample. Cytokine concentrations were determined using protein arrays. The bars represent the average of three replicate samples that each contain 5 replicate features per cytokine assay and the lines represent the standard deviation among the three replicates.

Marked differences in induced cytokine patterns between B. anthracis and Yersinia exposures were found. Also, the levels of induction of these cytokines differed among the different bacteria. For example, Yersinia species induced much higher cytokine response than B. anthracis for IL-1α, IL-1β, IL-6, and TNF-α (Figure 2). The two strains of B. anthracis bacteria induced different levels of IL-1β and TNF-α (Figure 2), including 2.2 times higher concentration of IL-1β and 1.6 times higher concentration of TNF-α for the Sterne strain than the Ames strain of B. anthracis. These differences were statistically significant (pairwise t-test p value = 0.0039 for IL-1β and 0.022 for TNF-α). To discriminate Y. pestis exposure from near neighbors, IL-10 levels can be used, showing cytokine concentrations following Y. enterocolitica exposure and Y. pseudotuberculosis exposure that are on average 5-fold higher and 2-fold higher, respectively, than after Y. pestis exposure (Figure 2). IL-10 differential expression was specific to the Yersinia spp. because exposure to B. anthracis strains showed comparable IL-10 levels to that in unexposed control.

The HOPACH algorithm estimated the number of clusters as five, and grouped the samples based on their host cytokine expression profiles as follows: 1) Y. pestis (KIM5 D27, India/P, and NYC), 2) Y. pseudotuberculosis, 3) Y. enterocolitica, 4) B. anthracis (Ames and Sterne), and 5) Control (Figure 3). The closer the pathogen-exposed samples are within the tree on the left, the more similar they are. Height of the branches indicates the distance between the successive nodes in the clustering. The method separated the B. anthracis and Yersinia infected blood samples. In addition, the cytokine profile of the mock-exposed control was more similar to the pattern produced by B. anthracis exposure than to the profile elicited by Yersinia.

Figure 3

Clustering result with HOPACH using the average linkage distance between clusters is shown. The eight pathogen-exposed samples are clustered according to the dendrogram on the left and cluster into five groups, 1) Y. pestis (KIM5, NYC, and India), 2) Y. pseudotuberculosis, 3) Y. enterocolitica, 4) B. anthracis (Ames and Sterne), and 5) Control. Sixteen cytokines (Eotaxin, IL-10, IL-12(p40), IL-15, IL-1α, IL-1β, IL-6, IL-8, IP-10, MCP-1, MIG, TNFα, TRAIL, sCD23, sCD95, and sICAM-1) are also reordered based on their correlations according to the dendrogram on the top. Clusters go from root at top to leaf node for each cytokine. Clusters in between are based on their agglomerative . The branch shows the similarity, the short the branch, the more similar. In addition, the eight rightmost proteins form a cluster that may involve inflammation-related cascades initiated by an innate immune response to these pathogen. Colors represent units of log10 [pg/ml], in ten equally spaced intervals increasing from white to dark red. A key showing the specific log10 values for each interval is shown in the figure.

Results of the hierarchical clustering when using the Euclidean distance between samples depended on the distance metric between clusters. The three methods for determining the distance between clusters (complete linkage, single linkage, and average linkage, see Materials and Methods) all established three major clusters: 1) Y. pestis and near neighbors, 2) B. anthracis, and c) Control. Differences between the results occurred when the Yersinia cluster was further divided. The average linkage method, consistent with Figure 3, formed a subgroup of the three Y. pestis strains, then grouped them first with Y. pseudotuberculosis followed by Y. enterocolitica. Complete and single linkage methods, however, first grouped the attenuated virulent strain of Y. pestis (India/P) with the more virulent strain (NYC), both clinical isolates from human plague cases, and then clustered them with Y. pseudotuberculosis, followed by the attenuated Y. pestis (KIM5 D27), and lastly with Y. enterocolitica. This is interesting from an evolutionary perspective because it has been proposed that Y. pestis evolved from Y. pseudotuberculosis within the last 10,000 years, and thus these two pathogens are more closely related [11].

When using hierarchical clustering with the correlation distance between the samples, the final clusters were independent of the distance metric between clusters, and agreed with the tree structure in Figure 3. The complete, single, and average linkage methods all resulted in the following major clusters: 1) Yersinia, 2) B. anthracis, and 3). Control. Within the Yersinia cluster, Y. pestis (NYC) was closest to Y. pestis (India/P), followed by Y. pestis (KIM5 D27), Y. pseudotuberculosis, and Y. enterocolitica.

Discussion

The HOPACH clustering method (Figure 3) produced five distinctly separated clusters: 1) Y. pestis (KIM5 D27, India/P, and NYC), 2) Y. pseudotuberculosis, 3) Y. enterocolitica, 4) B. anthracis (Ames and Sterne), and 5) Control. This result is consistent with the findings using the correlation distance and the Euclidean distance with average linkage. In addition, HOPACH estimated the optimal number of clusters as five. That is, the Yersinia subcluster is best if it is divided into the three clusters specified by 1) through 5) above. Y. enterocolitica forms its own cluster, and so does Y. pseudotuberculosis. Y. pestis (KIM5 D27), Y. pestis (India/P), and Y. pestis (NYC) are grouped into one cluster. Further subdivisions lead to an overall clustering with inferior quality.

In addition to clustering the cytokine expression profiles across bacterial treatments, Figure 3 also groups the cytokines themselves and clusters the proteins based on their similarities across the pathogen exposures and reorders them accordingly. Interestingly, the three pro-inflammatory cytokines IL-1β, TNFα, and IL-6 clustered closely, and so did the three chemokines MCP-1, IP-10, and IL-8. Although these 6 cytokines do not cluster as a single group, they do cluster at a branch further away from the leaf node, which includes IL-10 and sCD95, to make a larger group of 8 proteins. Several of these proteins are involved in inflammatory conditions, such as IL-1beta, TNFα, IL-6, [22] and have been shown to be upregulated in cell culture and animal model specifically exposed to biothreat agents [23]. Increased expression of IL-6 and TNFα clustered together in a study involving mouse splenic CD11b + cells following sub-lethal Y. enterocolitica infection [24]. In addition, several of the cytokines in this cluster, namely TNF-alpha, IL1-beta, IL-10, and MCP-1 are expressed higher in exposed whole blood as compared to control in this study and in whole blood exposure to LPS from several other gram negative bacterial pathogens [19]. In addition to expression differences, the absence of detected cytokine expression can also be helpful in discriminating pathogen exposure.

The multiplex detection of 30 cytokines in this study revealed the early phase cytokine expression profiles in human plasma following exposures to B. anthracis (Ames and Sterne), Y. pestis (KIM5 D27, NYC and India/P), Y. pseudotuberculosis, and Y. enterocolitica. The expression levels of 8 cytokines, IL-1α, IL-1β, IL-6, IL-8, IL-10, IP-10, MCP-1, and TNFα were significantly different from that of unexposed control (Figure 2). Although the focus of our work was to show that cytokine expression profiling can discriminate between different pathogen exposures in a human whole blood ex vivo model, these results also represent an initial attempt to characterize the full cytokine response to each individual pathogen. Our preliminary study using a single exposure protocol at a single time post-exposure will need to be supplemented with more thorough investigation in order to determine the usefulness of using cytokine levels for diagnosing pathogen exposure. However, the single time point chosen, 4 hours, is sufficient to detect proteomic changes and has been used in previous studies examining cytokine levels [25-27]. This time point represents a start towards a more complete temporal study, as has been done with gene expression patterns for two of the pathogens studied here [25,27]. In addition, studies that provide expression patterns for a single cytokine using multiple time points will also be needed to make the results of this paper clinically useful, such as has been done by, Cooper and coworkers, who examined IL-12p40 and IL-12p70 levels following different growth conditions and exposure levels for a time course of Y. pestis exposed dendritic cells [28]. The results of the current work shows a similar expression pattern trend to this previous work, in which, Y. pestis induces IL-12p40 and at a substantially higher level than IL-12p70.

Our results showed that the expression levels of 3 chemokines, IL-8, MCP-1 and IP-10, were induced by both Yersinia and B. anthracis exposures. No significant differences were found for these cytokines between Yersinia and B. anthracis exposures. IL-8, MCP-1 and IP-10 are chemokines that enable the migration of leukocytes from the blood to the site of inflammation. IL-8 is a key chemokine regulating neutrophil recruitment [29]. The essential involvement of IL-8 in acute inflammation was demonstrated by neutralizing IL-8 with its antibody. When highly specific antibody against IL-8 was administered in acute inflammatory reactions induced by several stimuli including lipopolysaccharide, neutrophil infiltration was blocked [30]. MCP-1 is known for its ability to act as potent chemoattractant and activator of monocytes/macrophages as well as NK cells but not neutrophils [31,32] . IP-10 has no chemotactic activity for neutrophils but attracts monocytes, NK, and T cells to the site of infection and regulates T cell maturation [33,34]. It was reported previously that elevated IL-8 and MCP-1 were secreted by human epithelial cells after Y. enterocolitica infection, but not IP-10 [35,36]. Human dendritic cells, infected with B. anthracis spores, secreted high level of IL-8 at 7.5 hours [16]. In our study, the fold increase of IL-8 was much greater than MCP-1 and IP-10 (Figure 2). For example, the induction of IL-8 by Ames strain of B. anthracis was 41 fold, while MCP-1 was 2 fold and IP-10 was 2.5 fold (Figure 2). This result may indicate that IL-8 is a dominant chemokine in early response (4 hours exposure in our study) and neutrophils are the major player in early inflammatory response.

Here we compared cytokines induced by B. anthracis and Yersinia exposures. Overall, Yersinia exposure induced higher levels of IL-1α, IL-1β, IL-6, IL-10 and TNFα than B. anthracis exposure, suggesting these cytokines could be used to develop an assay for discriminating Yersinia spp. from B. anthracis exposures. The vaccine strain (Sterne) of B. anthracis induced higher levels of IL-1β and TNFα than the virulent strain (Ames) (Figure 2), suggesting these cytokines can contribute to a biomarker panel to discriminate if a particular isolate of B. anthracis is virulent. There was also a difference in induction of IL-10 between Y. pestis and near neighbors (Figure 2), suggesting this cytokine is a candidate biomarker for discriminating the virulence of Yersinia species. These data regarding IL-10 expression following Yersinia spp. exposure are in agreement with published literature that shows Y. enterocolitica and Y. pestis can elicit statistically different levels of IL-10 expression [37]. Differences in IL-10 induction may be due to differences in the lcrV protein among Yersinia spp.[38]. The different cytokine profiles induced by B. anthracis and Yersinia here may be partially due to different surface antigens on the outermost part of these pathogens and the manner in which these bacteria were grown. Lipopolysaccharide (LPS), the main constituent of the outer membrane of Gram-negative bacteria, and peptidoglycan (PGN), the major cell wall component of Gram-positive bacteria, have been reported to elicit markedly different immune responses [39]. However, virulence factors, such as B. anthracis lethal toxin and Yersinia virulence antigen, LcrV, may also play important roles in differential cytokine induction. This view is supported by numerous reports that B. anthracis toxin and virulence factors of Yersinia bacteria (Yops, invasin, LcrV) modulate host cytokine responses [40-51].

While the various clustering methods resulted in slightly different final hierarchies, all were consistent in separating the unexposed control from the samples exposed to B. anthracis or to the Y. pestis and near neighbors. Agreement on this level among the various clustering procedures lends more confidence to the overall results. On a more detailed level, the methods grouped slightly differently the samples exposed to the Y. pestis and near neighbors, which indicates that these samples cannot be unequivocally separated based on the current data and additional biomarkers or a larger sample set would be needed. The most advanced HOPACH method estimated the optimal number of clusters in the data as five, corresponding to the unexposed control, and the four species: B. anthracis, Y. pseudotuberculosis, Y. enterocolitica, and Y. pestis (avirulent and virulent) (Figure 3).

Information gained from the targeted protein array data for host response complements genomic [52-56], and other proteomic studies [57-60] of host-pathogen interactions. The success of the WEEM and computational method to distinguish pathogen exposure, based on host response in this initial study, is encouraging and suggests a number of possibilities for future studies to refine the findings. Comparative analysis, such as the current work, can potentially reveal the critical pathogenic mechanism(s) and host innate immune responses during infection as was previously shown for Y. pestis and Y. pseudotuberculosis[61]. Opportunities include using statistical hypothesis tests based on analysis of variance to assess the significance of the observed differences among the host-pathogen cytokine concentration profiles, as well as performing follow-up studies to focus more on the Y. pestis and near neighbor cluster. In addition, the methods can be extended to investigate host responses to diverse pathogens in multiple host model systems to cross validate the significance of the biomarkers to distinguish pathogen exposures.

Conclusion

Results from this study suggest that cytokine arrays coupled with statistical clustering methods can distinguish exposures to pathogens, including multiple strains of Y. pestis, Y. pseudotuberculosis, Y. enterocolitica, and B. anthracis. These methods differentiate both near neighbors and distant evolutionary microbes based on host response data. The distinct cytokine profiles also provide insight into both the host response and virulence mechanisms of diverse pathogens. In summary, characterization of host responses based on cytokine profiles has translational application, potentially providing the identification of infectious diseases and leading toward the ultimate goal of presymptomatic detection via sentinel surveillance of pathogen exposure and appropriate treatment.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

BC participated in the design of the study, conducted experiments, and drafted and finalized the manuscript. KM conducted experiments. IF participated in the design of the study, performed the statistical analysis, and helped draft the manuscript. PL helped draft the manuscript. SM conceived of the study, obtained funding, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.