Comparison of Transcriptional Signatures of Three Staphylococcal Superantigenic Toxins in Human Melanocytes

Staphylococcus aureus, a gram-positive bacterium, causes toxic shock through the production of superantigenic toxins (sAgs) known as Staphylococcal enterotoxins (SE), serotypes A-J (SEA, SEB, etc.), and toxic shock syndrome toxin-1 (TSST-1). The chronology of host transcriptomic events that characterizes the response to the pathogenesis of superantigenic toxicity remains uncertain. The focus of this study was to elucidate time-resolved host responses to three toxins of the superantigenic family, namely SEA, SEB, and TSST-1. Due to the evolving critical role of melanocytes in the host’s immune response against environmental harmful elements, we investigated herein the transcriptomic responses of melanocytes after treatment with 200 ng/mL of SEA, SEB, or TSST-1 for 0.5, 2, 6, 12, 24, or 48 h. Functional analysis indicated that each of these three toxins induced a specific transcriptional pattern. In particular, the time-resolved transcriptional modulations due to SEB exposure were very distinct from those induced by SEA and TSST-1. The three superantigens share some similarities in the mechanisms underlying apoptosis, innate immunity, and other biological processes. Superantigen-specific signatures were determined for the functional dynamics related to necrosis, cytokine production, and acute-phase response. These differentially regulated networks can be targeted for therapeutic intervention and marked as the distinguishing factors for the three sAgs.


Introduction
Staphylococcus aureus (S. aureus) is widely circulated in nature and carried by 25-33% of normal individuals in the anterior nares and skin [1,2]. The extreme penetrance of this bacteria and its ability to colonize skin, open wounds, and other surfaces makes it a serious threat in facilities that provide health care [3,4]. The myriad of exotoxins synthesized and secreted by S. aureus include the Streptococcal enterotoxins (SEs), such as SEA-SEE, SEG-SEI, SEK-SET, and SEY, and the toxic shock syndrome toxin (TSST-1). As SEA is the most common toxin in food poisoning, SEB is recognized for its potent toxicity as a biological Normal human epidermal melanocytes (NHEM) and the reagents required for culturing the cells were purchased from Clonetics ® (Lonza, Walkersville, MD, USA). Cells were maintained in Melanocyte Growth Medium (MGM) BulletKit ® according to the supplier's instructions (Lonza, Walkersville, MD, USA). Cell cultures were established at the recommended starting cell density of 10,000 cells per cm 2 and maintained in 150 cm 2 flasks at 37 • C and 5% CO 2 in a humidified incubator.
SEA, SEB, and TSST-1 were purchased from Toxin Technology (Toxin Technology, Sarasota, FL, USA). The toxins were diluted from the stock solution to 25 µg/mL in the MGM growth media (Lonza, Walkersville, MD, USA). On the day of the assay, the cells were treated with the appropriate amount of toxins to reach a final concentration of 200 ng/mL. The toxins were inactivated by adding TRIzol (Invitrogen, Carlsbad, CA, USA) at 0.5, 2, 6, 12, 24, and 48 h post-exposure (p.e.). As controls, untreated melanocytes were grown in parallel and harvested at the same time points. Each time point for each toxin was represented by a single culture.

RNA Isolation
Total RNA was isolated with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) using the manufacturer's procedure, followed by a cleanup procedure using the RNeasy MinElute Cleanup kit (QIAGEN, Germantown, MD, USA). The integrity of the extracted RNA was assessed using the 2100 Bioanalyzer instrument (Agilent Technologies, Santa Clara, CA, USA), and RNA integrity number (RIN) values were recorded.

Transcriptomic Assay and Analysis
The dual dye microarray hybridization was carried out using the SurePrint 4 × 44 K v2 Microarray Kit (Agilent Technologies, Inc., Santa Clara, CA, USA) following the vendor's protocol. Cy-5-labelled 200 ng of purified RNA was co-hybridized with Cy-3-labelled reference RNA (Agilent Technologies, Inc., Santa Clara, CA, USA) and bound to Agilent 4 × 44 k slides (Design ID: 026652). These arrays contain 41,000 unique probes targeting 27,958 Entrez gene RNAs. Following standard protocol, overnight hybridization at 55 • C was followed by a series of washes. The slides were scanned with an Agilent DNA microarray scanner and the features were extracted using the default setting of the Feature Extraction software (Feature Extraction software v.10.7, Agilent, CA, USA). The genes that displayed transcriptomic expressions at a fold change higher than 2 (fold change ≥ 2) were selected for further analysis.

Gene Expression Validation by Nanostring Assays
A custom NanoString panel (NanoString Technologies, Seattle, WA, USA) was designed for genes deemed functionally important for the current study. The results and discussion section justify our choice of genes listed in Table S1. Six genes-GIGYF2, INO80, USF2, WDR89, PPIA, and EIF2B1-were selected as housekeeping genes based on their stable expression levels in melanocytes [27]. We followed the standard nCounter instructions [28], a master-mix containing hybridization buffer, Reporter ProbeSet, and Capture ProbeSet (volume:volume ratio of 1:1:0.5) was prepared, of which 25 µL was added to 5 µL target RNA. The GEN2 Prep Station incubation time was set at the higher sensitivity setting (3 h) and 280 fields of view (FOV) were routinely captured. Analysis and normalization of the raw NanoString data was conducted using nSolver Analysis Software v3.0 (NanoString Technologies, Seattle, WA, USA).

Genomic Responses to the Three Toxins Are Characterized by Unique Host Expression Patterns
Principal component analysis (PCA) of transcriptomic expression data showed timeresolved clustering patterns of melanocytes exposed to three toxins for six treatment sequels ( Figure 1). PC1 and PC2 represented 21.7% and 16.1% of the total variance; thus, together, PC1 and PC2 represent nearly 38% of the total variance. Within the transcriptomic variance defined by PC1 and PC2, we found three distinct clusters for each of the toxin types. These time points emerged clustered following longitudinal trends. For example, 30 min and 2 h SEA p.e. time points clustered together, and this combination was labelled as the early treatment phase. The early treatment phase was distantly located in the PCA plot from the middle treatment phase and was defined by 6 h and 12 h SEA p.e. time points. Finally, the late treatment phase was defined by 24 h and 48 h SEA p.e. time points, which were juxtaposed in the PCA landscape and distally located from the middle phase. A hypothetical line connecting these three treatment phases showed a potential temporal trend. A very similar picture emerged from TSST-1. The genes responding to SEB treatment, however, showed a different clustering pattern, which was more apparent in the late treatment phase of SEB. A considerable Euclidian distance was observed between 24 h and 48 h SEB p.e. Therefore, unlike SEA and TSST-1, we included 24 h SEB p.e. in the middle treatment phase along with its original two members, namely 6 h and 12 h SEB p.e. This arrangement automatically labelled 48 h SEB p.e. as the sole candidate of the SEB late-treatment phase. Interestingly, the middle-to-late treatment phases (12 h, 24 h, and 48 h) of SEA p.e. clustered closely to the middle treatment phases (6 h, 12 h, 24 h) of SEB p.e. example, 30 min and 2 h SEA p.e. time points clustered together, and this combin was labelled as the early treatment phase. The early treatment phase was distantly loc in the PCA plot from the middle treatment phase and was defined by 6 h and 12 h p.e. time points. Finally, the late treatment phase was defined by 24 h and 48 h SEA time points, which were juxtaposed in the PCA landscape and distally located from middle phase. A hypothetical line connecting these three treatment phases showed tential temporal trend. A very similar picture emerged from TSST-1. The genes resp ing to SEB treatment, however, showed a different clustering pattern, which was m apparent in the late treatment phase of SEB. A considerable Euclidian distance wa served between 24 h and 48 h SEB p.e. Therefore, unlike SEA and TSST-1, we include h SEB p.e. in the middle treatment phase along with its original two members, nam h and 12 h SEB p.e. This arrangement automatically labelled 48 h SEB p.e. as the sole didate of the SEB late-treatment phase. Interestingly, the middle-to-late treatment ph (12 h, 24 h, and 48 h) of SEA p.e. clustered closely to the middle treatment phases (6 h, 24 h) of SEB p.e. Since neither of the time point experiments have technical or biological replicates present strategy of grouping time sequela into the early, middle, and late treatment ph essentially enhanced the statistical confidence of the overall results. Using the longi nal patterns of transcriptomic expressions, we sub-grouped the genes in three sets: ( 'Early' gene group, in which the transcriptomic fold changes were greater than |2| f least one of the two time points (30 m and 2 h p.e.) of the early treatment phase; (ii 'Consistent' gene group, in which the transcriptomic fold changes were greater than in all time points, and (iii) the 'Late' gene group, in which the transcriptomic fold cha were greater than |2| for at least one of the two time points (24 h and 48 h p.e.) of the treatment phase. The exception was the SEB treatment, for which the late treatment p included only 48 h p.e. Next, we combined (i) and (ii) to form 'Early-Consistent' Since neither of the time point experiments have technical or biological replicates, the present strategy of grouping time sequela into the early, middle, and late treatment phases essentially enhanced the statistical confidence of the overall results. Using the longitudinal patterns of transcriptomic expressions, we sub-grouped the genes in three sets: (i) the 'Early' gene group, in which the transcriptomic fold changes were greater than |2| for at least one of the two time points (30 m and 2 h p.e.) of the early treatment phase; (ii) the 'Consistent' gene group, in which the transcriptomic fold changes were greater than |2| in all time points, and (iii) the 'Late' gene group, in which the transcriptomic fold changes were greater than |2| for at least one of the two time points (24 h and 48 h p.e.) of the late treatment phase. The exception was the SEB treatment, for which the late treatment phase included only 48 h p.e. Next, we combined (i) and (ii) to form 'Early-Consistent' gene groups; similarly, (ii) and (iii) were combined to form 'Late-Consistent' gene groups. These gene groups were used for functional analysis.

Differences in Transcriptional Regulation in Response to the Three Toxins
In agreement with the PCA trend, the number of genes showing altered transcription varied greatly in response to the three toxins ( Figure S4). Comparisons of early and late genomic responses to each of the toxins showed differences that were at their maximum after SEB treatment in both up-and downregulated genes. The largest number of genes responding with fold change (FC > |2|, nearly 1100 genes) were observed in SEB-L, whereas nearly 100 genes showed FC > |2| in SEB-E. In contrast, SEA-E and SEA-L comprised the least number of genes with FC > |2|. Nearly 450 genes showed transcriptomic modulations at early time points and nearly 600 genes were modulated during the late time points. Treatment with TSST-1 toxin elicited a response somewhat like SEA p.e. Interestingly, there was a common trend among all three toxins: the number of perturbed genes increased with the progression of treatment time, indicating the transcriptomic storm typically augmented by this family of sAgs [29][30][31][32] ( Figure S4).

Biological Networks and Functions That Were Differentially Regulated by the Three Toxins
Functional analysis was performed using the genes listed under SEA-E, SEA-L, SEB-E, SEB-L, TSST-1-E, and TSST-1-L, respectively, to elucidate the time-dependent, toxin-specific enrichment profiles of biological and canonical functions. Table S2 lists the top biological functions (p < 0.001) and canonical networks (p < 0.01) associated with the three early treatment categories, SEA-E, SEB-E, and TSST-1-E. The list was filtered to include only those biological processes which were significantly enriched and functionally relevant to cell survival and the defense and maintenance of skin cells. In a similar fashion, genes belonging to the late treatment phase were probed to generate a list of significant biological and canonical processes that were enriched due to the prolonged toxin exposure (Table S3). Table 1 lists the top biological pathways (p < 0.001) and canonical functions (p < 0.01) that represent melanocytes' dendritic cell-like (DC like) or macrophage-like property. 'Antigen presentation pathways', 'dendritic cell maturation', 'IL17 signalling', and 'chemokine signalling', among others, emerged as the top functions that are related to melanocytes' immunogenicity.
ILK signalling emerged as a significant network that was conserved between the early and late treatment phases in response to all three sAg. Functional annotation of the 36 genes (Table S4) associated with the ILK signalling pathway demonstrated association with two cellular processes, namely the cell death and tight junction signalling. Other networks that responded in common to at least two toxins and were conserved throughout the time-course of the study include acute phase response signaling, the antigen presentation pathway, the complement system, and agranulocyte adhesion and diapedesis.
The Venn diagram in Figure S5A elucidated those biological networks that were common among as well as exclusive to SEA-E, SEB-E, and TSST-1-E. Nine networks related to cell survival and maintenance were affected by all three toxins. SEA-E and TSST-1-E shared the largest number (28) of networks, including those, which were associated with endometriosis, proliferation of connective tissue cells, and angiogenesis. SEA-E and SEB-E shared the smallest group of networks (2), which were related to skin disorders such as chronic skin disorder and chronic psoriasis.
A Venn diagram of the functional annotation enriched by the three late treatment phases, namely SEA-L, SEB-L, and TSST-1-L, (Figure S5B), demonstrated a cohesive picture of the early treatment phase ( Figure S5A). The number of overall annotated networks was greater for the late phase (87 as compared to 66 networks for the early treatment phase), as described in Tables S2 and S3. The largest number of networks was shared between SEA-L and TSST-1-L as in the early treatment phase, with similar enriched networks, namely endometriosis and proliferation of connective tissue cells. A total of 19 networks were commonly enriched for SEA-L and SEB-L; hence, the late treatment phase was associated with a higher number of significantly enriched gene networks than those associated with the early treatment phase. All three sAgs induced responses highly enriched for three biological processes: necrosis, skin diseases, and inflammation. Separate hierarchical clustering was performed using three gene sets, namely 217 genes from the necrosis network (Figure 2), 53 genes from the inflammation network ( Figure S6), and 167 genes from the skin diseases network ( Figure S7). The clustering analysis in Figure 2 identified four distinct groups of genes (indicated within yellow borders and labeled as groups A-D in Figure 2), which could be exclusive necrosis markers for TSST-1-L, TSST-1-E, SEA-L, and SEA-E, respectively. This hierarchical analysis failed to mine any exclusive signatures for SEB-E and SEB-L, respectively. Both the inflammation ( Figure S6) and the skin diseases ( Figure S7) clusters were mined as a single set, each under SEB-L (labelled group A in Figures S6 and S7, respectively). The complete list of all six gene sets is compiled in Table S5.
Biomedicines 2022, 10, x FOR PEER REVIEW 8 of 15 ( Figure S7). The clustering analysis in Figure 2 identified four distinct groups of genes (indicated within yellow borders and labeled as groups A-D in Figure 2), which could be exclusive necrosis markers for TSST-1-L, TSST-1-E, SEA-L, and SEA-E, respectively. This hierarchical analysis failed to mine any exclusive signatures for SEB-E and SEB-L, respectively. Both the inflammation ( Figure S6) and the skin diseases ( Figure S7) clusters were mined as a single set, each under SEB-L (labelled group A in Figures S6 and S7, respectively). The complete list of all six gene sets is compiled in Table S5.

Confirmation of Expression Pattern for Select Genes from the Necrosis Clusters
We performed validation of gene expression levels by NanoString nCounter ® technology. Table S1 lists the top thirteen highly perturbed genes (up-and downregulated) grouped under necrosis. This list is limited to genes responding only to SEA and TSST-1 for two reasons: first, none of the genes responding in the SEB-E phase were grouped in the three clusters discussed above (Table S5), and second, a lack of sufficient RNA samples for the SEB 48 h treatment point forced us to exclude genes that belong to the SEB-L treatment phase.
Overall, a positive correlation was observed between the NanoString and microarrays results. Of the 13 genes tested, 12 genes followed the same directionality of fold changes for the NanoString and the microarray results ( Figure S8) with the exception of one gene (PLCB1).

Discussion
The present study investigated in vitro host gene expression patterns induced by SEA, SEB, and TSST-1 during six time points ranging from 0.5 h to 48 h post-toxin exposure. A less frequently tested human skin cell type, but a major component of skin cell-mediated immunology, namely melanocytes, were selected as the target cells. The hybrid character of melanocytes was highlighted as we mined those biofunctions which were linked to the dendritic cell activities and/or the macrophage-based immune responses. This study could have benefited from incorporating additional time points to enhance the resolution of sequential biological events. For instance, our data suggested that the dosages of SEA and TSST-1 used for melanocyte treatments were potentially exhaustive within 24 h p.e.; in this context, extended time points could be highly informative. Furthermore, additional replicates in this study would result in better statistically significant gene identification. To mitigate this drawback to some extent, we mined the networks that met the cut off p < 0.05 using hypergeometric tests.

Distinct Temporal Trend of Pathogenesis Initiated by sAgs
The three toxins SEA, SEB, and TSST-1 of the sAg family are distinct in their structural, functional, and mechanistic properties [7,8,33]. Present literature not only lacks an understanding of molecular pathogenesis underlying the sAgs' toxicity, but also fails to fully comprehend the role of melanocytes in response to sAgs. The melanocytes' dendriticlike nature and their strategic location in the superficial layers of skin qualify them to be excellent mediators of initial immune defense against the sAgs [16][17][18][19][20][21][22]. We presented a whole genome-level investigation to compare the melanocytes temporal responses to SEA, SEB, and TSST-1.
A striking observation when comparing SEA and TSST-1 was the similarity in their gene expression patterns across the p.e. time course. Although SEA and TSST-1 share weak overall structural homology, TSST-1 can be displaced by SEA due to shared MHC class II binding sites [33]. This sheds light on the similarities in their mode of action as evidenced by the maximum number of shared networks for both early and late treatment phases.
Compared to SEA and TSST-1, the magnitude of transcriptional response perturbed by SEB was relatively smaller during the early treatment phase. However, the number of genes perturbed by SEB sequentially ramped up. This sort of delayed response is typical for any tissue that is not enriched with lymphocytes, as they are not the direct cellular targets of SEB [14]. Subsequently, SEB caused considerable genomic perturbations between 24 h and 48 h p.e. This trend is to be expected, as SEB typically causes a rapid neutrophil cell death accompanied by vascular congestion and leakage 24 h p.e., causing a shift to a predominantly adaptive immune response [30]. A perturbation in eNOS signalling pathways, potent vasodilators, was reported in the current study.
Another important observation was the temporal differences between SEA-and SEBinduced pathogenesis, particularly during their middle-to-late treatment phases (Figure 1). Nevertheless, a certain cohesiveness emerged between these two sAgs at the functional level. There were 11 and 19 networks that were synchronously enriched by both SEA and SEB at the early and late treatment phases ( Figure S3A,B). This fact may demonstrate an underlying similarity in their mode of pathogenesis. Early pathogenesis caused by SEAand SEB-perturbed genes manifested in skin disorders. In concurrence, SEB exclusively targeted genes linked to T-lymphocytes and their related functions, whereas SEA targeted glucose and protein metabolism networks. The consequences may include dysregulation of immune functions, apoptosis and cell death.
All three toxins enriched several networks related to cell death at early exposure phase and this response continues throughout the time course of the study. This response could be attributed to the moderately high doses of toxins used in the present study. Even though the three toxins perturbed the similar networks during the early exposure phase, as time progressed, each toxin had its unique mode of action in achieving the outcome manifested by cell death and apoptosis. One of the networks that was consistently perturbed by all three toxins across the p.e. time-course was ILK signalling. ILK functions as a kinase and signal transmitter or as a scaffold protein to facilitate cell-matrix interactions, cell signalling, and cytoskeletal organization [34]. These signals control processes related to survival, proliferation, differentiation, adhesion, migration, contractility, and neovascularization. Inhibition of ILK arrests the cell cycle and promotes apoptosis [35]. This is a key observation to support the following argument.
Early perturbation of genes associated with superoxide radical degradation in SEA indicates an oxidative stress-driven early onset of cell death [36]. TSST also perturbed this mechanism at later time points. Treatment with SEA and TSST down regulated the transcriptional levels of SOD1 and TYRP1, which potentially diminished the synthesis of different isoforms of superoxide dismutase (SODs). The potential loss of SODs highlighted the onset of oxidative stress initiated by the toxins [37], ultimately leading to onset of apoptosis during the late phase p.e.
Additional aspects of the apoptotic network, such as ERK/MAPK, were enriched by SEA and SEB at early p.e. phases, which appears to show a SE-induced apoptotic pathway distinct from that induced by TSST-1 [38,39]. We observed increased expression levels of FOS and NFAT genes during early p.e. SEA and SEB treatments. The FOS gene encodes the proto-oncogene c-FOS protein and NFATs, which are known widely for their cytokine gene expression properties and have been increasingly shown to regulate other genes related to cell cycle progression, cell differentiation, and apoptosis [39,40]. Late phase, SEB p.e. up-regulated genes that encode oncoproteins, such as Rho GTPase, which is also linked with ERK/MAPK [41]. Consequently, the G2/M DNA Damage Checkpoint Regulation, a critical biofunction closely linked with apoptosis, was highly perturbed. At late p.e. phase, SEA cross activated PI3K/AKT signaling, a critical pathway which affects many intracellular processes, including cell survival, growth, and migration.

Late Phase SEB Is Associated with Certain Dermatological Disorders
sAgs have long been implicated in the development of various inflammatory skin diseases such as psoriasis, atopic dermatitis, Kawasaki Syndrome, etc. [42,43] We observed that all three toxins modulated genes associated with the pathogenesis of psoriasis and chronic psoriasis starting from the early treatment phase. Psoriasis is often associated with functions like cell death, inflammation, autoimmune syndrome, and the production of ROS and nitric oxide [44]. From early to late treatment phases, SEA and TSST-1 shifted the expression of the gene enriching networks that are linked to lichen planus and endometriosis. During the late treatment phase, SEB regulated two unique set of genes that are closely linked to psoriasis and dermatomyositis, respectively. These genes are listed under their respective disease names in Table S3. Concurrent enrichment of oxidative stress networks could be related to the NRF2-mediated oxidative stress response and eNOS signaling pathways. Together these networks typically compromise the host's antioxidant defense mechanisms, a hallmark indicator of psoriasis [45].

Several Genes of Immunological Networks Are Differentially Modulated by Toxins
The skin exhibits a highly specialized innate immune response to invading pathogens and external stimuli. The major immune players-keratinocytes, Langerhans cells, dendritic cells, resident T-cells, and innate lymphoid cells-act in a coordinated fashion, from sensing the external stimuli to communicating through inflammatory signalling cascades, to ultimately regulating immune homeostasis [46,47]. Accumulating evidence uncovered a hybrid role of melanocytes in regulating innate and adaptive immunity [16][17][18][19][20][21][22]24,[48][49][50][51][52][53][54]. Similar to keratinocytes, melanocytes express several types of toll-like receptors (TLRs) and have the ability to produce several pro-inflammatory cytokines and chemokines [48,52,54]. Melanocytes also regulate the adaptive immunity through their functional similarities to lysosomes, such as capability to phagocytose and their antigen presentation and processing aptitudes [20,48,55]. In this context, we listed those networks (Table 1) which are associated with melanocytes' hybrid role in responding to sAgs.
All of the toxin-induced adaptive immune responses could be attributed to the networks associated with leukocyte (granulocyte/agranulocyte) adhesion, a marker for second tier responses to inflammation induced by infection. Although all toxins contributed to adaptive immunity simulation, the patterns of cytokine production and acute-phase responses differed among the three toxins. For instance, during the early treatment phases of both SEA and SEB, the cytokine and chemokine signalling networks were comprised of CXCL1, CXCL12, and PLCB1, which control leukocyte trafficking; CCL2 and CCL7, which are involved in monocyte migration and macrophage recruitment; and CFL1, which regulates cell morphology and cytoskeletal organization. Early host responses to SEB and TSST-1 included an acute phase response signal that typically triggers non-specific inflammation, leukocytosis, complement activation, protease inhibition, clotting, etc. These responses persisted until 48 h p.e.
All three toxins perturbed IL-17 signalling, a pro-inflammatory signal that bridges innate and adaptive immune responses by playing critical roles in T-cell activation and in promoting the expansion and recruitment of innate immune cells, such as neutrophils [56]. The IL-17 signalling pathway was implicated in response to toxins via alterations of the transcription of several genes in this network, including CXCL1, CXCL5, CXCL8, CCL2, CCL20, and MAP2K6.

Conclusions
To our knowledge, this is the first mRNA-level study describing the temporal response of human melanocytes to three staphylococcal superantigenic toxins, namely SEA, SEB, and TSST-1. We observed distinct temporal patterns of transcriptomic regulation for the three individual toxins. The majority of the identified networks were related to necrosis and inflammation, in agreement with previous publications [38][39][40], although most of the past studies targeted different cells than melanocytes. Pathways related to innate immunity, such as the patterns of cytokine production and acute-phase response, showed toxin-specific regulation. The time-resolved response to SEB assault took a more differential pattern than SEA and TSST-1. In conclusion, these three toxins followed distinguishable pathways to achieve a common endpoint manifested by the cell death coordinated with apoptosis and necrosis. Hence, the temporal knowledge of their pathogenesis could be the key to customized intervention. There are 2 networks commonly perturbed by SEA and SEB. Likewise, 19 and 3 networks are commonly perturbed by SEA and TSST-1, and SEB and TSST-1, respectively. There are 9 non-canonical networks that were perturbed by all three sAgs. All of these networks are listed in the diagram. (B) Common and unique non-canonical pathways enriched by the three sAgs at late phases of pathogenesis. For instance, there are 24, 30 and 3 networks are uniquely perturbed by SEA, SEB and TSST-1 at early time points. There are 10 networks commonly perturbed by SEA and SEB. Likewise, 11 and 0 networks are commonly perturbed by SEA and TSST-1, and SEB and TSST-1, respectively. There are 9 non-canonical networks that were perturbed by all three sAgs. All of these networks are listed in the diagram; Figure S6: Hierarchical clustering analysis using of 53 genes with log2 fold change > |2| enriching the inflammation pathway. Euclidian algorithm is used to sort both conditions and genes. Each block represents one gene, and its color code is at the bottom right. The conditions from left to right are named as TSST-1-L, TSST-1-E, SEA-L, SEA-E, SEB-L and SEB-E, which represent TSST-1-L at late time point, TSST-1-E at early time point, SEA-L at late time point, SEA-E at early time point, SEB-L at late time point and SEB-E at early time point, respectively; Figure S7: Hierarchical clustering analysis using of 167 genes with log2 fold change > |2| enriching the skin disease pathway. The Euclidian algorithm is used to sort both conditions and genes. Each block represents one gene, and its color code is at the bottom right. The conditions from left to right are named as TSST-1-L, TSST-1-E, SEA-L, SEA-E, SEB-L and SEB-E; Figure S8: Targeted gene expression analysis using the NanoString platform to validate the microarray data. Validation study includes a set of 13 genes from three different conditions, namely SEA at early time point, SEA at late time point and TSST-1 at late time point. The color code profile is at the bottom left; Table S1: Top 13 highly perturbed genes grouped under the necrosis cluster; Table S2: Top biological functions and diseases (p < 0.001) and canonical functions (p < 0.01) identified through IPA for early post-exposure SEA, SEB, and TSST-1 treatments; Table S3: Top biological functions and diseases (p < 0.001) and canonical functions (p < 0.01) identified through IPA for late post-exposure SEA, SEB, and TSST-1 treatments (p < 0.001); Table S4: A list of 36 genes highly perturbed by one of the three sAgs during the early and late post-exposure phases; Table S5: List of significantly different genes that enrich necrosis, inflammation, and the skin diseases pathways, respectively. Genes are sorted by their fold changes, that is, log 2 transformed. The necrosis network is perturbed by TSST-1 at the early (TSST-1-E) and late (TSST-1-L) time points, and SEA at the early (SEA-E) and late (SEA-L) time points. Likewise, networks linked to skin disease and inflammation, respectively, are perturbed by SEB at the late time point (SEB-L).