RNA-Seq Analysis of the Host Response to Staphylococcus aureus Skin and Soft Tissue Infection in a Mouse Model

Staphylococcus aureus is a leading cause of skin and soft tissue infections (SSTI), which are primarily self-limiting. We conducted a comprehensive analysis of the host transcriptome during a S. aureus SSTI to provide insight on the protective mechanisms that thwart these infections. We utilized a murine SSTI model in which one ear is epicutaneously challenged while the other is not. We then harvested these infected and uninfected ears, as well as ears from naïve mice, at one, four, and seven days post-challenge, and performed RNA sequencing (RNA-seq) using the Illumina platform. RNA-seq data demonstrated a robust response at the site of infection. Comparison of gene expression profiles between infected ears and the non-infected ears of challenged mice defined the local response to infection, while comparisons of expression profiles of non-infected ears from challenged mice to ears of naïve mice revealed changes in gene expression levels away from the site indicative of a systemic response. Over 1000 genes exhibited increased expression locally at all tested time points. The local response was more robust than the systemic response. Through evaluation of the RNA-seq data using the Upstream Regulator Analytic as part of the Ingenuity Pathway Analysis software package, we found that changes in the activation and inhibition of regulatory pathways happen first locally, and lag behind systemically. The activated pathways are highly similar at all three time points during SSTI, suggesting a stable global response over time. Transcript increases and pathway activation involve pro- and anti-inflammatory mediators, chemotaxis, cell signaling, keratins, and TH1/TH17 cytokines. Transcript decreases and pathway inhibition demonstrate that metabolic genes and anti-inflammatory pathways are repressed. These data provide insight on the host responses that may aid in resolution of this self-limited S. aureus infection, and may shed light on potential immune correlates of protection for staphylococcal SSTI.


Introduction
Staphylococcus aureus is a normal colonizer of the human nose that is found in almost onethird of the population [1]. S. aureus features a broad and highly redundant repertoire of virulence factors that allow it to colonize and damage the host, as well as evade host immune responses and cause a variety of diseases in all areas of the body. As a result, this bacterium is a leading cause of both nosocomial and community-acquired infections in the United States [2,3]. In recent years, S. aureus infections have become more prevalent in the community, infecting patients with no predisposing risk factors [3].
Most community-acquired S. aureus infections occur in the skin and soft tissue and include boils, impetigo, cellulitis, folliculitis, and abscesses, often caused by isolates of the USA300 subtype [4]. While these infections are largely self-limiting in nature, severe invasive illness can result [5]. Therefore, an understanding of the protective immune response in the skin is important in order to elucidate potential new ways of combatting S. aureus while it remains localized in this area.
Investigators are beginning to understand some facets of the host response to S. aureus SSTI. The skin itself acts as an immunologic barrier to infection, with its surface maintaining an unsuitable temperature and pH for bacterial growth [6]. Antimicrobial peptides such as β-defensin, as well as skin commensals such as Staphylococcus epidermidis and their secreted products (e.g., phenol-soluble modulins [7]) also work to inhibit S. aureus growth [6]. If an infection does take hold, neutrophil recruitment to the site of infection has been demonstrated as crucial for subsequent protection [8], mediated by pro-inflammatory cytokines such as IL-1β [9,10] and by IL-17, the product of a TH17 adaptive response [11].
Much still remains unknown about the host response to S. aureus SSTI, and, specifically, what correlates with resolution of these infections. While microarray analyses have been performed examining the murine gene expression profile in S. aureus-infected skin [9], limitations to microarray studies can lead to incomplete results. Next generation sequencing, such as RNA-seq, can provide greater sensitivity and eliminates issues with hybridization and nonspecific detection [12]. Therefore, in this study we utilized RNA-seq to generate transcriptional profiles of skin from the ears of mice infected with S. aureus sub-type USA300 using an epicutaneous infection model [13,14]. We compared the differential gene expression of skin from infected ears to that from uninfected ears from the same challenged mice over time in order to evaluate, for the first time, the local response to infection directly at the site of SSTI. We also compared the non-infected ears from challenged mice to naïve mice in order to elucidate the systemic response to SSTI challenge. The data generated provide further insight into the host response to S. aureus over the course of a self-limited infection, both locally and systemically, which may be useful in determining potentially novel pathways that are important for clearance of the pathogen.

Ethics Statement
For all animal studies, protocols were reviewed and approved by the Institutional Animal Care and Use Committees (IACUC) of the Center for Biologics Evaluation and Research (Silver Spring, MD; permit number 2010-04). All surgeries were performed under ketamine/xylazine anesthesia, and all efforts were made to minimize animal suffering. Animals were sacrificed at the time points indicated below using CO 2 inhalation.
(Illumina) per the manufacturer's protocol. Sequencing reads were annotated and aligned to the UCSC BALB/c mouse reference genome using TopHat [16]. For analysis of the infection data, the alignment files from TopHat were used to generate read counts for each gene and a statistical analysis of differential gene expression was performed using the DESeq package from Bioconductor [17]. We made comparisons between naïve uninfected mice and unchallenged ears from infected mice, as well as between challenged and unchallenged ears from the same infected mouse, at one, four, and seven days post-challenge. Three naïve and three infected mice were evaluated at each time point. A gene was considered differentially expressed if the false discovery rate (FDR) for differential expression was less than 0.01 and the fold change was at least two-fold (Log 2 fold change (LFC) of 1).

Selection of top differentially expressed genes
All significantly differentially expressed genes for all time points are presented in the S1 Table. For simplicity of analysis, we chose to present only the 50 genes demonstrating the greatest increase or decrease in transcript levels from each time point within the manuscript. Under most tested conditions, we note that a proportion of differentially expressed genes show an up-regulation level of "inf". These values occur when readings from the comparator ear (the non-infected ear from challenged mice when examining the local response, or the naïve ear from unchallenged mice when examining the systemic response) are below the limit of detection (i.e. no reads mapped to this gene in these samples). However, it is unrealistic biologically that the readings are truly zero. Because any level of reads in the test ear (i.e., the infected ear when examining the local response, or the uninfected ear from challenged mice when evaluating the systemic response) will then give a LFC value of inf, the magnitude of the actual difference is lost; we are unable to differentiate whether this difference is a few reads or thousands of reads. Therefore, for genes with LFC values of inf, we looked at the raw read data to determine the actual read levels in the test ear. In order to report a gene with a LFC of inf as a top differentially expressed gene (Tables 1-4), we chose to include only those genes that had a read level in the test ear that was higher than the smallest number of reads seen for any gene within the top 50 list that had a LFC of a numerical value. Alternatively, we also included genes with lower read values than this minimum, but a LFC of inf at multiple time points. We hypothesize that low read numbers consistently detected in test ears with undetectable reads in the comparator ears in different cohorts of mice, at multiple time points, likely reflect a real result. We also included genes that had a LFC of inf at one time point, but a significant differential expression with a numerical LFC at another time point. Through these criteria, we believe we have narrowed the genes with a LFC of inf to most likely reflect meaningful expression changes.

Upstream Regulator Analysis
In order to identify which signaling pathways were potentially activated or repressed during S. aureus SSTI, we used the Upstream Regulator Analytic (URA) of IPA (Ingenuity Systems, www.ingenuity.com). This software evaluates the overlap among experimentally derived gene lists and the Ingenuity Knowledge Base, a comprehensive database of genetic data containing approximately 5 million findings manually selected from the scientific literature or third party databases [18]. URA calculates the statistical significance of the overlap between the experimental dataset and the database's regulatory pathways and reports a p value using Fisher's Exact Test. URA also analyzes the direction of the differential gene expression to predict the activation or repression of these pathways (for example, if Regulator X is reported in the literature to cause up-regulation of genes A, B, C, and D, and these four genes are found experimentally to have significantly increased transcript levels, URA will predict the pathway regulated by X is activated; if genes A, B, C, and D have significantly decreased transcript levels, URA will predict the pathway is inhibited) and reports an Activation Z-score. Activation Z-scores 2 are considered "activated" and scores 2 are considered "inhibited".

Sequencing Data Access
All raw sequencing reads have been submitted to the NCBI Sequence Read Archive (http:// www.ncbi.nlm.nih.gov/sra) under the accession number SRP040121. The processed gene

RNA-Seq of murine skin during S. aureus SSTI
We compared cDNA profiles of infected and uninfected murine ears at one, four, and seven days post-challenge. Sequencing reads were aligned to the mouse genome and were analyzed as described in the Materials and Methods. We defined differentially expressed genes as having a minimum of two-fold change in expression (log 2 fold change 1, P < 0.01) when comparing infected ears to uninfected ears from challenged mice to evaluate the local response at the site of infection, or when comparing uninfected ears from challenged mice to naïve mice to examine the systemic response.  We found 1477, 1298, and 1351 genes with significantly increased transcript levels in infected ears compared to non-infected ears from the same challenged mice at days one, four, and seven, respectively. There were 585 genes that exhibited significantly decreased transcript levels at the site of infection (comparing infected ears compared to uninfected ears from challenged mice) at Day 1, 195 at Day 4, and 143 at Day 7. The complete list is available in S1 Table. In order to evaluate the systemic response, we also compared the RNA-seq profiles of the uninfected ears from challenged mice to those of naïve mice. In this analysis we found 123 genes with significantly increased transcript levels at Day 1, 66 at Day 4, and 38 at Day 7. We found 21 genes with significantly decreased transcript levels in uninfected ears from challenged mice compared to naïve mice at Day 1, 9 at Day 4, and 56 at Day 7 (S1 Table). These results demonstrate that the challenged mice exhibit a higher level of differential gene expression locally at the infection site than at more distal sites.
Because of the large amount of data generated in these studies, the top 50 differentially expressed genes at each tested time point are presented in Tables 1-4. These tables were compiled by sorting the gene expression data for each comparison (infected ears to non-infected ears from challenged mice for the local response, and non-infected ears from challenged mice to naïve mice for the systemic response) at each time point (Day 1, 4, and 7), and listing the 50 genes for each that had the highest LFC for genes with increased transcript levels, or the lowest LFC for genes that had decreased transcript levels. For each of the top 50 genes, we then determined the LFC at the other time points and added these values to the tables to show the gene's change in expression over time. Many genes were highly altered in transcript levels at more than one time point; thus, none of Tables 1-4 contain 150 genes (50 for each of the time points). Table 1 contains the top genes with the highest increases in transcript levels for each time point during the local response, where we compared infected ears at each time point to the corresponding uninfected ears from the same mice. Of these top genes, 39 encode proteins that have known functions in the immune response, including 17 cytokines, chemokines, and chemokine-like proteins, as well as Cxcr1, which encodes an IL-8 receptor. Several genes encoding cytolytic granule proteins (Olfm4, Mcpt1, and Gzmk) were also top up-regulated immune genes at least at one tested time point during SSTI. Overall, the majority of genes with a known immune function that have the greatest increases at the site of infection appear to have a role in the innate response. Table 1 also shows that 22 genes encoding proteins involved in metabolism, cell proliferation, and regulation are among the genes with the greatest transcript increases at the site of infection. We also found six genes encoding epidermal proteins, and 39 uncharacterized genes, in the top 50 list of genes with the highest increases in transcript levels. We also note that the majority of the genes in Table 1 had significantly increased transcript levels at more than one tested time point, suggesting that the response is consistent over time. This is reflected in the overall list of genes with locally significant increases in transcript levels (S1 Table), with 60% of these genes demonstrating significant increases in transcript levels at multiple times (Fig 1A). Table 2 contains the genes that demonstrated the greatest increases in transcript levels systemically at each tested time point; these genes were differentially expressed when comparing the uninfected ears from challenged mice to those of naïve mice. Of these genes, 20 encode proteins with an immune function. The immune genes that are common between the local and systemic responses are Cxcl2, Cxcl3, Irg1, S100a9, Saa3, Il-1β, Trem1, and Clec4e. For each of these genes, the LFC in the local response is higher than it is in the systemic response (i.e. 8.97 vs 2.29 for the LFC of Cxcl3 at Day 1 locally and systemically, respectively). Also, while transcript levels for each of these common genes are consistently increased at each time point in the local response, only Cxcl2, Cxcl3, and Irg1 are significantly increased at all three time points when looking at the systemic response. The higher number of top genes that encode proteins that function in the immune response, as well as their greater and longer transcript level increases, suggest that the immune response is focused at the site of infection rather than systemically during staphylococcal SSTI. Table 2 indicates that 49 genes with the greatest increases in transcript levels during the systemic response (comparing non-infected ears from challenged mice to naïve mice) are of unknown function. Of the characterized genes, 31 encode proteins that are involved in metabolism, cell proliferation, and regulation, and 20 in epithelial cell and/or hair formation. However, unlike in the local response, where many of the top genes had consistently increased transcript levels throughout the infection, most of the systemically up-regulated genes showed significant differential expression in the challenged mice compared to naïve, control mice at only one time point. Like with the local response, this trend was reflected in the entire list of systemically up-regulated genes (S1 Table and Fig 1C). We also saw that some genes had among the highest LFC at one time point, but were significantly down-regulated at another (indicated with italics, Table 2). These observations suggest that genetic changes happening systemically are more transient and that change occurs more quickly than at the site of infection, where genes with increased expression levels generally appear to exhibit more stable expression over time. Table 3 contains the 50 genes with the greatest decreases in transcript levels at the site of infection, comparing infected ears to the uninfected ears from the same challenged mice. Again, many of these genes encode proteins with unknown function; those that have a known function appear important for cellular proliferation, differentiation, and metabolism. Only five of the genes in Table 3 encode proteins with an immune function. Fcer2a, which encodes an Fc receptor for IgE, has consistently decreased transcript levels at Days 1, 4, and 7, and Serpinb1c (which encodes a neutrophil proteinase inhibitor), is down-regulated starting at Day 4 postchallenge. Two genes that encode proteins important for binding to bacterial products (Bpifb2 and Lyg2, which encode proteins that bind LPS and peptidoglycan, respectively) have significantly decreased transcript levels at Day 7. Forty-six of the genes that have the greatest decreases in transcript levels at the site of infection encode proteins that are involved in metabolism/cell proliferation/regulation, with the majority of these genes demonstrating transcript decreases early during infection (Day 1 or Day 4 post-challenge). We also found that, by day 7, 16 of the genes with the greatest decreases in transcript levels in response to local infection encode proteins that are involved in epidermal formation. However, the majority of genes with the greatest decreases in transcript level at each time point encode proteins of unknown function, with 62 of the top genes falling into this category. Unlike the genes that had significantly increased transcript levels, where a high proportion of genes had transcript increases at more than one time point, genes with significant decreases in transcript levels were mostly not shared over multiple times, which was observed both with the genes with greatest decrease in transcription (shown in Table 3) and in the entire list of genes with significant transcript decreases at the site of infection (S1 Table and Fig 1B). Out of 769 genes that had significant decreases in transcript levels in infected ears compared to non-infected ears from challenged mice, only 126 were shared over multiple time points, suggesting that the decrease in transcript levels of locally-expressed genes in our SSTI model is also more transient than transcript increases. Table 4 contains all of the genes that had significant decreases in transcript levels systemically as there were fewer than 50 for each time point. Twelve of these genes encode proteins with an immune function. Ten of these twelve are decreased only at Day 7 post-challenge, and of these, six are involved in the adaptive response (Ighg2c and Igkc, which encode the gamma and kappa chains of immunoglobulin; Pax5, which encodes a protein involved in B cell differentiation [19]; Ms4a1, which encodes CD20, a regulator of B cell activation and proliferation [20], Cd22, which encodes a protein that functions in B cell signaling [21], and Cd99, which encodes a protein that plays a role in T cell migration into inflamed skin [22]). The majority of genes whose transcripts were decreased in non-infected ears from challenged mice compared to ears from naïve mice encode proteins that have functions in cellular metabolism, proliferation, or regulation, or lack a known function. Several genes that exhibited significant transcript decreases at one time point showed significant transcript increases at another (highlighted in italics in Table 4). Because the RNA preps from challenged mice used for this analysis are the same ones used in other analyses (the non-infected ear RNA used here as the test RNA is the same RNA used as the control RNA for examining at the local response) where we see very consistent expression of genes over time, we do not believe these temporal differences in gene expression are aberrant. The LFC values for systemic gene expression would also be affected by the number of reads in ears from naïve mice in these comparisons. However, because we used matched control mice that were kept in the same conditions and sacrificed at the same time as the challenged mice, we feel that any changes in gene expression due to environmental factors will be similar among all mice within a time point cohort. Therefore, we hypothesize the genes transcribed during the systemic response may have a transient expression pattern with rapid increases and decreases in transcript levels as the infection progresses.
In order to evaluate whether protein production mirrors gene expression, we performed Western blotting on homogenized ear lysates using antibodies specific for three proteins encoded by genes that were increased locally during infection. Western blots indicated that protein production for IL-1β, S100A8, and S100A9 is significantly increased upon infection (Fig 2,  lanes denoted by "I"), which confirms the RNA-seq data for the genes encoding these proteins and further validates the RNA-seq data as a whole.

Upstream Regulator Analysis of murine gene expression during SSTI
In order to funnel down the transcriptomic data, we next used the Upstream Regulator Analytic (URA) tool from the Ingenuity Pathway Analysis software package (www.ingenuity.com) to identify potential pathways that may be activated or repressed during S. aureus SSTI in the epicutaneous challenge model. This analysis software compares a curated database of genes and pathways to experimentally obtained transcriptomic data. It then uses statistical analysis to determine what percent of the genes within a pathway are differentially expressed in the data, and hypothesizes whether that pathway is activated or inhibited. If a high proportion of genes within a pathway are expressed in a manner that the database indicates signifies the pathway is functioning, URA will specify that the pathway is activated. Conversely, if a significant proportion of genes within a pathway are expressed in a manner opposite to what would be expected if the pathway is active, URA will indicate that the pathway is likely inhibited. Therefore, while URA is hypothetical, the analysis provides evidence-based suggestions of pathways that may be activated or inhibited based on transcriptomic data. It will also highlight activated or inhibited pathways that have upstream regulators that are not differentially expressed. These two factors allow for a global view of the transcriptomic data and for investigation into genes as pathway regulators that could otherwise be missed.
On Day 1 post-challenge, URA of the RNA-seq data for the local response (comparing infected ears to non-infected ears from the same challenged mice) indicated that 153 regulatory pathways were potentially activated at the site of infection. However, at this same time point, URA of the non-infected ears from challenged mice compared to naïve mice (systemic response) showed no activated pathways. At Day 4 post-challenge, there were 133 potentially activated pathways at the site of infection (local response), whereas URA of the systemic response showed only two activated pathways. At Day 7 post-challenge, URA of the RNA-seq data for the infected vs uninfected ears from the same mice (local response) showed indicated 127 possibly activated pathways, and URA of the RNA-seq data from the non-infected ears from challenged mice vs naïve mice (systemic response) showed 10 potentially activated pathways.
When evaluating inhibition of regulatory pathways, URA indicated 29 potentially inhibited pathways at the site of infection Day 1 post-challenge. On Day 4, 24 pathways were inhibited at the site of infection. When comparing uninfected ears from challenged mice and naïve mice to evaluate the systemic response, no pathways were inhibited on either Day 1 or Day 4. On Day 7 post-challenge, 23 pathways were potentially inhibited when analyzing the RNA-seq data from infected vs uninfected ears from challenged mice (local response). Systemically, three pathways were possibly inhibited when examining RNA-seq data from uninfected ears from challenged mice compared to naïve mice.

Evaluation of upstream regulator activation during the course of SSTI
We next examined the URA analyses to categorize the types of pathways that were indicated as potentially activated during SSTI, both at the site of infection and at distal skin sites. For simplification, we are presenting only the top activated and inhibited pathways here (Tables 5  and 6; full URA data are available in S2-S7 Tables). In order to choose these top pathways, we sorted the URA data based on overlap p value, with the most significant p values listed first (signifying pathways that have the most significant overlap with pathways listed in the IPA Knowledge Base). We then further sorted for top pathways by choosing the most significant pathways (based on p value) that had an Activation Z score 2 (for activated pathways) or -2 (for inhibition). When more than 10 pathways met these criteria, we chose the 10 with the most significant p values.
The top activated pathways are listed in Table 5. The majority of top pathways activated in the local response (comparing the infected ears to non-infected ears from challenged mice) were shared across all three tested time points and several were activated systemically by Day 7 post-challenge. The TNF pathway is activated locally on Days 1, 4, and 7, and systemically (comparing non-infected ears from challenged mice to naïve mice) on Days 4 and 7. The NFκB and IL-1β pathways are activated locally at all three time-points, and systemically on Day 7. The IL-1α pathway was among the top activated pathways locally on Day 1, both locally and systemically on Day 4, and systemically on Day 7. It was also activated locally on Day 7, though it was not a top pathway under this condition (Z score 5.045, P value 1.04 x10 -19 ; S6 Table). The RELA pathway was a top activated pathway locally at Days 1 and 4, and systemically on Day 7. This pathway was activated locally at Day 7 also, but again was not among the top pathways (Z score 3.822, P value 3.27 x 10 -23 ; S6 Table). Other top pathways were activated only locally at the site of infection, including IFNγ, TCR, STAT3 (not on the top pathway list at Day 4; Z score 2.459, P value 1.72 x 10 -19 ; S4 Table), TREM1 (not on the top pathway list at day 7; Z score 2.242, P value 8.72 x 10 -21 ; S6 Table), JUN (not on the top pathway list at day 7; Z score 2.910, P value 9.77 x 10 -16 ; S6 Table), and IL-27 (not on the day 1 top activated pathway list; Z score 2.168, P value 3.6 x 10 -16 ; S2 Table). The top activated pathways at all three tested time points are all involved in the immune response. In particular, these pathways suggest the importance of the TH1 (IFNγ, TREM1, IL-27), TH17 (STAT3, IL-1β, IL-17A), and overall pro-inflammatory responses (NFκB, IL-1α, RELA, JUN, ERK1/2, SELPLG). Pathway inhibition was similar across all three time points as well, with 24 out of 43 inhibited pathways shared at Days 1, 4, and 7 (S2-S7 Tables). As with the activated pathways, we also listed the top inhibited pathways in Table 6. Of these top pathways, two (JAG2 and Mir-155-5p) were inhibited locally at all three tested time points, and systemically at Day 7, and five were among the top locally inhibited pathways at all three time points but unaffected systemically (CD3, IL-1RN, TAB1, CD28, and MAPK1). Two pathways were on the top inhibited pathways lists at Days 4 and 7, and were also inhibited on Day 1, but did not make the top list for this time point (Mir-146a-5p, Z score -3.075, P value 9.5 x 10 -7 ; SOCS1, Z score -2.959, P value 5.41 x 10 -6 ; S2 Table). IL-37, which was a top locally inhibited pathway on Day 7, was also inhibited on Days 1 (Z score -2.408, P value 8.77 x 10 -7 ; S2 Table) and 4 (Z score -2.408, P value 1.22 x 10 -7 ; S4 Table) but was not a top inhibited gene at these times. The SOCS3 pathway was a top locally inhibited pathway at Day 7, and was also locally inhibited at Day 4 (Z score -2.619, P value 2.21 x 10 -7 ; S4 Table) but was not inhibited at day 1. These data further support the similarity of the response over time during S. aureus SSTI. As with the top activated pathways, the top inhibited pathways all function in the immune response. These pathways generally appear to be involved in the anti-inflammatory response (e.g., SOCS1 and SOCS3, which work to suppress cytokine production, and IL-37, which inhibits innate immunity [23]). Inhibition of these anti-inflammatory pathways may augment a pro-inflammatory response. Top upstream regulators selected based on p value of overlap and activation Z score. Data are ordered by p value of overlap. If more than 10 pathways were indicated as activated for a condition, the top 10 regulators were chosen as the 10 most significant regulators (based on p value) that had a Z score 2.
2 Z score infers the activation states of predicted regulators based on expression of the downstream genes within the pathway; a score 2 indicates activation. 3 P value of overlap evaluates whether there is a statistically significant overlap between differentially expressed genes and the genes that are regulated by the upstream regulator.
doi:10.1371/journal.pone.0124877.t005 Top upstream regulators selected based on p value of overlap and activation Z score. Data were ordered by p value of overlap. If more than 10 pathways were indicated as inhibited for a condition, the top 10 regulators were chosen as the 10 most significant regulators (based on p value) that had a Z score -2.
2 Z score infers the activation states of predicted regulators based on expression of the downstream genes within the pathway; A score -2 indicates inhibition. 3 P value of overlap evaluates whether there is a statistically significant overlap between differentially expressed genes and the genes that are regulated by the upstream regulator. doi:10.1371/journal.pone.0124877.t006 We also found that the CD3 and CD28 pathways, both critical to T cell receptor signaling [24], were consistently inhibited at all three time points. Though we did see a higher level of variability in the RNA-seq data in terms of differential expression of individual genes at the tested time points (refer to Tables 1-4, which demonstrate that a number of the genes that are highly differentially expressed at one time are not differentially expressed at other times), the URA suggests that global pathway activation and repression is similar over the course of infection.

Discussion
In order to better understand the effects of staphylococcal infection on the host, we performed a comprehensive genetic analysis of the mouse transcriptome using high-throughput sequencing of cDNA (RNA-Seq) prepared from RNA enriched from mouse ears using an epicutaneous infection model [13]. In this study, we chose to compare infected ears to uninfected ears from the same challenged mice in order to evaluate the local response at the site of infection. We also compared cDNA from the uninfected ears from challenged mice to cDNA from ears of naïve mice in order to determine differential gene expression systemically over the time course of infection. Through these analyses, we determined that a large number of genes are differentially expressed at the site of infection both early and late during SSTI, while a smaller number of genes are affected systemically. When we categorized these differentially expressed genes into pathways using the IPA Upstream Regulator Analytic (URA), we found that all of the pathways that showed the highest level of significance both locally and systemically were involved in the immune response.
URA bundles transcriptomic data together into potentially activated or inhibited pathways, which can give hints as to what is happening within the cellular environment to lead to the observed expression data results. URA can predict activation or inhibition of pathways whether or not the gene encoding the upstream regulator itself is differentially expressed, providing added benefit over examining gene expression data alone. The majority of the predicted activated pathways, in fact, did not have a differentially expressed upstream regulator (S2-S7 Tables). For example, both the IL-12 and IL-18 pathways were predicted to be activated at every time point (S2-S7 Tables), and while Il-12a and Il-12b show significantly increased transcript levels, the Il-18 gene itself was not differentially expressed (S1 Table). If one were only examining the gene expression data, the potential importance of the IL-18 pathway could be missed. These analyses broaden avenues of future study of staphylococcal SSTI using a rational, statistically-based categorization of gene expression data. However, it is important to note that the URA is theoretical, and follow-up experiments are necessary to confirm the activation or inhibition of these pathways.
The immune-associated genes that demonstrated the highest differential expression at the site of infection (Table 1) suggest a significant local inflammatory response. Our data largely agree with the microarray analysis performed by Cho and colleagues [9,11]; all 20 top up-regulated genes at four hours post-infection in the previous dataset had significantly increased transcript levels in our SSTI model at Day 1 post-infection. The activation of the IL-1β pathway at all tested time points in our data, demonstrated through URA, also supports Cho et al.'s contention that this cytokine plays an important role in SSTI. Our URA predicts that the IL-1α pathway was also activated at the site of infection at all tested time points; this pathway was activated systemically both at Day 4 and Day 7. Several pathways that are important for the TH17 response were activated at the tested time points, including the TNF, IL-6, TGF-β1, IL-21, and GM-CSF pathways. Cho and colleagues also demonstrated the importance of IL-17 in resolution of S. aureus skin infection [11]. Il-17a was also one of the top locally up-regulated genes Day 1 post-challenge in our model, and was IL-17 was activated globally as a pathway regulator (meaning that a significant number of downstream genes were differentially expressed) at all three time points.
We feel that the high degree of similarity between our Day 1 data and Cho et al.'s four hour post-infection data [9] validate the results we obtained through RNA-seq. Future studies will follow up on the biological significance of several genes that we have found to be differentially expressed in our model.
Our data expand on Cho and colleagues' analyses through testing more time points and through use of a more sensitive method of studying differential gene expression. We also evaluated both local and systemic responses to S. aureus SSTI instead of focusing only on the response at the site of infection. In order to best elucidate protective immune responses, we believe it is important to fully understand gene expression both at the site of infection and at distal sites.
While a significant proportion of differentially expressed genes were shared between Days 1, 4, and 7 (Fig 1), we also found many genes that are differentially expressed early and not late, or vice versa. For example, Saa2, which encodes an acute phase protein, has significantly increased transcript levels at the site of infection on Days 1 and 4, but is not differentially expressed at Day 7. SAA2 is important in the initial response to infection and its expression is activated by pro-inflammatory cytokines [25]. Csf2, which encodes a cytokine involved in granulocyte production, and the gene encoding Ficolin B, a pattern recognition receptor [26], both demonstrate significant transcript increases only at Day 1 at the site of infection. Conversely, the gene encoding granyzme K, a cytolytic granule found in CTLs and NK cells [27], had signfiicant increases in transcript levels at the site of infection only at Day 7. These genes' expression profiles give a sense of how the host response is changing over time.
We also found a large number of uncharacterized genes that were differentially expressed early during infection but not late, or vice versa. These genes could be potentially interesting as unknown mediators of any number of pathways, host responses, cellular metabolism, etc. Because SSTI is self-limiting [13], we hypothesize that that the genes that are differentially expressed later but not early are potentially most important for direct resolution of the infection, and may provide insight into the response that leads to clearance by 14 days post-challenge. However, those genes that are expressed only early in the infection might be important for helping to stimulate later responses, or for holding bacterial levels in check, allowing for future resolution.
Through examining gene expression over time, we also were able to demonstrate that activation of systemic immune responses lags behind that of the local response. Of the immune genes with the most significantly systemically increased transcript levels only on Day 7 (see Table 2), all except for two (Il-20 and Defb6) had significant increases in transcripts locally by one day after infection (Table 1 and S1 Table). The kinetics of infection, along with the patterns of gene expression, may help aid in dissection of temporal changes that lead to clearance.
Our gene expression data suggest that a TH1 response is occurring during staphylococcal SSTI. Ifnγ had significant increases in transcript levels at the site of infection on Days 4 and 7 (LFC = 4.15 and 4.43, respectively; S1 Table), and its pathway was also predicted through URA as activated locally at Days 1, 4, and 7 after challenge because a large number of genes affected by IFNγ demonstrated expression patterns indicative of pathway activation (S2, S4 and S6 Tables). This cytokine is important for a variety of immune functions, including activating innate immune cells and skewing the adaptive response toward TH1 [28]. IL-12 and IL-18 synergistically amplify IFNγ production [29], increasing the TH1 response. Both genes encoding IL-12 had significantly increased transcript levels at the site of infection in our model (Il-12a Day 1 LFC 2.84; Day 4 LFC 2.00; Il-12b Day 4 LFC 2.23, Day 7 LFC 4.29; S1 Table). We also saw significant local transcript increases of genes encoding the IFNγ-inducible TH1 chemokines CXCL9 (Table 1), CXCL10 (LFC 4.21, 4.59, and 4.87, respectively), and CXCL11 (LFC 2.09, 3.75, and 5.21, respectively) at Days 1, 4, and 7 [30]. Whether this response is beneficial is debatable, as Montgomery and colleagues demonstrated that the TH1/IFNγ pathway prevented protective immunity in their model of SSTI [31], and Nippe et al. demonstrated in a different subcutaneous infection model that mice skewed toward a TH1 response exhibited increased swelling at the infection site and higher bacterial loads by one week post-challenge [32]. The pathways controlled by IL-4, the hallmark cytokine of a TH2 response [28], as well as TH2 cytokine IL-13, were not activated in our model. In fact, the IL-13 pathway is predicted to be inhibited systemically by Day 7 post-challenge (S7 Table). These results are suggestive that, at least in our epicutaneous model of SSTI, the TH1 response may be augmented over the TH2 response. Previous work demonstrated that anti-alpha toxin antibody levels can correlate at least partially to protection in this model [33], but it remains unclear what mechanism of protection is most important in natural infection. The TH1 response augments the production of complement-fixing and opsonizing antibodies [34], which could be important in the response to S. aureus. Further studies will be required to dissect TH1/TH2 responses in the epicutaneous SSTI model.
Besides the mediators discussed above, a substantial number of other genes encoding proteins that are important for the innate response had significantly increased transcript levels, or their pathways were predicted as activated based on the transcriptional data, in our SSTI model. These include a variety of other cytokines and chemokines, factors important for leukocyte adhesion, components of the complement cascade, proteases found in neutrophil granules, and mast cell proteases. Since the infection is self-limiting, remains localized, and clears by day 14, it is likely that the innate response plays an important role in containing S. aureus SSTI. Neutrophils have been implicated for their importance in controlling these infections [8,9,13]. Our results suggest that other components of the innate response, such as macrophages stimulated through the TH1 response, may also play significant roles, and future work will begin to dissect these mechanisms.
Keratinocytes are now recognized for their importance in regulating the immune response in the skin [35]. These cells produce cytokines and chemokines as well as express cytokine and chemokine receptors. A considerable number of keratin genes had increased transcript levels both locally and systemically during SSTI, especially early in the infection (S1 Table). Besides their importance as major cytoskeletal proteins in the epithelia [36], some keratins have been identified as important regulators of the immune response. While keratin genes have been examined for their roles in inflammatory diseases such epidermolytic icthyosis [37], their role in the response to staphylococcal pathogenesis has to this point remained undescribed. These genes may provide novel means of targeting these infections. The gene encoding KRT17, which promotes epithelial cell proliferation as well as a TH1/TH17 response [38], had significantly increased transcript levels at the site of infection at Days 1 and 4 post-challenge (LFC 1.76 and 1.39, respectively), and Krt16 was highly expressed at the site of infection at all three time points, with LFCs greater than 5. KRT16 was identified for its importance in innate immune regulation during epithelial infection, where it may provide an important checkpoint in the pro-inflammatory feedback loop [39]. KRT1 is involved in negative regulation of the proinflammatory response in the epithelia [40], and Krt1 was significantly up-regulated at Days 4 and 7 at the site of infection in our model (LFC 1.36 and 1.77, respectively). Roth and colleagues established that Krt1 -/mice have markedly higher levels of the pro-inflammatory cytokine IL-33 [40]. We noted increased transcript levels of Il-33 in infected ears at Day 4 (LFC 1.38, S1 Table). However, by Day 7, Il-33 transcripts are no longer increased, while Krt1 transcripts remain increased. This suggests that KRT1 may play a role in dampening inflammation during S. aureus SSTI, possibly by modulating IL-33 levels. By Day 7, 50 keratin-associated genes had significantly decreased transcript levels at the site of infection (S1 Table). We hypothesize that keratins may be part of the immune balance in the host during S. aureus SSTI, through rapid early increases in keratin gene expression, perhaps to increase cell number/epithelial thickness, as well as to augment and control the innate immune response; as the infection continues, tissue damage worsens and the adaptive response begins to take over. At this point, the keratin genes are down-regulated. KRT8 has been postulated as a potential binding site for S. aureus adhesin ClfB [41]. Interestingly, the Krt8 gene is significantly down-regulated at the site of infection in our model at Day 4 (LFC -1.31, S1 Table). This lowered expression may help the host to decrease bacterial adherence in an effort to limit infection. Overall, the keratins are a largely uninvestigated area of the S. aureus/host interaction that deserves greater attention.
While the current literature studying responses to S. aureus SSTI have focused on inflammatory responses, our data suggest that genes that help to keep inflammation in check may also be important in the host response. A number of differentially expressed genes at the site of infection encode proteins with immune modulating functions. These include Nlrp12 (Table 1) [42], Nlrp6 (significantly increased at the site of infection on Days 1 and 4; S1 Table) [43], Irg1 (Table 1) [44], Olfm4 (Table 1) [45], Soscs3 (locally increased transcript levels on Days 1, 4, and 7; S1 Table) [46], zinc finger protein Zc3h12a (locally increased transcript levels on days 1 and 4; S1 Table) [47], and several Cd300 family members (both increased and decreased transcript levels, S1 Table) [48]. We also note that a significant proportion of inhibited pathways are involved in attenuation of inflammation (Table 6 and S2-S7 Tables). We hypothesize that the proteins encoded by these genes may help prevent dysregulated inflammation that could exacerbate skin damage; however, it is likely that a number of anti-inflammatory pathways must also be inhibited in order to obtain a level of inflammation necessary to clear S. aureus. This immune balance could be critical for generation of a successful response to the pathogen. Further research into these pathways is required to determine their actual role in the response to SSTI.
In summary, we present a comprehensive transcriptomic analysis of the time course of murine gene expression for up to one week after S. aureus challenge in an epicutaneous skin infection model. We examined differential gene expression both at the site of infection as well as systemically at three time points by comparing infected ears to uninfected ears from the same challenged mice, or comparing uninfected ears from challenged mice to naïve mice. We also used computational analysis of the transcriptomic data to generate predictions on potentially activated and inhibited pathways during both the local and systemic response to infection over time. Through these evaluations we found that the systemic response lags behind the local response in terms of pathway activation and inhibition. We determined that a majority of the differentially expressed immune response genes are critical for the innate, TH17, and TH1 responses. We also found that a large number of keratin-associated genes are differentially expressed over time, which could give insight into tissue damage, remodeling, and the potential immune function of these genes. The high level of differential gene expression, both at the site of infection and systemically, provides a great deal of potential avenues for further study into the response to this infection. RNA-seq does not provide information as to which gene(s) are responsible for resolution of infection in our model; instead, the data generated from this work allow us to develop testable hypotheses that may further our knowledge in this area. These data may help us understand how the host can contain the infection in this self-limiting model. Further study using other animal models of SSTI, including humanized mice, will help to determine if the differences in gene expression we elucidated may be applicable to human diseases.
Supporting Information S1 Table. Log fold change (Log 2 ) for all RNA-seq data. (XLS) S2 Table. URA of RNA-seq data from local response, Day 1. The data used for this analysis compared infected ears to non-infected ears from challenged mice. (XLS) S3 Table. URA of RNA-seq data from systemic response, Day 1. The data used for this analysis compared non-infected ears from challenged mice compared to naïve mice. (XLS) S4 Table. URA of RNA-seq data for local response, Day 4. The data used for this analysis compared infected ears to non-infected ears from challenged mice. (XLS) S5 Table. URA of RNA-seq data for systemic response, Day 4. The data used for this analysis compared non-infected ears from challenged mice compared to naïve mice. (XLS) S6 Table. URA of RNA-seq data for local response, Day 7. The data used for this analysis compared infected ears to non-infected ears from challenged mice. (XLS) S7 Table. URA of RNA-seq data for systemic response, Day 7. The data used for this analysis compared non-infected ears from challenged mice compared to naïve mice. (XLS)