Gene Expression Profile Analyses of the Skin Response of Balb/c-Nu Mice Model Injected by Staphylococcus aureus.

BACKGROUND: Pathogenesis and persistence of many skin diseases are related to Staphylococcus aureus (S. aureus) colonization. S. aureus infection can cause varying degrees of changes in cell gene expression, resulting in complex changes in cell phenotype and finally changes in cell life activities.
MATERIALS AND METHODS: The transcriptomes of healthy and Staphylococcus aureus (S. aureus)-infected murine skin tissues were analyzed. We identified 638 differentially expressed genes (DEGs) in the infected tissues compared to the control samples, of which 324 were upregulated and 314 were downregulated, following the criteria of P < 0.01 and |log2FC| > 3. The DEGs were functionally annotated by Gene Ontology (GO), KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway and the protein-protein interaction (PPI) network analyses.
RESULTS: The upregulated DEGs were mainly enriched in GO terms, such as response to stimulus, immune system process and signal transduction, as well as in the complement and coagulation cascade pathway. Thus, S. aureus infection likely activates these pathways to limit the influx of neutrophils and prevent skin damage. Four clusters were identified in the PPI network, and the major hubs were mainly related to cell cycle and proliferation, and mostly downregulated. The expression levels of Nox4, Mmrn1, Mcm5, Msx1 and Fgf5 mRNAs were validated by qRT-PCR and found to be consistent with the RNA-Seq data, confirming a strong correlation between the two approaches.
CONCLUSION: The identified genes and pathways are potential drug targets for treating skin inflammation caused by S. aureus and should be investigated further.

Introduction

Staphylococcus aureus (S. aureus) is a conditional pathogen that forms part of normal intestinal flora of mammals but can infect other organs, especially the skin, nose and pharynx.1 Historically, S. aureus infections have been obtained within hospital environments, but the incidence or occurrence of S. aureus infections is increasing in community settings.2 Normally, S. aureus does not cause harm to human carriers; while, it can cause both invasive and noninvasive infections. Recurrent S. aureus skin infections have been correlated to atopic dermatitis (AD),3–5 and almost 80% of the individuals with AD show cutaneous S. aureus colonization.6–9 Based on its susceptibility to beta-lactams, S. aureus is commonly described as methicillin-susceptible S. aureus (MSSA) or methicillin-resistant S. aureus (MRSA).

A series of host responses will be occurred induced by both MSSA and MRSA. After S. aureus reaches the sub epidermal space, it is either locally controlled or spread in the dermis, resulting in soft tissue infection. These different processes depend not only on the virulence of S. aureus, but also on the defense ability of the host. Inflammation is the key pathological aspect of dermatological S. aureus infections, although the skin itself acts as an immune barrier, because of its unsuitable surface temperature and pH, existing antimicrobial peptides, and skin commensals (Staphylococcus epidermidis and their secreted products). However, if infection does exist, the body produces a series of pro- and anti-inflammatory mediators, chemotaxis, cell signaling, keratins, and TH1/TH17 cytokines to mediate neutrophil recruitment, and other acute and systemic immune response. Besides, S. aureus modulates the host environment through super-antigens that stimulate pro-inflammatory cytokine production by the keratinocytes, and disrupt the skin barrier. Polymorphonuclear neutrophils (PMNs) play a particularly important role in effectively eliminating S. aureus infections, and they contribute to the massive inflammation typically accompanying S. aureus infections.

The rapid evolution of whole transcriptome profiling using next-generation high-throughput RNA sequencing (RNA-Seq) has shed light on the complex mechanisms and pathways of skin infections.10 Investigators are beginning to understand some facets of the host response to S. aureus infection using high-throughput sequencing technology, especially the responses induced by MRSA strains, which are often resistant to multiple antibiotics.

Malachowa et al established a soft tissue model of infection caused by USA300, an epidemic community-associated methicillin-resistant S. aureus (CA-MRSA) strain in a rabbit to investigate host–pathogen interactions. Of the host genes tested, they found that transcripts encoding IL-8, IL1β, oncostatin M-like, CCR1, CXCR1 (IL8RA), CCL4 (MIP-1β) and CCL3 (MIP1α)-like proteins were up-regulated during S. aureus abscess formation.11 Brady et al established a murine SSTI model infected with S. aureus sub-type USA300 (MRSA) in ear epicutaneously and compared the infected and non-infected skin from challenged mice. Besides, comparisons of expression profiles of non-infected ears from challenged mice to ears of naïve mice were also discussed. In their study, the activation and inhibition pathways happen first locally, and lay behind systemically. However, they did not focus on the comparison of the infected skin to the normal skin from different individuals, which also contains a lot of information during the host response.12 In addition, transcriptomic profiling of healthy and inflamed skin biopsies has revealed significant differences in the expression levels of genes related to inflammation, type 2 immune response pathways, cornification and lipid biosynthesis.13,14

However, published information is still limited. The differences of host responses between MSSA and MRSA infections should be discussed. The understanding of the molecular mechanism induced by MSSA will be added in the existing knowledge. To this end, we established an intradermal S. aureus infection model in mice and compared the transcriptomes of the inflamed versus normal skin tissues using the RNA-Seq technology. The differentially expressed genes (DEGs) were functionally annotated by GO mRNA expression of some randomly selected genes was validated by qRT-PCR, and showed consistency with the RNA-Seq data.

Materials and Methods

Experimental Materials

S. aureus S1 strain was purchased from China General Microbiological Culture Collection Center (CGMCC Number: 1.4519) and identified as S. aureus (ATCC Number: 25923). Following overnight culture on a nutrient agar (NA) at 37°C, single colonies were picked and inoculated into a nutrient broth (NB). The inoculum was cultured overnight at 37°C with constant shaking at 200 rpm, followed by 2% subculture at 37°C for 2h to obtain mid-logarithmic phase bacteria. The broth culture was centrifuged, and the bacterial cell pellet was washed and resuspended in normal saline at the density of 3×107 CFU/100μL.15 The absorbance at 600 nm was measured to estimate the number of CFUs, which was verified by overnight culture on NA.

Establishment of Inflammation Model

Ten Balb/c-NU 6–8 weeks old male nude mice weighing ~20 g were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China). They were housed in the animal facility under standard conditions of light and temperature and fed standard laboratory chow and water ad libitum. After being fed adaptively for a week, the mice were then randomly divided into the control and model groups. Intradermal S. aureus infection was established as previously described.16 The mice were injected subcutaneously into their right hind limb with 100μL S. aureus suspension (3×107 CFU), and the control mice received the same volume of sterile normal saline. The following day, the mice were anesthetized by intraperitoneal injection of tribromoethanol (98.5%; J & K Scientific Ltd) at the dose of 240 mg/kg BW. The inflamed and normal skin were peeled, and divided into two groups. Samples in one group were fixed with paraformaldehyde and embedded in paraffin. The tissue sections were stained with hematoxylin and eosin (H & E), and the thickness of stratum spinosum and stratum epidermis were measured using the image analysis software program Case Viewer. The other skin samples in another group were taken for RNA isolation.

RNA Isolation

Total RNA was extracted from the skin samples using TRIzol® Reagent (Sigma Technologies, St Louis, USA) according to the manufacturer’s instructions, and genomic DNA was removed with DNase I (TaKara Bio, Japan). RNA quality was determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara CA, USA) and quantified using NanoDrop 2000 (Thermo Scientific™ Technologies, Wilmington, DE, USA). Only high-quality RNA samples (OD260/280=1.8~2.2, OD260/230≥2.0, RIN≥6.5, 28S:18S≥1.0, >2μg) were used to construct the sequencing library (see supplementary material for technical details).

Screening for Differentially Expressed Genes (DEGs)

The critical information like platform, quality control, bioinformatics tools and the RNA sequencing data analysis (quality control, assembly, and mapping) were exhibited in detail in the supplementary file. in the supplementary file showed the summary of trimming and read mapping results of the sequences generated from the skin samples with or without the injection of Staphylococcus aureus. A P-value and the fold change (FC) for each gene were calculated to denote its expression difference between libraries.

The DEGs between the control and inflamed skin samples were screened using P-value < 0.01 and |log2FC| > 3 as the thresholds. The DEGs were functionally annotated using Database for Annotation, Visualization and Integrated Discovery (DAVID ver 6.8, ). The significant gene ontology (GO) terms for biological process (PB), molecular function (MF), and cellular component (CC) for the up- and down-regulated genes were identified on the basis of enrichment score >2.17 In addition, the significantly enriched pathways were identified from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using the P < 0.05 as the criterion.

Construction of Protein–Protein Interaction Network

The PPI network was constructed using Cytoscape ver. 3.5.1. Weighted PPIs of human proteins were downloaded from HPRD (Human Protein Reference Database, , Release 8).18,19 The DEGs were converted from mouse to human using the HCOP: Orthology Predictions Search in HGNC database tool ().20–22

Quantitative RT-PCR

The RNA was reverse transcribed using the FastQuant RT Kit (with gDNase, TIANGEN Biotech, Beijing, China) in a 20µL reaction volume according to the manufacturer’s instructions. QRT-PCR for Nox4, Mmrn1, Mcm5, Msx1, Fgf5 and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was performed using SuperReal PreMix Plus (SYBR Green, TIANGEN Biotech, Beijing, China) and gene-specific primers () on AmpliTaq Gold (PerkinElmer, USA). All reactions were run in triplicates and GAPDH was used as the internal control. The cycling parameters were as follows: initial denaturation for 2 min at 95°C, followed by 40 cycles of 15s at 95°C, 20s at 57°C and 30s at 72°C, and a final extension at 95°C for 15s.

Data Availability Statement

All raw sequencing reads have been submitted to the NCBI Sequence Read Archive () under the BioProject accession number PRJNA647241.

Statistics

Statistical analysis was conducted using GraphPad Prism software version 7.0. All data were presented as mean value where error bars indicated standard deviation (SD). Two-tailed Student’s test was used to compare two groups, and P < 0.05 was considered statistically significant.

Results

Histological Evaluation of the Skin

Localized S. aureus infection significantly altered the epidermal morphology. Innate immune recognition of S. aureus promotes proinflammatory signals, leading to neutrophil recruitment and abscess formation, which is a hallmark of S. aureus infection. As shown in Figure 1A and B, the surface of the infected skin turned grey and rough. Histological examination showed complete epidermis in the control group, with neatly arranged epidermal cells, and clear reticular ridges and dermal papillae. In addition, the healthy dermis had wavy fibrous tissues, intact and uniformly distributed hair follicles and sebaceous glands, and no infiltration of inflammatory cells (shown in Figure 1C). In contrast, the S. aureus-infected skin displayed significantly thickened granular and spinous layers, uneven epidermal thickness, discontinuity between the epidermis and dermis, aberrant hair follicles and sebaceous glands, and massive neutrophil infiltration into the dermal and subcutaneous tissues (shown in Figure 1D). The thickness of the stratum spinosum and epidermis were measured, and were significantly greater in the infected versus the control skin (P-value <0.001, shown in Figure 1E and F).

Figure 1

Histological evaluation of the skin (n = 5 mice for each group). (A and B), representative images of normal untreated control mouse (A) and the treated mouse with injection of S. aureus (B). (C and D), representative images of the hematoxylin-eosin staining performed to observe the morphological structure of skin in the control group (C) and treated group (model) (D). (E and F) were thicknesses quantifications of stratum spinosum and stratum epidermis from (C) and (D) (n=5). P-value from t-test. *P <0.05. Error bars indicate mean ± SD.

Identification of DEGs

A total of 638 genes were identified in the inflamed skin samples relative to normal skin (|log2 fold change|> 3 and P-value <0.01), of which 324 were up-regulated and 314 were down-regulated (shown ). As shown in Table 1, 38 of the top 100 upregulated genes encode immune response-related proteins, of which five (C7, C5ar1, F13a1, C6 and Cfh) are involved in complement and coagulation cascades, six (Serpina3n, Saa2, Saa1, Saa3, Lbp and Cd163) are related to the acute-phase response, five (Saa1, Saa2, Saa3, Cxcl13 and Ccl6) regulate chemotaxis, and four (Oas1a, Oas1f, Oas1g and Tlr8) are involved in virus-specific immune response. Twenty-one genes with the highest increase in transcript levels encode proteins that have known functions in metabolism and cell proliferation, whereas 8 are functionally uncharacterized. In addition, 14 upregulated genes encode structural proteins of the dermis, including tenascin N (Tnn), osteomodulin (Omd), cartilage intermediate layer protein (Clip), myocilin (Myoc), and multimerin 1 (Mmrn1). Three metallopeptidase-encoding genes (Mmp8, Mmp11 and Lvrn) and seven genes encoding redox proteins also showed significant increase in their transcript levels.

Table 1

Top Genes with Significantly Increased Transcript Levels

Symbol	Gene Name	Log2FC	Function
Oas1a	2’-5’ oligoadenylate synthetase 1A	7.8666	Immune function
Oas1f	2’-5’ oligoadenylate synthetase 1F	7.7425	Immune function
Saa2	Serum amyloid A 2	7.5192	Immune function
Saa1	Serum amyloid A 1	7.03	Immune function
Oas1g	2’-5’ oligoadenylate synthetase 1G	6.5824	Immune function
C7	Complement component 7	6.4889	Immune function
Saa3	Serum amyloid A 3	6.2966	Immune function
Serpina3k	Serine (or cysteine) peptidase inhibitor, clade A, member 3K	6.0593	Immune function
Ankrd1	Ankyrin repeat domain 1 (cardiac muscle)	6.0418	Immune function
Serpina3m	Serine (or cysteine) peptidase inhibitor, clade A, member 3M	5.8389	Immune function
Cxcl13	Chemokine (C-X-C motif) ligand 13	5.6559	Immune function
Ms4a8a	Membrane-spanning 4-domains, subfamily A, member 8A	5.5616	Immune function
Fcer1a	Fc receptor, IgE, high affinity I, alpha polypeptide	5.2366	Immune function
Ctla2a	Cytotoxic T lymphocyte-associated protein 2 alpha	5.1228	Immune function
Vsig4	V-set and immunoglobulin domain containing 4	5.0285	Immune function
Fcgr4	Fc receptor, IgG, low affinity IV	4.862	Immune function
Ccl6	Chemokine (C-C motif) ligand 6	4.8081	Immune function
Serpina3n	Serine (or cysteine) peptidase inhibitor, clade A, member 3N	4.7055	Immune function
Cpa3	Carboxypeptidase A3, mast cell	4.5542	Immune function
Musk	Muscle, skeletal, receptor tyrosine kinase	4.4989	Immune function
Cfh	Complement component factor h	4.3727	Immune function
F13a1	Coagulation factor XIII, A1 subunit	4.3551	Immune function
Ctxn3	Cortexin 3	4.3057	Immune function
Ms4a2	Membrane-spanning 4-domains, subfamily A, member 2	4.1518	Immune function
Itgam	Integrin alpha M	4.145	Immune function
Lyve1	Lymphatic vessel endothelial hyaluronan receptor 1	4.1333	Immune function
Vnn1	Vanin 1	4.0442	Immune function
Il1f6	Interleukin 1 family, member 6	4.0035	Immune function
C5ar1	Complement component 5a receptor 1	3.9952	Immune function
Cd5	CD5 antigen	3.9773	Immune function
Adgre1	Adhesion G protein-coupled receptor E1	3.9568	Immune function
Igkc	Immunoglobulin kappa constant	3.945	Immune function
Il1rl1	Interleukin 1 receptor-like 1	3.9318	Immune function
Tlr8	Toll-like receptor 8	3.9235	Immune function
Cd163	CD163 antigen	3.8919	Immune function
C6	Complement component 6	3.8262	Immune function
Lbp	Lipopolysaccharide binding protein	3.8179	Immune function
Cd209e	CD209e antigen	3.7556	Immune function
Glp1r	Glucagon-like peptide 1 receptor	7.1474	Metabolism/Cell proliferation/Regulation
Capn11	Calpain 11	5.8636	Metabolism/Cell proliferation/Regulation
Pdk4	Pyruvate dehydrogenase kinase, isoenzyme 4	5.6197	Metabolism/Cell proliferation/Regulation
Calca	Calcitonin/calcitonin-related polypeptide, alpha	5.5922	Metabolism/Cell proliferation/Regulation
Inmt	Indolethylamine N-methyltransferase	5.5074	Metabolism/Cell proliferation/Regulation
Trim63	Tripartite motif-containing 63	5.4875	Metabolism/Cell proliferation/Regulation
Ankrd2	Ankyrin repeat domain 2 (stretch responsive muscle)	4.6329	Metabolism/Cell proliferation/Regulation
4930486L24Rik	RIKEN cDNA 4930486L24 gene	4.5061	Metabolism/Cell proliferation/Regulation
Capn9	Calpain 9	4.3767	Metabolism/Cell proliferation/Regulation
Mustn1	Musculoskeletal, embryonic nuclear protein 1	4.3373	Metabolism/Cell proliferation/Regulation
Ubtd1	Ubiquitin domain containing 1	4.2869	Metabolism/Cell proliferation/Regulation
Pamr1	Peptidase domain containing associated with muscle regeneration 1	4.2856	Metabolism/Cell proliferation/Regulation
Galnt15	UDP-N-acetyl-alpha-D-galactosamine:polypeptide N-acetylgalactosaminyltransferase 15	4.2751	Metabolism/Cell proliferation/Regulation
Ifi207	Interferon activated gene 207	3.9666	Metabolism/Cell proliferation/Regulation
Fbxo32	F-box protein 32	3.9597	Metabolism/Cell proliferation/Regulation
Retnlg	Resistin like gamma	3.9338	Metabolism/Cell proliferation/Regulation
Clec10a	C-type lectin domain family 10, member A	3.9336	Metabolism/Cell proliferation/Regulation
Plac8	Placenta-specific 8	3.9177	Metabolism/Cell proliferation/Regulation
Doc2b	Double C2, beta	3.7793	Metabolism/Cell proliferation/Regulation
Ifi205	Interferon activated gene 205	3.7661	Metabolism/Cell proliferation/Regulation
Fgf23	Fibroblast growth factor 23	3.7351	Metabolism/Cell proliferation/Regulation
Prg4	Proteoglycan 4 (megakaryocyte stimulating factor, articular superficial zone protein)	7.8227	Structure dermal components
C1qtnf3	C1q and tumor necrosis factor related protein 3	4.5314	Structure dermal components
Gdf5	Growth differentiation factor 5	4.4774	Structure dermal components
Fcna	Ficolin A	4.1605	Structure dermal components
Cilp	Cartilage intermediate layer protein, nucleotide pyrophosphohydrolase	4.0797	Structure dermal components
Reln	Reelin	4.0422	Structure dermal components
Mmp11	Matrix metallopeptidase 11	4.0318	Structure dermal components
Myoc	Myocilin	4.0125	Structure dermal components
Wisp2	WNT1 inducible signaling pathway protein 2	3.9736	Structure dermal components
Lvrn	Laeverin	3.9249	Structure dermal components
Omd	Osteomodulin	3.8767	Structure dermal components
Mmrn1	Multimerin 1	3.819	Structure dermal components
Mmp8	Matrix metallopeptidase 8	3.7604	Structure dermal components
Tnn	Tenascin N	3.7583	Structure dermal components
Mrgprg	MAS-related GPR, member G	7.451	Transportation
Mrgprb3	MAS-related GPR, member B3	4.9451	Transportation
Reg1	Regenerating islet-derived 1	4.3582	Transportation
Mrgprb1	MAS-related GPR, member B1	4.1727	Transportation
Adap2	ArfGAP with dual PH domains 2	4.0426	Transportation
Best2	Bestrophin 2	3.9707	Transportation
Hspb7	Heat shock protein family, member 7 (cardiovascular)	3.9061	Transportation
Slc28a2	Solute carrier family 28 (sodium-coupled nucleoside transporter), member 2	3.8161	Transportation
Apoc3	Apolipoprotein C-III	4.4993	Lipid metabolism
Pcolce2	Procollagen C-endopeptidase enhancer 2	3.7774	Lipid metabolism
Abca6	ATP-binding cassette, sub-family A (ABC1), member 6	3.7789	Lipid metabolism
Gldc	Glycine decarboxylase	4.8255	Oxidation-reduction process
Nox4	NADPH oxidase 4	4.6939	Oxidation-reduction process
Cyp8b1	Cytochrome P450, family 8, subfamily b, polypeptide 1	4.6826	Oxidation-reduction process
Gpx3	Glutathione peroxidase 3	4.0838	Oxidation-reduction process
Nat8f3	N-acetyltransferase 8 (GCN5-related) family member 3	3.9506	Oxidation-reduction process
Lox	Lysyl oxidase	3.941	Oxidation-reduction process
Tdo2	Tryptophan 2,3-dioxygenase	3.7899	Oxidation-reduction process
Serpinb6e	Serine (or cysteine) peptidase inhibitor, clade B, member 6e	6.5088	Function unknown
4932443I19Rik	RIKEN cDNA 4932443I19 gene	4.583	Function unknown
Fam124a	Family with sequence similarity 124, member A	4.438	Function unknown
Mir3966	MicroRNA 3966	4.2942	Function unknown
Fam180a	Family with sequence similarity 180, member A	4.0238	Function unknown
Klhl38	Kelch-like 38	3.8075	Function unknown
Sirpb1c	Signal-regulatory protein beta 1C	3.7748	Function unknown
Prom1	Prominin 1	3.7664	Function unknown

The top 100 downregulated genes are listed in Table 2, of which 24 have unknown function. Only two genes encode immune response-related proteins, including PDZ binding kinase (Pbk) and transient receptor potential cation channel (Trpm1). Thirty-two genes regulate metabolism and cell proliferation, of which seven encode cell cycle proteins (Ccnb1, Cdk1, Ccnb2, Plk1, Cdc20, Cdc25c and Ccna2), and four encode negative regulators of cell proliferation (Hist1h2ag, Hist1h2ai, Nppc and Hist1h2an). The proteins encoded by Hist1h2ag and Hist1h2ai also play important roles in DNA repair. In addition, genes encoding DNA damage response proteins, including Exo1, Uhrf1, Ticrr, Kif22, Pclaf, Tgm6 and Pif1, were also significantly downregulated by S. aureus infection. Furthermore, 16 genes encode dermal proteins, of which 14 genes encode proteins related to keratin or keratin filament. Other major downregulated genes were those encoding proteins involved lipid metabolism (5) and transportation (10).

Table 2

Top Genes with Significantly Decreased Transcript Levels

Symbol	Gene Name	Log2FC	Function
Pbk	PDZ binding kinase	−4.10114	Immune Function
Trpm1	Transient receptor potential cation channel, subfamily M, member 1	−3.85416	Immune Function
Sct	Secretin	−7.16723	Metabolism/Cell proliferation/Regulation
Gprc5d	G protein-coupled receptor, family C, group 5, member D	−5.06897	Metabolism/Cell proliferation/Regulation
Spink12	Serine peptidase inhibitor, Kazal type 12	−4.9971	Metabolism/Cell proliferation/Regulation
Adh6b	Alcohol dehydrogenase 6B (class V)	−4.74925	Metabolism/Cell proliferation/Regulation
Plk1	Polo-like kinase 1	−4.74748	Metabolism/Cell proliferation/Regulation
Avp	Arginine vasopressin	−4.74139	Metabolism/Cell proliferation/Regulation
Ky	Kyphoscoliosis peptidase	−4.70848	Metabolism/Cell proliferation/Regulation
Mylk4	Myosin light chain kinase family, member 4	−4.67611	Metabolism/Cell proliferation/Regulation
Nkx1-2	NK1 transcription factor related, locus 2 (Drosophila)	−4.5013	Metabolism/Cell proliferation/Regulation
Ube2c	Ubiquitin-conjugating enzyme E2C	−4.40454	Metabolism/Cell proliferation/Regulation
Cdk1	Cyclin-dependent kinase 1	−4.39539	Metabolism/Cell proliferation/Regulation
Scrg1	Scrapie responsive gene 1	−4.34082	Metabolism/Cell proliferation/Regulation
Ccnb2	Cyclin B2	−4.22004	Metabolism/Cell proliferation/Regulation
Mxd3	Max dimerization protein 3	−4.189	Metabolism/Cell proliferation/Regulation
Cdc20	Cell division cycle 20	−4.17061	Metabolism/Cell proliferation/Regulation
Cdca5	Cell division cycle associated 5	−4.15544	Metabolism/Cell proliferation/Regulation
Cdc25c	Cell division cycle 25C	−4.15077	Metabolism/Cell proliferation/Regulation
Dlx4	Distal-less homeobox 4	−4.11843	Metabolism/Cell proliferation/Regulation
Tagln3	Transgelin 3	−4.11208	Metabolism/Cell proliferation/Regulation
Ska3	Spindle and kinetochore associated complex subunit 3	−4.09096	Metabolism/Cell proliferation/Regulation
Tgm6	Transglutaminase 6	−4.06916	Metabolism/Cell proliferation/Regulation
Ccnb1	Cyclin B1	−4.04623	Metabolism/Cell proliferation/Regulation
Ccna2	Cyclin A2	−4.04422	Metabolism/Cell proliferation/Regulation
Cdca3	Cell division cycle associated 3	−3.98818	Metabolism/Cell proliferation/Regulation
Hist1h2an	Histone cluster 1, H2an	−3.97682	Metabolism/Cell proliferation/Regulation
Hoxc13	Homeobox C13	−3.92584	Metabolism/Cell proliferation/Regulation
Nppc	Natriuretic peptide type C	−3.8806	Metabolism/Cell proliferation/Regulation
Cdca8	Cell division cycle associated 8	−3.87644	Metabolism/Cell proliferation/Regulation
Kif11	Kinesin family member 11	−3.87194	Metabolism/Cell proliferation/Regulation
E2f8	E2F transcription factor 8	−3.85788	Metabolism/Cell proliferation/Regulation
Mybl2	Myeloblastosis oncogene-like 2	−3.83805	Metabolism/Cell proliferation/Regulation
Jakmip2	Janus kinase and microtubule interacting protein 2	−3.81494	Metabolism/Cell proliferation/Regulation
Ncapg	Non-SMC condensin I complex, subunit G	−3.80713	Metabolism/Cell proliferation/Regulation
Krtap8-1	Keratin associated protein 8-1	−7.57011	Structure dermal components
Krtap11-1	Keratin associated protein 11-1	−6.15673	Structure dermal components
Fgf5	Fibroblast growth factor 5	−5.55214	Structure dermal components
Krt31	Keratin 31	−5.53606	Structure dermal components
Krt90	Keratin 90	−5.40147	Structure dermal components
Krt34	Keratin 34	−5.24946	Structure dermal components
Krt36	Keratin 36	−5.19916	Structure dermal components
Krtap13-1	Keratin associated protein 13-1	−5.18137	Structure dermal components
Krtap16-1	Keratin associated protein 16-1	−4.70055	Structure dermal components
Krt84	Keratin 84	−4.57811	Structure dermal components
Krt33a	Keratin 33A	−4.33963	Structure dermal components
Krtap24-1	Keratin associated protein 24-1	−4.15059	Structure dermal components
Krtap2-4	Keratin associated protein 2-4	−3.97812	Structure dermal components
Krt81	Keratin 81	−3.96064	Structure dermal components
Krt83	Keratin 83	−3.91786	Structure dermal components
Mmp12	Matrix metallopeptidase 12	−3.8559	Structure dermal components
Idi2	Isopentenyl-diphosphate delta isomerase 2	−7.48319	Lipid metabolism
Htr5b	5-hydroxytryptamine (serotonin) receptor 5B	−5.06777	Lipid metabolism
Gykl1	Glycerol kinase-like 1	−4.96037	Lipid metabolism
Birc5	Baculoviral IAP repeat-containing 5	−4.03313	Lipid metabolism
Melk	Maternal embryonic leucine zipper kinase	−3.92305	Lipid metabolism
Hist1h2ag	Histone cluster 1, H2ag	−5.21662	DNA repair
Pclaf	PCNA clamp associated factor	−4.22398	DNA repair
Tgm6	Transglutaminase 6	−4.06916	DNA repair
Hist1h3g	Histone cluster 1, H3g	−4.02951	DNA repair
Exo1	Exonuclease 1	−4.01017	DNA repair
Kif22	Kinesin family member 22	−3.9821	DNA repair
Pif1	PIF1 5’-to-3’ DNA helicase	−3.91843	DNA repair
Ticrr	TOPBP1-interacting checkpoint and replication regulator	−3.86324	DNA repair
Hist1h2ai	Histone cluster 1, H2ai	−3.86315	DNA repair
Uhrf1	Ubiquitin-like, containing PHD and RING finger domains, 1	−3.83378	DNA repair
Cacng4	Calcium channel, voltage-dependent, gamma subunit 4	−5.02411	Transportation
Gja3	Gap junction protein, alpha 3	−4.46609	Transportation
Slc13a5	Solute carrier family 13 (sodium-dependent citrate transporter), member 5	−4.42664	Transportation
Mup12	Major urinary protein 12	−4.31564	Transportation
Oca2	Oculocutaneous albinism II	−4.28906	Transportation
Slc5a5	Solute carrier family 5 (sodium iodide symporter), member 5	−4.14539	Transportation
Mup19	Major urinary protein 19	−4.10512	Transportation
Kif20a	Kinesin family member 20A	−3.97648	Transportation
Slc24a5	Solute carrier family 24, member 5	−3.91393	Transportation
Cnih2	Cornichon family AMPA receptor auxiliary protein 2	−3.8604	Transportation
Wfdc18	WAP four-disulfide core domain 18	−6.27637	Function unknown
Tchhl1	Trichohyalin-like 1	−5.66721	Function unknown
AC154502.1		−5.16798	Function unknown
Otop2	Otopetrin 2	−5.15936	Function unknown
AC125444.4		−4.61759	Function unknown
Pinlyp	Phospholipase A2 inhibitor and LY6/PLAUR domain containing	−4.44108	Function unknown
Slc45a2	Solute carrier family 45, member 2	−4.28805	Function unknown
Kif18b	Kinesin family member 18B	−4.24566	Function unknown
BC030867	cDNA sequence BC030867	−4.19876	Function unknown
Sptssb	Serine palmitoyltransferase, small subunit B	−4.14632	Function unknown
Wdr95	WD40 repeat domain 95	−4.142	Function unknown
Fer1l4	Fer-1-like 4 (C. elegans)	−4.14156	Function unknown
Lrr1	Leucine rich repeat protein 1	−4.08969	Function unknown
B3galt5	UDP-Gal:betaGlcNAc beta 1,3-galactosyltransferase, polypeptide 5	−4.08608	Function unknown
Ska1	Spindle and kinetochore associated complex subunit 1	−4.06129	Function unknown
Tmem28	Transmembrane protein 28	−4.03851	Function unknown
Adamts20	A disintegrin-like and metallopeptidase (reprolysin type) with thrombospondin type 1 motif, 20	−3.9805	Function unknown
Spc24	SPC24, NDC80 kinetochore complex component, homolog (S. cerevisiae)	−3.96965	Function unknown
Fam64a	Family with sequence similarity 64, member A	−3.96887	Function unknown
Sapcd2	Suppressor APC domain containing 2	−3.95029	Function unknown
Spag5	Sperm associated antigen 5	−3.93997	Function unknown
AC108401.1		−3.90024	Function unknown
Tmem61	Transmembrane protein 61	−3.88278	Function unknown
2700016F22Rik	RIKEN cDNA 2700016F22 gene	−3.81196	Function unknown

GO and KEGG Pathway Analysis

The top 20 GO terms of up- and down-regulated DEGs are shown in Figure 2. The top 5 GO terms for the up-regulated genes were regulation of response to stimulus, extracellular space, defense response, regulation of response to external stimulus and immune system process (shown in Figure 2A), and that for the down-regulated DEGs were keratin filament, chromosome segregation, mitotic cell cycle process, cell cycle process and intermediate filament (shown in Figure 2B). On the other hand, the most enriched KEGG pathways of the up-regulated genes were complement and coagulation cascades and Staphylococcus aureus infection (shown in Figure 2C), and that of the down-regulated genes were cell cycle and oocyte meiosis (shown in Figure 2D). GO clustering enrichment analysis further showed that the up-regulated DEGs were mainly involved in acute-phase response (enrichment score ES = 3.06, gene count GC = 8), chemotaxis (ES = 2.86, GC = 18), negative regulation of inflammation factor production (ES = 2.12, GC=6) and lipid metabolism (ES = 2, GC = 10) (shown in Figure 3A). The GO term of signal transduction was enriched in most up-regulated genes (39). As shown in Figure 3B, the down-regulated DEGs were associated with cell cycle (ES = 39.15, GC = 70), chromosome segregation (ES = 20.58, GC = 59), structural molecule (ES = 9.10, GC = 22), microtubule-based movement (ES = 8.72, GC = 32), melanin biosynthetic process (ES = 4.42, GC = 7), and DNA metabolism (ES = 2.21, GC = 8). The detailed information is presented in and .

Figure 2

The top 20 GO terms (A) and KEGG pathways (B) of different expressed genes enrichment in RNA-sequencing (P-value ≤ 0.05). Gene Counts: number of target genes in each GO term or KEGG pathway. Rich factor: the ratio of the number of target genes divided by the number of all the genes in each GO term or KEGG pathway. The size of pot indicated the gene counts, and the color reflected the different P-value.

Figure 3

GO clustering results of DEGs with a criterion of Enrichment Score > 1. (A) Up-regulated. (B) Down-regulated. Counts represent the gene number enriched in each term. The Enrichment score was obtained according to the built-in program in DAVID, which ranks the biological significance of gene groups based on overall P-value of all enriched annotation terms. The higher the enrichment score is, the more important the term is.

PPI Network of DEGs and Core Genes in the PPI Network

The PPI network contained 1498 nodes (DEGs) and 6701 edges (interactions between the DEGs). Using the twofold median value of node degree in this network as a cutoff point,23 1104 nodes were identified as hubs. The hub interaction network was constructed (), and topological features, such as “Average Shortest Path Length”, “Closeness Centrality”, “Betweenness Centrality”, “Neighborhood Connectivity”, and “Degree” were calculated for each hub (shown in ). Clustering analysis further revealed 11 groups based on the physiological processes (shown in ), and the hub genes in each cluster were mainly associated with cell growth and death, cellular community, cell cycle, organismal systems, immune system, infectious diseases, signal transduction and apoptosis (shown in Table 3). A total of 295 major hubs (degree higher than 2-fold of all hubs; median value = 6) were screened and 4 clusters were built (shown in Figure 4). As shown in Table 4, the top 5 pathways in Cluster 1 were cell cycle (P- value = 1.17E-30), p53 signaling pathway (P-value = 1.76E-12), DNA replication (P-value = 7.53E-09), viral carcinogenesis (P-value = 3.20E-07), and homologous recombination (P-value = 2.31E-06), that in Cluster 2 were ErbB signaling pathway (P-value = 1.57E-18), Natural killer cell mediated cytotoxicity (P-value = 4.44E-17), PI3K-Akt signaling pathway (P-value = 2.55E-16), pathways in cancer (P-value = 8.16E-16) and proteoglycans in cancer (P-value = 3.22E-15), in Cluster 3 were pathways in cancer (P-value = 1.29E-27), chronic myeloid leukemia (P-value = 1.52E-21), hepatitis B (P-value = 1.13E-20), pancreatic cancer (P-value = 1.40E-19), and prostate cancer (P-value = 1.57E-19), and in Cluster 4 were oocyte meiosis (P-value = 0.005), legionellosis (P-value = 0.05), and salmonella infection (P-value = 0.08). Furthermore, the top 4 high-degree hub nodes in each cluster included CDK1, MCM10, CD5, PDGFRA, BRCA1, LEF1, AURKA, BUB1 etc. (shown in Table 5), and were the most down-regulated genes.

Table 3

Main GO and Pathway Enrichment Analysis Relating to All Hubs

Clusters in Network	GO and Pathway Enrichment Analysis
1	Cell growth and death
2	Cellular community
3	Infectious diseases and cancers
4	Cell cycle
5	Organismal Systems
6	Immune system
7	Infectious diseases
8	Signal transduction
9	Apoptosis
10	——
11	——

Abbreviation: GO, Gene ontology.

Figure 4

The top clusters of model with relatively high degrees selected from the protein–protein interaction network. The hubs were selected following the criteria of degree higher than 2-fold of all hubs, and median value=6. The clusters were built using Cytoscape ver. 3.5.1.

Table 4

Pathway Enrichment Analysis of PPI Clusters Related to Major Hubs (Top 5)

Clusters	Pathway Enrichment Analysis (KEGG)‡	Count†	P-value
1	hsa04110: Cell cycle	26	1.17E-30
	hsa04115: p53 signaling pathway	12	1.76E-12
	hsa03030: DNA replication	8	7.53E-09
	hsa05203: Viral carcinogenesis	12	3.20E-07
	hsa03440: Homologous recombination	6	2.31E-06
2	hsa04012: ErbB signaling pathway	19	1.57E-18
	hsa04650: Natural killer cell mediated cytotoxicity	20	4.44E-17
	hsa04151: PI3K-Akt signaling pathway	28	2.55E-16
	hsa05200: Pathways in cancer	29	8.16E-16
	hsa05205: Proteoglycans in cancer	22	3.22E-15
3	hsa05200: Pathways in cancer	46	1.29E-27
	hsa05220: Chronic myeloid leukemia	22	1.52E-21
	hsa05161: Hepatitis B	27	1.13E-20
	hsa05212: Pancreatic cancer	20	1.40E-19
	hsa05215: Prostate cancer	22	1.57E-19
4	hsa04114: Oocyte meiosis	3	0.00514
	hsa05134: Legionellosis	2	0.053695
	hsa05132: Salmonella infection	2	0.081498

Note: †Count represents the gene number.

Abbreviation: ‡KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 5

The Top 5 DEGs with Higher Degrees in the Clusters Related to Major Hubs

Clusters	Gene	Expression Changes	Degree	Log2FC
1	Cdk1	Down	16	−4.40
	Mcm10	Down	14	−3.78
	Cdc6	Down	12	−3.01
	Plk1	Down	10	−4.75
	Mcm5	Down	10	−3.08
2	Cd5	Up	12	3.98
	Pdgfra	Down	11	3.36
	Cdc25c	Down	7	−4.15
	Agtr1	Up	6	3.62
	Racgap1	Down	5	−3.28
3	Brca1	Down	35	−3.10
	Lef1	Down	12	−3.70
	Igfbp3	Up	10	3.03
	Top2a	Down	9	−3.58
	Msx1	Down	9	−3.50
4	Aurka	Down	2	−3.70
	Bub1	Down	2	−3.80
	Myoc	Up	2	4.01
	Cdc20	Down	2	−4.17
	Ndc80	Down	1	−3.59

qRT-PCR Confirmation of Differential mRNA Expression

Expression of Nox4, Mmrn1, Mcm5, Msx1, and Fgf5 was further validated by qRT-PCR and were coincident with RNA-Seq data (shown in Figure 5). Thus, both approaches showed a strong correlation.

Figure 5

Relative expression of each selected gene validated by qRT-PCR using GAPDH as the housekeeping gene. The expression of each gene relative to GAPDH was calculated using a modified method of delta Ct(2−ΔCt). Data represent the mean ± SD of triplicate experiments. ***P value <0.001. (A) The relative expression of Nox4. (B) The relative expression of Mmrn1. (C) The relative expression of Msx1. (D) The relative expression of Fgf5. (E) The relative expression of Mcm5. (F) Log2FC values of the selected genes from the RNA-Seq data.

Discussion

The Gram-positive bacterium S. aureus is the most common cause of skin infections in humans,15 although the molecular mechanisms underlying its pathological process are largely unknown. In this study, we established a mouse model of intradermal S. aureus infection, and observed the cardinal signs of localized inflammation, including epidermal thickening and massive infiltration of neutrophils in the dermis and epidermis. Miller et al summarized the host defense against S. aureus cutaneous infections and pointed again that one hallmark of S. aureus infections is the formation of neutrophil abscess, which is necessary for bacterial clearance.24 The recruitment of neutrophils to a site of S. aureus infection in the skin involves two steps, one is the innate immune process, recognition pathogen-associated molecular patterns of invading microorganisms, and the other is the subsequent initiation of pro-inflammatory immune responses, such as the secretion or production of cytokines, chemokines and adhesion molecules which promote neutrophil rolling, adhesion and diapedesis.24,25 While the inflammatory response combats pathogen invasion, uncontrolled infiltration of granulocytes can cause extensive tissue damage and epidermal dysfunction.26–28 In fact, several regulatory mechanisms limit the influx of neutrophils and prevent collateral damage due to a persistent pro-inflammatory state.29

To reveal the gene expression changes behind the above observations we carried out a transcriptome analysis. A total of 638 DEGs were identified in the study, of which 324 and 314 genes were, respectively, up- and down-regulated in the inflamed skin relative to the normal skin. The up-regulated DEGs were mainly enriched in GO terms like response to stimulus, immune system process and signal transduction, and pathways related to complement and coagulation cascades and S. aureus infection. Furthermore, 38 genes with the highest increase in transcript levels encoded immune response-related proteins like cytokines and chemokines, and the complement and coagulation cascades were enriched in most upregulated genes. Consistent with our findings, Brady et al12 reported 39 highly upregulated genes associated with the immune system in a mouse model of S. aureus soft tissue infection. A major part of these 39 genes (17 genes) encoding proteins included cytokines, chemokines, and chemokine-like proteins, as well as Cxcr1. Three genes encoding cytolytic granule proteins were included according to Brady’s study.12 While, we did not find related genes encoding cytolytic granule proteins in our study.

The complement system is a protein hydrolysis cascade that opsonizes pathogens and augments the effect of antibodies and phagocytic cells.30 Thirteen complement-related genes were differentially expressed in the inflamed skin, of which C7, C5ar1, F13a1, C6 and Cfh were highly upregulated. In addition, the top 4 enriched GO clusters following S. aureus infection were acute-phase response, chemotaxis, negative regulation of inflammation factor production and lipid metabolism.

Pathogenic infection also triggers the release of pro-inflammatory cytokines, chemokines and adhesion molecules by keratinocytes and the resident skin immune cells (such as Langerhans cells and γδ T cells in the epidermis, and dendritic cells, macrophages, fibroblasts, mast cells, B and T cells, plasma cells and natural killer cells in the dermis).25,28,29,31 These molecules promote the neutrophils recruitment which plays an important role in bacterial clearance.

In addition, several cytokines that regulate keratinocyte proliferation and differentiation are often dysregulated during inflammation.32 The chemokines in particular play an important part in the recruitment, activation and differentiation of inflammatory cells.33 For instance, Boeckmans et al detected increased expression of inflammatory cytokines such as CCL2, CCL5, CCL7, CCL8, CXCL5, CXCL8, IL1a, IL6 and IL11 in a human skin stem cell-derived model of non-alcoholic steatohepatitis (NASH).34 Consistent with this, 38 immune response-related genes (cytokines and chemokines, and the complement and coagulation cascades) above mentioned were upregulated in the S. aureus-infected skin, Cxcl13, Cxcl5, Ccl6 and Ccl9, for instance. There still existed other genes associated with neutrophil recruitment showing an increase in the expression. Several adhesion molecules, such as cell adhesion molecule3 (Cadm3), selectins (Sell and Selp) and adhesion molecule interacting with CXADR antigen 1(Amica1), playing an important part in neutrophil rolling, adhesion and diapedesis, increased in the inflamed skin.

The production and secretion of pro-inflammatory mediators depends on the activation of sensors like the toll-like receptors (TLRs), interleukin-1 receptor (IL-1R) that recognize pathogen-associated molecular patterns (PAMPs), and trigger downstream signaling pathways. The TLR or IL-1R pathway is associated with the sensitivity of various tissues and organs to S. aureus infection. The defects in these pathways can result in an increased susceptibility to the infections.35,36 The TLR or IL-1R family members signal through myeloid differentiation primary response protein 88 (MYD88), which subsequently interacts with IL-1R-associated kinase 4 (IRAK4) to activate TNFR-associated factor 6(TRAF6).24 In this study, both TLR and IL-1R family members were detected at high transcript levels. The TLR pathway members, such as toll-like receptor 4(Tlr4), toll-like receptor 7(Tlr7), toll-like receptor 8(Tlr8), cathepsin K (Ctsk), and lipopolysaccharide binding protein (Lbp), were up-regulated in the inflamed skin. It has been shown that TLRs are activated by recognizing S. aureus lipopeptides and lipoteichoic acid.37 These genes subsequently participate in the regulation of interleukin biosynthetic process, together with the IL-1R family members like Il1r2 and Il1rl1.

Miller et al reported that the activation of IL-1R is due to the stimulation of IL-1α and IL-1β.24 They found that different establishment of model results in differences in immune response. Both IL-1α and IL-1β contribute to host defense during a superficial S. aureus skin infection, whereas IL-1β has a predominant role during a deeper S. aureus skin infection.38 In our study, the subcutaneous model was established. Meanwhile, we did not detect the increase of IL-1β, which might be due to the differences in mice and the strains we used between the literature38 and our study.

Inflammasomes are molecular platforms activated in response to cell infection or stress that trigger the maturation of proinflammatory cytokines to participate in innate immune defense.39 Miller et al reported that inflammasome-mediated production of IL-1β is required for neutrophil recruitment against S. aureus in vivo.40 There are several types of inflammasome, such as NLRP1, NLRP3, NLRC4, and AIM2. In our study, we detected the gene expressions participating in the inflammasome complex assembly, although these genes were not included as DEGs by the subjective criteria used. The increased expression of Nlrp3 (NLR family, pyrin domain containing 3, Log2FC = 1.57, P-Adjust=0.001), Gsdmd (gasdermin D, Log2FC = 1.30, P-Adjust=5.01E-14), and Cd36 (CD36 antigen, Log2FC = 1.59, P-Adjust=1.48E-21) were detected, which all participated in assembly of NLRP3 inflammasome, and the secretion of IL-1β. Nlrc4 (NLR family, CARD domain containing 4, Log2FC = 1.15, P-Adjust=0.021) participating in the IL-1β secretion, also played an important part in the forming of NLRC4 inflammasome complex. In the study, Casp12 (caspase 12, Log2FC = 1.58, P-Adjust=3.59E-20) was increased, taking part in three types of inflammasome complex (NLRP3, AIM2 and NLRC4).

People believed that TH17 cells have a role in host defense against S. aureus at cutaneous site of infection.24,41 TH17 cells produce IL-17, IL-17F, and IL-22, thereby inducing a massive tissue reaction owing to the broad distribution of the IL-17 and IL-22 receptors.41 In the study, genes participating in IL-17 signaling pathway were accelerated, including lipocalin 2(Lcn2), chemokine (C-X-C motif) ligand 5(Cxcl5), S100 calcium binding protein A8 (calgranulin A) (S100a8), and S100 calcium binding protein A9 (calgranulin B) (S100a8). Inflammatory responses involving IL-17 probably contribute to arthritis, asthma, skin immune reactions and autoimmune disorders. Meanwhile, interleukin 17 receptor A (Il17ra), chemokine (C-C motif) ligand 7 (Ccl7), chemokine (C-C motif) ligand 12 (Ccl12) were also detected in a relatively high level, though they were not included as DEGs by the subjective criteria used. The activation of IL-17 receptors (IL-17RA-IL-17RC) on epidermal keratinocytes results in the production of antimicrobial peptides that have bacteriostatic and bactericidal activity against S. aureus.24 Antimicrobial peptides such as β-defensins and cathelicidins are induced production subsequently to come into play bacteriostatic or bactericidal activity against S. aureus.25 In the study, two genes participating in cathelicidins related process were changed in the expression. Neutrophilic granule protein (Ngp) was accelerated. Kininogen 2 (Kng2) was also increased, which plays a role in peptide secretion. While, cathelicidin antimicrobial peptide (Camp) decreased in the study due to the different period of inflammatory process.

In addition, several growth factors promote the accumulation and proliferation of inflammatory cells in infected tissues, and the fibroblast growth factor genes Fgf14 and Fgf23 were activated in the infected skin. Fgf14 is the receptor of TNF-related weak inducer of apoptosis (TWEAK), which triggers the release of pro-inflammatory factors like IL-6, IL-8, MMP-1, IFNγ-inducible protein 10 and prostaglandin E2 from fibroblasts and synoviocytes.42 Fgf23 is a direct stimulator of systemic inflammation, and promotes the secretion of IL-6 and C-reactive protein by liver cells.43 Prostaglandins (PG), including PG2, prostaglandin D2 (PGD2), prostaglandin F2α (PGF2α) and prostacyclin (PGI2), play a key role in the development of acute inflammation. Genes encoding for PGD2 synthase and PGI2 synthase were up-regulated in the model group, which is consistent with the high levels of PGD2 and PGI2 observed in inflamed tissues.44 The PGF receptor, a member of the prostanoid receptor subfamily, was also upregulated in the infected skin, which may be related to increased function of PGF2α.

The aryl hydrocarbon receptor (AhR) plays an important role in B cell maturation, and activation of macrophages, dendritic cells and neutrophils following lipopolysaccharide challenge or influenza virus infection. It is known to increase the expression levels of cytochrome P450 (Cyp), which is a key regulator of xenobiotic metabolism and the generation of intracellular reactive oxygen species (ROS).45 The S. aureus-infected tissues showed high transcript levels of AhR, and the Cyp2b23, Cyp2e1, Cyp4a10, Cyp4b1 and Cyp8b1 genes. ROS act as second messengers in redox-sensitive signal transduction, inflammation and apoptosis.46 Furthermore, several antioxidants maintain intracellular redox homeostasis by neutralizing end products of oxidation, and mitigate tissue damage due to oxidative stress. Genes encoding glutathione peroxidases (Gpx3 and Gpx7), glutathione S-transferases (Gsta3, Gstp2, Gstt1 and Gstt2) and other antioxidants (for eg, Nox4, Gpx3 and Lox) were significantly upregulated in the model group.

Consistent with the morphological changes seen in the epidermis following S. aureus infection, the expression levels of genes encoding extracellular matrix proteins were significantly altered. Fourteen of the most upregulated genes encode proteins related to the dermis structure. For instance, Cilp encodes cartilage intermediate layer protein that prevents fibrosis during pressure overload-induced fibrotic remodeling.47 In addition, Mmrn1 is involved in platelet adhesion and can bind to multiple proteins, which is related to the occurrence and development of many diseases.48 Four genes (Mmp2, Mmp8, Mmp11 and Lvrn) encoding different MMPs were significantly upregulated in the model group. The MMPs are a family of structurally related matrix-degrading enzymes that accelerate tissue disintegration during various pathological processes, including inflammation.49 They are normally expressed at low levels, and rapidly induced by various extracellular stimuli including growth factors, cytokines, oncoproteins, and UV and infrared radiations.49 MMP8 or collagenase-2 (neutrophil collagenase) can cleave collagens, elastins, fibronectins, gelatins and laminins in the extracellular matrix (ECM). MMP2 or gelatinase A, degrades denatured collagens or gelatins, stromelysins or proteoglycanases. Consistent with the high transcript levels of MMPs in the inflamed tissues, ECM components including the keratin family, keratin associated proteins, and cytoskeletal proteins were also downregulated in the model group.

S. aureus infection also increased transcript levels of filaggrin (Flg) and the filaggrin family member 2 (Flg2). Flg overexpression is associated with the disruption of keratin and other cytoskeletal elements, loss of desmosomal proteins from intercellular junctions, and cell cycle arrest.50 Tight junction proteins maintain the integrity of the epidermal barrier and prevent transepidermal water loss (TEWL).51,52 We found that genes encoding small proline-rich protein 2D (Sprr2d) and Keratin 78 (Krt78) were increased in the model group. Likewise, Brady et al found that Sprr2a1, Sprr2j-ps, Sprr3, Sprr2k and Krt6b, which have known functions in epidermal formation, were upregulated in S. aureus-infected tissues.30 Angiopoietins 1–4 (Angpt1-4) are growth factors belonging to the angiopoietin/TIE (tyrosine kinase with Ig and EGF homology domains) family that are involved in angiogenesis.53 In our study, we observed that Angpt1, Angpt2 and angiopoietin-like (Angptl) 1, 2, 4 and 7 transcripts were upregulated following S. aureus infection. Both Angpts and Angptls are activated by inflammatory mediators and other stimuli,54 and exacerbate inflammation-related tissue damage by increasing endothelial permeability and pericyte shedding from the basement membrane. These changes induce vascular leakage and promote immune or/and trans-endothelial migration of cancer cells.54

The down-regulated genes were mainly linked to GO terms like cell cycle, cell division and cytoskeletal parts, and the cell cycle (mmu:04110), oocyte meiosis (mmu:04114) and homologous recombination (mmu:03440) pathways. Two genes encoding immunity-related proteins were downregulated, of which Pbk promotes biosynthesis of pro-inflammatory factors like TNF-α, IL-1 and cyclooxygenase-2 via the p38 MAPK pathway.55 Most downregulated genes were enriched in cell cycle and chromosome segregation GO clusters, and 32 of the top 100 downregulated genes participate in the processes of metabolism and cell proliferation. These results indicate that genes related to cell growth and death are repressed by inflammation. The GO term structural molecule was the third most enriched among these DEGs, and related to the increased expression of MMPs. Fourteen of the most downregulated genes are related to keratin or keratin filaments, which modulate hair formation.30 Furthermore, keratinocytes also regulate the immune barrier function of the skin by secreting cytokines and chemokines. However, since we used the Balb/c nude mice, altered gene expression related to hair formation did not translate to any phenotypic changes. In addition, skin inflammation in the model group was accompanied by reduced expression of the junction protein claudin, which is consistent with the findings reported before.51,52

There are four enriched clusters in the PPI network, and most hub genes were down-regulated and mainly related to cell proliferation. CDK1, CD5, IGFBP3 and CDC20 were identified as the major hub proteins. CDK1 is a member of the cyclin-dependent kinase that regulates cell cycle progression. In the present study, Cdks (including Cdk1, Cdk2 and Cdk4), cyclins (Cyclin A2, Cyclin B1 and Cyclin B2) and their downstream transcription factors (E2f1, E2f2, and E2f3) were down-regulated, which can lead to cell cycle arrest. Other proliferation-related genes including Cdc20, Ccna2, Cdc6, Bub1, Bub1b, Ccnb1 and Ccnb2 were also downregulated. This is consistent with the observation that eukaryotic cells respond to DNA damage by inducing cell cycle arrest and DNA repair. However, in a previous study using USA300 strain-injected mice, cell cycle-related genes were not among the top 50 DEGs. In their study, genes related to immune function, signal transduction and epidermal formation account for a large proportion.12 Perhaps, it is because the experimental designs of the two studies are quite different, such as the selection of animals, the type of Staphylococcus aureus, the conditions of model establishment, and the way of obtaining samples, etc. On the other hand, it also shows the complexity of the organisms, for there are differences in occurrence time in the biological processes directly and indirectly related to the inflammatory process. As for our study, CD5 was up-regulated. It is reported that CD5 can be considered an immune checkpoint inhibitor used to potentiate anti-tumoral T cell responses or promote immunity against infections and modulate autoimmune reactions. Increased CD5 expression has been shown in anergized T cells after DC antigen targeting, whereas CD5 downregulation allows cytotoxic T cells to overcome tumor evasion and promotes cell death of autoreactive T cells and B-CLL cells.56 Furthermore, the pathogen-recognition receptors CD209d and CD209e that mediate pathogen endocytosis and lysosomal degradation57 were also upregulated. This is expected since the increased expression of pathogenic antigens is a major determinant of the immune system’s efficacy. The insulin-like growth factor binding protein (IGFBP) is a family of 6 proteins that bind to IGF and regulate insulin signaling. IGFBP3 induces apoptosis of retinal endothelial cells by stimulating PKA.58 While the Igfbp3 gene was up-regulated in the S. aureus-infected skin, Pka transcript levels were unaffected, indicating that protein IGFBP3 likely does not interact with PKA in the inflamed skin. IGFBP also regulates EMT,59 which is an essential cellular program in embryogenesis, tumor metastasis and wound healing,60 and is triggered by growth factors, inflammatory cytokines, chemokines, transcription factors and other proteins.61 These results further confirm that cell proliferation is inhibited by S. aureus-induced inflammation.

Based on its susceptibility to beta-lactams, S. aureus is commonly described as methicillin-susceptible S. aureus (MSSA) or methicillin-resistant S. aureus (MRSA). Moran et al found that S. aureus USA300 subtype (MRSA) is the cause of most community-acquired skin infections.62 More importantly, the impact of MSSA on infection still needs to be fully addressed. Although, we used the ATCC 25923 strain (one kind of MSSA) to induce model in the study. The transcriptomic changes in genes controlling the immune system, cell proliferation and dermal structure were more or less similar with that observed in USA300-induced soft tissue infection.12 The results were fully discussed above.

There are still limitations in our research. It would be better if we considered the genes related to S. aureus infection reported publicly when performing gene expression verification by qRT-PCR. Uninfected skin tissues from the same challenged mice should be added in the research to evaluate the local response to infection. Besides, inflammation is a process of change and development, and our research pays more attention to the genetic changes in the period of acute inflammation.

Nevertheless, our results still provide some theoretical support in the study of the pathogenesis of S. aureus and the screening of anti-inflammatory or anti- S. aureus drugs.